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
30char 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
61namespace Lorene {
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
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
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
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
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.