LORENE
et_bin_bhns_extr_excurv.C
1 /*
2  * Method of class Et_bin_bhns_extr to compute the extrinsic curvature tensor
3  * in the Kerr-Schild background metric or in the conformally flat one
4  * with extreme mass ratio
5  *
6  * (see file et_bin_bhns_extr.h for documentation).
7  *
8  */
9 
10 /*
11  * Copyright (c) 2004-2005 Keisuke Taniguchi
12  *
13  * This file is part of LORENE.
14  *
15  * LORENE is free software; you can redistribute it and/or modify
16  * it under the terms of the GNU General Public License version 2
17  * as published by the Free Software Foundation.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 char et_bin_bhns_extr_excurv_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_excurv.C,v 1.4 2014/10/13 08:52:55 j_novak Exp $" ;
31 
32 /*
33  * $Id: et_bin_bhns_extr_excurv.C,v 1.4 2014/10/13 08:52:55 j_novak Exp $
34  * $Log: et_bin_bhns_extr_excurv.C,v $
35  * Revision 1.4 2014/10/13 08:52:55 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.3 2014/10/06 15:13:07 j_novak
39  * Modified #include directives to use c++ syntax.
40  *
41  * Revision 1.2 2005/02/28 23:13:25 k_taniguchi
42  * Modification to include the case of the conformally flat background metric
43  *
44  * Revision 1.1 2004/11/30 20:49:13 k_taniguchi
45  * *** empty log message ***
46  *
47  *
48  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_excurv.C,v 1.4 2014/10/13 08:52:55 j_novak Exp $
49  *
50  */
51 
52 // C headers
53 #include <cmath>
54 
55 // Lorene headers
56 #include "et_bin_bhns_extr.h"
57 #include "etoile.h"
58 #include "coord.h"
59 #include "unites.h"
60 
61 namespace Lorene {
62 void Et_bin_bhns_extr::extrinsic_curv_extr(const double& mass,
63  const double& sepa) {
64 
65  using namespace Unites ;
66 
67  if (kerrschild) {
68 
77  // Components of shift_auto with respect to the Cartesian triad
78  // (d/dx, d/dy, d/dz) of the mapping :
79 
80  Tenseur shift_auto_local = shift_auto ;
81  shift_auto_local.change_triad( mp.get_bvect_cart() ) ;
82 
83  // Gradient (partial derivatives with respect to the Cartesian
84  // coordinates of the mapping)
85  // dn(i, j) = D_i N^j
86 
87  Tenseur dn = shift_auto_local.gradient() ;
88 
89  // Return to the absolute reference frame
91 
92  // Trace of D_i N^j = divergence of N^j :
93  Tenseur divn = contract(dn, 0, 1) ;
94 
95  // Computation of quantities coming from the companion (K-S BH)
96  // ------------------------------------------------------------
97 
98  const Coord& xx = mp.x ;
99  const Coord& yy = mp.y ;
100  const Coord& zz = mp.z ;
101 
102  Tenseur r_bh(mp) ;
103  r_bh.set_etat_qcq() ;
104  r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
105  r_bh.set_std_base() ;
106 
107  Tenseur xx_con(mp, 1, CON, ref_triad) ;
108  xx_con.set_etat_qcq() ;
109  xx_con.set(0) = xx + sepa ;
110  xx_con.set(1) = yy ;
111  xx_con.set(2) = zz ;
112  xx_con.set_std_base() ;
113 
114  Tenseur xsr_con(mp, 1, CON, ref_triad) ;
115  xsr_con = xx_con / r_bh ;
116  xsr_con.set_std_base() ;
117 
118  Tenseur msr(mp) ;
119  msr = ggrav * mass / r_bh ;
120  msr.set_std_base() ;
121 
122  Tenseur lapse_bh2(mp) ; // lapse_bh * lapse_bh
123  lapse_bh2 = 1. / (1.+2.*msr) ;
124  lapse_bh2.set_std_base() ;
125 
126  // Computation of some auxiliary functions
127  // ---------------------------------------
128 
129  shift_auto_local.change_triad( ref_triad ) ;
130 
131  Tenseur tmp1(mp, 2, CON, ref_triad) ;
132  Tenseur tmp2(mp, 2, CON, ref_triad) ;
133  Tenseur tmp3(mp, 2, CON, ref_triad) ;
134  tmp1.set_etat_qcq() ;
135  tmp2.set_etat_qcq() ;
136  tmp3.set_etat_qcq() ;
137 
138  for (int i=0; i<3; i++) {
139  for (int j=0; j<3; j++) {
140  tmp1.set(i, j) = -2.*lapse_bh2()%msr()%xsr_con(i)%xsr_con(j) ;
141 
142  tmp2.set(i, j) = -3.*lapse_bh2()%xsr_con(i)%xsr_con(j)
143  -4.*lapse_bh2()*msr()%xsr_con(i)%xsr_con(j) ;
144 
145  tmp3.set(i, j) = xsr_con(i)%shift_auto_local(j) ;
146  }
147  }
148 
149  tmp1.set_std_base() ;
150  tmp2.set_std_base() ;
151  tmp3.set_std_base() ;
152 
153  Tenseur tmp4(mp) ;
154  tmp4.set_etat_qcq() ;
155  tmp4.set() = 0 ;
156  tmp4.set_std_base() ;
157 
158  for (int i=0; i<3; i++)
159  tmp4.set() += xsr_con(i) % shift_auto_local(i) ;
160 
161  tmp4.set_std_base() ;
162 
163  // Computation of contraction
164  // --------------------------
165 
166  Tenseur tmp1dn = contract(tmp1, 1, dn, 0) ;
167 
168  // Computation of A^2 A^{ij}
169  // -------------------------
171 
172  for (int i=0; i<3; i++) {
173  for (int j=i; j<3; j++) {
174  tkij_auto.set(i, j) = dn(i, j) + dn(j, i)
175  + tmp1dn(i, j) + tmp1dn(j, i)
176  + 2.*lapse_bh2()%msr()/r_bh()%( tmp3(i, j) + tmp3(j, i)
177  + tmp4() % tmp2(i, j) )
178  - double(2)/double(3) * tmp1(i, j)
179  * (divn() - lapse_bh2() % msr() / r_bh() % tmp4()) ;
180  }
181  tkij_auto.set(i, i) -= double(2) /double(3)
182  * (divn() - lapse_bh2() % msr() / r_bh() % tmp4()) ;
183  }
184 
185  tkij_auto = - 0.5 * tkij_auto / nnn ;
186 
188 
189  // Computation of A^2 A_{ij} A^{ij}
190  // --------------------------------
191 
192  Tenseur xx_cov(mp, 1, COV, ref_triad) ;
193  xx_cov.set_etat_qcq() ;
194  xx_cov.set(0) = xx + sepa ;
195  xx_cov.set(1) = yy ;
196  xx_cov.set(2) = zz ;
197  xx_cov.set_std_base() ;
198 
199  Tenseur xsr_cov(mp, 1, COV, ref_triad) ;
200  xsr_cov = xx_cov / r_bh ;
201  xsr_cov.set_std_base() ;
202 
203  Tenseur tmp5(mp, 2, COV, ref_triad) ;
204  Tenseur tmp6(mp, 2, CON, ref_triad) ;
205  tmp5.set_etat_qcq() ;
206  tmp6.set_etat_qcq() ;
207 
208  for (int i=0; i<3; i++) {
209  for (int j=0; j<3; j++) {
210  tmp6.set(i, j) = 0 ;
211  }
212  }
213  tmp6.set_std_base() ;
214 
215  for (int i=0; i<3; i++) {
216  for (int j=0; j<3; j++) {
217  tmp5.set(i, j) = 2.*msr()%xsr_cov(i)%xsr_cov(j) ;
218  }
219  }
220 
221  for (int i=0; i<3; i++) {
222  for (int j=0; j<3; j++) {
223  for (int k=0; k<3; k++) {
224  tmp6.set(i, j) += tkij_auto(i,k) % tkij_auto(j,k) ;
225  }
226  }
227  }
228 
229  tmp5.set_std_base() ;
230  tmp6.set_std_base() ;
231 
232  Tenseur tmp7(mp) ;
233  tmp7.set_etat_qcq() ;
234  tmp7.set() = 0 ;
235  tmp7.set_std_base() ;
236 
237  for (int i=0; i<3; i++) {
238  for (int j=0; j<3; j++) {
239  tmp7.set() += tmp5(i,j) % tmp6(i,j) ;
240  }
241  }
242  tmp7.set_std_base() ;
243 
244  Tenseur tmp8(mp) ;
245  tmp8.set_etat_qcq() ;
246  tmp8.set() = 0 ;
247  tmp8.set_std_base() ;
248 
249  for (int i=0; i<3; i++) {
250  for (int j=0; j<3; j++) {
251  tmp8.set() += tmp5(i,j) % tkij_auto(i,j) ;
252  }
253  }
254  tmp8.set_std_base() ;
255 
257  akcar_auto.set() = 2.*tmp7() + tmp8() % tmp8() ;
259 
260  for (int i=0; i<3; i++) {
261  for (int j=0; j<3; j++) {
262  akcar_auto.set() += tkij_auto(i, j) % tkij_auto(i, j) ;
263  }
264  }
265 
267 
269 
270  }
271  else {
272 
282  // Components of shift_auto with respect to the Cartesian triad
283  // (d/dx, d/dy, d/dz) of the mapping :
284 
285  Tenseur shift_auto_local = shift_auto ;
286  shift_auto_local.change_triad( mp.get_bvect_cart() ) ;
287 
288  // Gradient (partial derivatives with respect to the Cartesian
289  // coordinates of the mapping)
290  // dn(i, j) = D_i N^j
291 
292  Tenseur dn = shift_auto_local.gradient() ;
293 
294  // Return to the absolute reference frame
295  dn.change_triad(ref_triad) ;
296 
297  // Trace of D_i N^j = divergence of N^j :
298  Tenseur divn = contract(dn, 0, 1) ;
299 
300  // Computation of A^2 A^{ij}
301  // -------------------------
303 
304  for (int i=0; i<3; i++) {
305  for (int j=i; j<3; j++) {
306  tkij_auto.set(i, j) = dn(i, j) + dn(j, i) ;
307  }
308  tkij_auto.set(i, i) -= double(2) /double(3) * divn() ;
309  }
310 
311  tkij_auto = - 0.5 * tkij_auto / nnn ;
312 
314 
315  // Computation of A^2 A_{ij} A^{ij}
316  // --------------------------------
317 
319  akcar_auto.set() = 0. ;
321 
322  for (int i=0; i<3; i++) {
323  for (int j=0; j<3; j++) {
324  akcar_auto.set() += tkij_auto(i, j) % tkij_auto(i, j) ;
325  }
326  }
327 
329 
331 
332  }
333 
334 }
335 }
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
bool kerrschild
Indicator of the background metric: true for the Kerr-Shild metric, false for the conformally flat on...
void extrinsic_curv_extr(const double &mass, const double &sepa)
Computes tkij_auto and akcar_auto from shift_auto , nnn and a_car .
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur 's...
Definition: etoile.h:828
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
Definition: etoile.h:889
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor generated by shift_auto .
Definition: etoile.h:925
Tenseur akcar_auto
Part of the scalar generated by shift_auto , i.e.
Definition: etoile.h:938
Tenseur nnn
Total lapse function.
Definition: etoile.h:509
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
Tenseur a_car
Total conformal factor .
Definition: etoile.h:515
Coord y
y coordinate centered on the grid
Definition: map.h:727
Coord z
z coordinate centered on the grid
Definition: map.h:728
Coord x
x coordinate centered on the grid
Definition: map.h:726
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition: map.h:791
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1542
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: tenseur.C:668
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.