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 
31 
32 /*
33  * $Id: et_bin_bhns_extr_excurv.C,v 1.5 2016/12/05 16:17:52 j_novak Exp $
34  * $Log: et_bin_bhns_extr_excurv.C,v $
35  * Revision 1.5 2016/12/05 16:17:52 j_novak
36  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37  *
38  * Revision 1.4 2014/10/13 08:52:55 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.3 2014/10/06 15:13:07 j_novak
42  * Modified #include directives to use c++ syntax.
43  *
44  * Revision 1.2 2005/02/28 23:13:25 k_taniguchi
45  * Modification to include the case of the conformally flat background metric
46  *
47  * Revision 1.1 2004/11/30 20:49:13 k_taniguchi
48  * *** empty log message ***
49  *
50  *
51  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_excurv.C,v 1.5 2016/12/05 16:17:52 j_novak Exp $
52  *
53  */
54 
55 // C headers
56 #include <cmath>
57 
58 // Lorene headers
59 #include "et_bin_bhns_extr.h"
60 #include "etoile.h"
61 #include "coord.h"
62 #include "unites.h"
63 
64 namespace Lorene {
65 void Et_bin_bhns_extr::extrinsic_curv_extr(const double& mass,
66  const double& sepa) {
67 
68  using namespace Unites ;
69 
70  if (kerrschild) {
71 
80  // Components of shift_auto with respect to the Cartesian triad
81  // (d/dx, d/dy, d/dz) of the mapping :
82 
83  Tenseur shift_auto_local = shift_auto ;
84  shift_auto_local.change_triad( mp.get_bvect_cart() ) ;
85 
86  // Gradient (partial derivatives with respect to the Cartesian
87  // coordinates of the mapping)
88  // dn(i, j) = D_i N^j
89 
90  Tenseur dn = shift_auto_local.gradient() ;
91 
92  // Return to the absolute reference frame
94 
95  // Trace of D_i N^j = divergence of N^j :
96  Tenseur divn = contract(dn, 0, 1) ;
97 
98  // Computation of quantities coming from the companion (K-S BH)
99  // ------------------------------------------------------------
100 
101  const Coord& xx = mp.x ;
102  const Coord& yy = mp.y ;
103  const Coord& zz = mp.z ;
104 
105  Tenseur r_bh(mp) ;
106  r_bh.set_etat_qcq() ;
107  r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
108  r_bh.set_std_base() ;
109 
110  Tenseur xx_con(mp, 1, CON, ref_triad) ;
111  xx_con.set_etat_qcq() ;
112  xx_con.set(0) = xx + sepa ;
113  xx_con.set(1) = yy ;
114  xx_con.set(2) = zz ;
115  xx_con.set_std_base() ;
116 
117  Tenseur xsr_con(mp, 1, CON, ref_triad) ;
118  xsr_con = xx_con / r_bh ;
119  xsr_con.set_std_base() ;
120 
121  Tenseur msr(mp) ;
122  msr = ggrav * mass / r_bh ;
123  msr.set_std_base() ;
124 
125  Tenseur lapse_bh2(mp) ; // lapse_bh * lapse_bh
126  lapse_bh2 = 1. / (1.+2.*msr) ;
127  lapse_bh2.set_std_base() ;
128 
129  // Computation of some auxiliary functions
130  // ---------------------------------------
131 
132  shift_auto_local.change_triad( ref_triad ) ;
133 
134  Tenseur tmp1(mp, 2, CON, ref_triad) ;
135  Tenseur tmp2(mp, 2, CON, ref_triad) ;
136  Tenseur tmp3(mp, 2, CON, ref_triad) ;
137  tmp1.set_etat_qcq() ;
138  tmp2.set_etat_qcq() ;
139  tmp3.set_etat_qcq() ;
140 
141  for (int i=0; i<3; i++) {
142  for (int j=0; j<3; j++) {
143  tmp1.set(i, j) = -2.*lapse_bh2()%msr()%xsr_con(i)%xsr_con(j) ;
144 
145  tmp2.set(i, j) = -3.*lapse_bh2()%xsr_con(i)%xsr_con(j)
146  -4.*lapse_bh2()*msr()%xsr_con(i)%xsr_con(j) ;
147 
148  tmp3.set(i, j) = xsr_con(i)%shift_auto_local(j) ;
149  }
150  }
151 
152  tmp1.set_std_base() ;
153  tmp2.set_std_base() ;
154  tmp3.set_std_base() ;
155 
156  Tenseur tmp4(mp) ;
157  tmp4.set_etat_qcq() ;
158  tmp4.set() = 0 ;
159  tmp4.set_std_base() ;
160 
161  for (int i=0; i<3; i++)
162  tmp4.set() += xsr_con(i) % shift_auto_local(i) ;
163 
164  tmp4.set_std_base() ;
165 
166  // Computation of contraction
167  // --------------------------
168 
169  Tenseur tmp1dn = contract(tmp1, 1, dn, 0) ;
170 
171  // Computation of A^2 A^{ij}
172  // -------------------------
174 
175  for (int i=0; i<3; i++) {
176  for (int j=i; j<3; j++) {
177  tkij_auto.set(i, j) = dn(i, j) + dn(j, i)
178  + tmp1dn(i, j) + tmp1dn(j, i)
179  + 2.*lapse_bh2()%msr()/r_bh()%( tmp3(i, j) + tmp3(j, i)
180  + tmp4() % tmp2(i, j) )
181  - double(2)/double(3) * tmp1(i, j)
182  * (divn() - lapse_bh2() % msr() / r_bh() % tmp4()) ;
183  }
184  tkij_auto.set(i, i) -= double(2) /double(3)
185  * (divn() - lapse_bh2() % msr() / r_bh() % tmp4()) ;
186  }
187 
188  tkij_auto = - 0.5 * tkij_auto / nnn ;
189 
191 
192  // Computation of A^2 A_{ij} A^{ij}
193  // --------------------------------
194 
195  Tenseur xx_cov(mp, 1, COV, ref_triad) ;
196  xx_cov.set_etat_qcq() ;
197  xx_cov.set(0) = xx + sepa ;
198  xx_cov.set(1) = yy ;
199  xx_cov.set(2) = zz ;
200  xx_cov.set_std_base() ;
201 
202  Tenseur xsr_cov(mp, 1, COV, ref_triad) ;
203  xsr_cov = xx_cov / r_bh ;
204  xsr_cov.set_std_base() ;
205 
206  Tenseur tmp5(mp, 2, COV, ref_triad) ;
207  Tenseur tmp6(mp, 2, CON, ref_triad) ;
208  tmp5.set_etat_qcq() ;
209  tmp6.set_etat_qcq() ;
210 
211  for (int i=0; i<3; i++) {
212  for (int j=0; j<3; j++) {
213  tmp6.set(i, j) = 0 ;
214  }
215  }
216  tmp6.set_std_base() ;
217 
218  for (int i=0; i<3; i++) {
219  for (int j=0; j<3; j++) {
220  tmp5.set(i, j) = 2.*msr()%xsr_cov(i)%xsr_cov(j) ;
221  }
222  }
223 
224  for (int i=0; i<3; i++) {
225  for (int j=0; j<3; j++) {
226  for (int k=0; k<3; k++) {
227  tmp6.set(i, j) += tkij_auto(i,k) % tkij_auto(j,k) ;
228  }
229  }
230  }
231 
232  tmp5.set_std_base() ;
233  tmp6.set_std_base() ;
234 
235  Tenseur tmp7(mp) ;
236  tmp7.set_etat_qcq() ;
237  tmp7.set() = 0 ;
238  tmp7.set_std_base() ;
239 
240  for (int i=0; i<3; i++) {
241  for (int j=0; j<3; j++) {
242  tmp7.set() += tmp5(i,j) % tmp6(i,j) ;
243  }
244  }
245  tmp7.set_std_base() ;
246 
247  Tenseur tmp8(mp) ;
248  tmp8.set_etat_qcq() ;
249  tmp8.set() = 0 ;
250  tmp8.set_std_base() ;
251 
252  for (int i=0; i<3; i++) {
253  for (int j=0; j<3; j++) {
254  tmp8.set() += tmp5(i,j) % tkij_auto(i,j) ;
255  }
256  }
257  tmp8.set_std_base() ;
258 
260  akcar_auto.set() = 2.*tmp7() + tmp8() % tmp8() ;
262 
263  for (int i=0; i<3; i++) {
264  for (int j=0; j<3; j++) {
265  akcar_auto.set() += tkij_auto(i, j) % tkij_auto(i, j) ;
266  }
267  }
268 
270 
272 
273  }
274  else {
275 
285  // Components of shift_auto with respect to the Cartesian triad
286  // (d/dx, d/dy, d/dz) of the mapping :
287 
288  Tenseur shift_auto_local = shift_auto ;
289  shift_auto_local.change_triad( mp.get_bvect_cart() ) ;
290 
291  // Gradient (partial derivatives with respect to the Cartesian
292  // coordinates of the mapping)
293  // dn(i, j) = D_i N^j
294 
295  Tenseur dn = shift_auto_local.gradient() ;
296 
297  // Return to the absolute reference frame
298  dn.change_triad(ref_triad) ;
299 
300  // Trace of D_i N^j = divergence of N^j :
301  Tenseur divn = contract(dn, 0, 1) ;
302 
303  // Computation of A^2 A^{ij}
304  // -------------------------
306 
307  for (int i=0; i<3; i++) {
308  for (int j=i; j<3; j++) {
309  tkij_auto.set(i, j) = dn(i, j) + dn(j, i) ;
310  }
311  tkij_auto.set(i, i) -= double(2) /double(3) * divn() ;
312  }
313 
314  tkij_auto = - 0.5 * tkij_auto / nnn ;
315 
317 
318  // Computation of A^2 A_{ij} A^{ij}
319  // --------------------------------
320 
322  akcar_auto.set() = 0. ;
324 
325  for (int i=0; i<3; i++) {
326  for (int j=0; j<3; j++) {
327  akcar_auto.set() += tkij_auto(i, j) % tkij_auto(i, j) ;
328  }
329  }
330 
332 
334 
335  }
336 
337 }
338 }
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur &#39;s...
Definition: etoile.h:831
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1186
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
Tenseur nnn
Total lapse function.
Definition: etoile.h:512
bool kerrschild
Indicator of the background metric: true for the Kerr-Shild metric, false for the conformally flat on...
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
Definition: etoile.h:892
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: tenseur.C:684
Map & mp
Mapping associated with the star.
Definition: etoile.h:432
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor generated by shift_auto .
Definition: etoile.h:928
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
void extrinsic_curv_extr(const double &mass, const double &sepa)
Computes tkij_auto and akcar_auto from shift_auto , nnn and a_car .
Tenseur a_car
Total conformal factor .
Definition: etoile.h:518
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
Coord y
y coordinate centered on the grid
Definition: map.h:745
Coord x
x coordinate centered on the grid
Definition: map.h:744
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:652
Tenseur akcar_auto
Part of the scalar generated by shift_auto , i.e.
Definition: etoile.h:941
Coord z
z coordinate centered on the grid
Definition: map.h:746
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1558