LORENE
blackhole_extr_curv.C
1/*
2 * Method of class Black_hole to compute the extrinsic curvature tensor
3 *
4 * (see file blackhole.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2005-2007 Keisuke Taniguchi
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char blackhole_extr_curv_C[] = "$Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_extr_curv.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $" ;
29
30/*
31 * $Id: blackhole_extr_curv.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $
32 * $Log: blackhole_extr_curv.C,v $
33 * Revision 1.4 2014/10/13 08:52:46 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.3 2014/10/06 15:13:02 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.2 2008/05/15 19:27:14 k_taniguchi
40 * Change of some parameters.
41 *
42 * Revision 1.1 2007/06/22 01:19:32 k_taniguchi
43 * *** empty log message ***
44 *
45 *
46 * $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_extr_curv.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $
47 *
48 */
49
50// C++ headers
51//#include <>
52
53// C headers
54#include <cmath>
55
56// Lorene headers
57#include "blackhole.h"
58#include "utilitaires.h"
59#include "unites.h"
60
61namespace Lorene {
63
64 // Fundamental constants and units
65 // -------------------------------
66 using namespace Unites ;
67
68 if (kerrschild) {
69
70 double mass = ggrav * mass_bh ;
71
72 Scalar rr(mp) ;
73 rr = mp.r ;
75 Scalar st(mp) ;
76 st = mp.sint ;
78 Scalar ct(mp) ;
79 ct = mp.cost ;
81 Scalar sp(mp) ;
82 sp = mp.sinp ;
84 Scalar cp(mp) ;
85 cp = mp.cosp ;
87
88 Vector ll(mp, CON, mp.get_bvect_cart()) ;
89 ll.set_etat_qcq() ;
90 ll.set(1) = st * cp ;
91 ll.set(2) = st * sp ;
92 ll.set(3) = ct ;
94
95
96 // Computation of \tilde{A}^{ij}
97 // -----------------------------
98
99 Scalar divshift(mp) ;
100 divshift = shift_rs(1).dsdx() + shift_rs(2).dsdy()
101 + shift_rs(3).dsdz() ;
102 divshift.std_spectral_base() ;
103
104 Sym_tensor flat_taij(mp, CON, mp.get_bvect_cart()) ;
105 flat_taij.set_etat_qcq() ;
106
107 for (int i=1; i<=3; i++) {
108 for (int j=1; j<=3; j++) {
109 flat_taij.set(i,j) = shift_rs(j).deriv(i)
110 + shift_rs(i).deriv(j)
111 - 2. * divshift * flat.con()(i,j) / 3. ;
112 }
113 }
114
115 flat_taij.std_spectral_base() ;
116
117 Scalar lapse_bh(mp) ;
118 lapse_bh = 1. / sqrt(1. + 2. * mass / rr) ;
119 lapse_bh.std_spectral_base() ;
120 lapse_bh.annule_domain(0) ;
121 lapse_bh.raccord(1) ;
122
123 Sym_tensor curv_taij(mp, CON, mp.get_bvect_cart()) ;
124 curv_taij.set_etat_qcq() ;
125
126 for (int i=1; i<=3; i++) {
127 for (int j=1; j<=3; j++) {
128 curv_taij.set(i,j) = -2. * lapse_bh * lapse_bh * mass
129 * (ll(i)*(shift_rs(j).dsdr()) + ll(j)*(shift_rs(i).dsdr())
130 - 2. * ll(i) * ll(j) * divshift / 3.) / rr ;
131 }
132 }
133
134 curv_taij.std_spectral_base() ;
135
136 Sym_tensor resi_taij(mp, CON, mp.get_bvect_cart()) ;
137 resi_taij.set_etat_qcq() ;
138
139 for (int i=1; i<=3; i++) {
140 for (int j=1; j<=3; j++) {
141 resi_taij.set(i,j) = 2. * lapse_bh * lapse_bh * mass
142 * ( ll(i) * shift_rs(j) + ll(j) * shift_rs(i)
143 + ( flat.con()(i,j)
144 - lapse_bh*lapse_bh*(9.+14.*mass/rr)*ll(i)*ll(j) )
145 * ( ll(1)*shift_rs(1)+ll(2)*shift_rs(2)
146 +ll(3)*shift_rs(3) )/ 3. )
147 / rr / rr ;
148 }
149 }
150
151 resi_taij.std_spectral_base() ;
152 resi_taij.inc_dzpuis(2) ;
153
154 taij_rs = 0.5 * pow(confo, 7.)
155 * (flat_taij + curv_taij + resi_taij) / lapconf ;
156
159
160 Sym_tensor taij_bh(mp, CON, mp.get_bvect_cart()) ;
161 taij_bh.set_etat_qcq() ;
162
163 for (int i=1; i<=3; i++) {
164 for (int j=1; j<=3; j++) {
165 taij_bh.set(i,j) = 2.*pow(lapse_bh,6.)*mass*(2.+3.*mass/rr)
166 *( (1.+2.*mass/rr) * flat.con()(i,j)
167 - (3.+2.*mass/rr) * ll(i) * ll(j) )
168 *pow(confo, 7.)/lapconf/3./rr/rr ;
169 }
170 }
171
172 taij_bh.std_spectral_base() ;
173 taij_bh.inc_dzpuis(2) ;
174 taij_bh.annule_domain(0) ;
175
176 taij = taij_rs + taij_bh ;
179
180 /*
181 Sym_tensor taij_ks_con(mp, CON, mp.get_bvect_cart()) ;
182 taij_ks_con.set_etat_qcq() ;
183
184 for (int i=1; i<=3; i++) {
185 for (int j=1; j<=3; j++) {
186 taij_ks_con.set(i,j) = 2.*pow(lap_bh2,2.5)*mass
187 * (2.+3.*mass/rr)
188 * ( (1.+2.*mass/rr)*flat.con()(i,j)
189 - (3.+2.*mass/rr)*ll(i)*ll(j) ) / 3. / rr / rr ;
190 }
191 }
192 taij_ks_con.std_spectral_base() ;
193 taij_ks_con.annule_domain(0) ;
194 taij_ks_con.inc_dzpuis(2) ;
195
196 cout << "taij(1,1) - taij_ks_con(1,1) : " << endl ;
197 cout << taij(1,1) - taij_ks_con(1,1) << endl ;
198 arrete() ;
199
200 cout << "taij_ks_con(1,1) : " << endl ;
201 cout << taij_ks_con(1,1) << endl ;
202 arrete() ;
203
204 cout << "taij(1,1) : " << endl ;
205 cout << taij(1,1) << endl ;
206 arrete() ;
207 */
208
209 // Computation of \tilde{A}^{ij} \tilde{A}_{ij}
210 // --------------------------------------------
211
212 Sym_tensor flat_dshift(mp, COV, mp.get_bvect_cart()) ;
213 flat_dshift.set_etat_qcq() ;
214
215 for (int i=1; i<=3; i++) {
216 for (int j=1; j<=3; j++) {
217 flat_dshift.set(i,j) = flat.cov()(j,1)*(shift_rs(1).deriv(i))
218 + flat.cov()(j,2)*(shift_rs(2).deriv(i))
219 + flat.cov()(j,3)*(shift_rs(3).deriv(i))
220 + flat.cov()(i,1)*(shift_rs(1).deriv(j))
221 + flat.cov()(i,2)*(shift_rs(2).deriv(j))
222 + flat.cov()(i,3)*(shift_rs(3).deriv(j))
223 - 2. * divshift * flat.cov()(i,j) / 3. ;
224 }
225 }
226
227 flat_dshift.std_spectral_base() ;
228
229 Sym_tensor curv_dshift(mp, COV, mp.get_bvect_cart()) ;
230 curv_dshift.set_etat_qcq() ;
231
232 for (int i=1; i<=3; i++) {
233 for (int j=1; j<=3; j++) {
234 curv_dshift.set(i,j) = 2. * mass
235 *( ll(j) *( ll(1)*(shift_rs(1).deriv(i))
236 + ll(2)*(shift_rs(2).deriv(i))
237 + ll(3)*(shift_rs(3).deriv(i)) )
238 + ll(i) *( ll(1)*(shift_rs(1).deriv(j))
239 + ll(2)*(shift_rs(2).deriv(j))
240 + ll(3)*(shift_rs(3).deriv(j)) )
241 - 2. * divshift * ll(i) * ll(j) / 3. ) / rr ;
242 }
243 }
244
245 curv_dshift.std_spectral_base() ;
246
247 Sym_tensor tmp1(mp, COV, mp.get_bvect_cart()) ;
248 tmp1.set_etat_qcq() ;
249
250 for (int i=1; i<=3; i++) {
251 for (int j=1; j<=3; j++) {
252 tmp1.set(i,j) = 2. * mass
253 *(ll(j)*( (flat.cov()(i,1)+2.*mass*ll(i)*ll(1)/rr)
254 *shift_rs(1)
255 + (flat.cov()(i,2)+2.*mass*ll(i)*ll(2)/rr)
256 *shift_rs(2)
257 + (flat.cov()(i,3)+2.*mass*ll(i)*ll(3)/rr)
258 *shift_rs(3)
259 )
260 + ll(i)*( (flat.cov()(j,1)+2.*mass*ll(j)*ll(1)/rr)
261 *shift_rs(1)
262 + (flat.cov()(j,2)+2.*mass*ll(j)*ll(2)/rr)
263 *shift_rs(2)
264 + (flat.cov()(j,3)+2.*mass*ll(j)*ll(3)/rr)
265 *shift_rs(3) )
266 ) / rr / rr ;
267 }
268 }
269 tmp1.std_spectral_base() ;
270 tmp1.inc_dzpuis(2) ;
271
272 Sym_tensor tmp2(mp, COV, mp.get_bvect_cart()) ;
273 tmp2.set_etat_qcq() ;
274
275 for (int i=1; i<=3; i++) {
276 for (int j=1; j<=3; j++) {
277 tmp2.set(i,j) = 2. * mass * lapse_bh * lapse_bh
278 * (ll(1)*shift_rs(1)+ll(2)*shift_rs(2)+ll(3)*shift_rs(3))
279 * (flat.cov()(i,j)
280 - (9.+28.*mass/rr+24.*mass*mass/rr/rr)*ll(i)*ll(j))
281 / 3. / rr / rr ;
282 }
283 }
284 tmp2.std_spectral_base() ;
285 tmp2.inc_dzpuis(2) ;
286
287 Sym_tensor taij_rs_down(mp, COV, mp.get_bvect_cart()) ;
288 taij_rs_down.set_etat_qcq() ;
289
290 taij_rs_down = 0.5 * pow(confo, 7.)
291 * (flat_dshift + curv_dshift + tmp1 + tmp2) / lapconf ;
292
293 taij_rs_down.std_spectral_base() ;
294 taij_rs_down.annule_domain(0) ;
295
296 Sym_tensor taij_bh_down(mp, COV, mp.get_bvect_cart()) ;
297 taij_bh_down.set_etat_qcq() ;
298
299 for (int i=1; i<=3; i++) {
300 for (int j=1; j<=3; j++) {
301 taij_bh_down.set(i,j) = 2.*pow(lapse_bh,4.)*mass*(2.+3.*mass/rr)
302 *pow(confo,7.)*(flat.cov()(i,j)-(3.+4.*mass/rr)*ll(i)*ll(j))
303 /lapconf/3./rr/rr ;
304 }
305 }
306
307 taij_bh_down.std_spectral_base() ;
308 taij_bh_down.inc_dzpuis(2) ;
309 taij_bh_down.annule_domain(0) ;
310
311 /*
312 Sym_tensor taij_ks_cov(mp, COV, mp.get_bvect_cart()) ;
313 taij_ks_cov.set_etat_qcq() ;
314
315 for (int i=1; i<=3; i++) {
316 for (int j=1; j<=3; j++) {
317 taij_ks_cov.set(i,j) = 2.*pow(lap_bh2,1.5)*mass * (2.+3.*mass/rr)
318 * ( flat.cov()(i,j)
319 - (3.+4.*mass/rr)*ll(i)*ll(j) ) / 3. / rr / rr ;
320 }
321 }
322 taij_ks_cov.std_spectral_base() ;
323 taij_ks_cov.annule_domain(0) ;
324 taij_ks_cov.inc_dzpuis(2) ;
325
326 cout << "taij_down(1,1) - taij_ks_cov(1,1) : " << endl ;
327 cout << taij_down(1,1) - taij_ks_cov(1,1) << endl ;
328 arrete() ;
329
330 cout << "taij_ks_cov(1,1) : " << endl ;
331 cout << taij_ks_cov(1,1) << endl ;
332 arrete() ;
333
334 cout << "taij_down(1,1) : " << endl ;
335 cout << taij_down(1,1) << endl ;
336 arrete() ;
337 */
338
339 Scalar taij_quad_rsrs(mp) ;
340 taij_quad_rsrs = 0. ;
341
342 for (int i=1; i<=3; i++) {
343 for (int j=1; j<=3; j++) {
344 taij_quad_rsrs += taij_rs_down(i,j) * taij_rs(i,j) ;
345 }
346 }
347 taij_quad_rsrs.std_spectral_base() ;
348
349 Scalar taij_quad_rsbh1(mp) ;
350 taij_quad_rsbh1 = 0. ;
351
352 for (int i=1; i<=3; i++) {
353 for (int j=1; j<=3; j++) {
354 taij_quad_rsbh1 += taij_rs_down(i,j) * taij_bh(i,j) ;
355 }
356 }
357 taij_quad_rsbh1.std_spectral_base() ;
358
359 Scalar taij_quad_rsbh2(mp) ;
360 taij_quad_rsbh2 = 0. ;
361
362 for (int i=1; i<=3; i++) {
363 for (int j=1; j<=3; j++) {
364 taij_quad_rsbh2 += taij_bh_down(i,j) * taij_rs(i,j) ;
365 }
366 }
367 taij_quad_rsbh2.std_spectral_base() ;
368
369 taij_quad_rs = taij_quad_rsrs + taij_quad_rsbh1 + taij_quad_rsbh2 ;
371
372 Scalar taij_quad_bh(mp) ;
373 taij_quad_bh = 8.*pow(lapse_bh,10.)*mass*mass*(2.+3.*mass/rr)
374 *(2.+3.*mass/rr)*pow(confo,12.)/3./pow(rr,4.)/lapconf/lapconf ;
375 taij_quad_bh.std_spectral_base() ;
376 taij_quad_bh.inc_dzpuis(4) ;
377
378 taij_quad = taij_quad_rs + taij_quad_bh ;
379
381
382 /*
383 Scalar taij_quad_ks(mp) ;
384 taij_quad_ks = 8. * pow(lap_bh2,3.) * mass * mass * (2.+3.*mass/rr)
385 * (2.+3.*mass/rr) / 3. / pow(rr, 4.) ;
386 taij_quad_ks.std_spectral_base() ;
387 taij_quad_ks.annule_domain(0) ;
388 taij_quad_ks.inc_dzpuis(4) ;
389
390 cout << "taij_quad - taij_quad_ks : " << endl ;
391 cout << taij_quad - taij_quad_ks << endl ;
392 arrete() ;
393 */
394 }
395 else { // Isotropic coordinates with the maximal slicing
396
397 // Computation of \tilde{A}^{ij}
398 // -----------------------------
399
400 Scalar divshift(mp) ;
401 divshift = shift(1).dsdx() + shift(2).dsdy() + shift(3).dsdz() ;
402 divshift.std_spectral_base() ;
403
404 Sym_tensor flat_taij(mp, CON, mp.get_bvect_cart()) ;
405 flat_taij.set_etat_qcq() ;
406
407 for (int i=1; i<=3; i++) {
408 for (int j=1; j<=3; j++) {
409 flat_taij.set(i,j) = shift(j).deriv(i) + shift(i).deriv(j)
410 - 2. * divshift * flat.con()(i,j) / 3. ;
411 }
412 }
413
414 flat_taij.std_spectral_base() ;
415
416 taij = 0.5 * pow(confo, 7.) * flat_taij / lapconf ;
417
420
421
422 // Computation of \tilde{A}^{ij} \tilde{A}_{ij}
423 // --------------------------------------------
424
425 Sym_tensor flat_dshift(mp, COV, mp.get_bvect_cart()) ;
426 flat_dshift.set_etat_qcq() ;
427
428 for (int i=1; i<=3; i++) {
429 for (int j=1; j<=3; j++) {
430 flat_dshift.set(i,j) = flat.cov()(j,1)*(shift(1).deriv(i))
431 + flat.cov()(j,2)*(shift(2).deriv(i))
432 + flat.cov()(j,3)*(shift(3).deriv(i))
433 + flat.cov()(i,1)*(shift(1).deriv(j))
434 + flat.cov()(i,2)*(shift(2).deriv(j))
435 + flat.cov()(i,3)*(shift(3).deriv(j))
436 - 2. * divshift * flat.cov()(i,j) / 3. ;
437 }
438 }
439
440 Sym_tensor taij_down(mp, COV, mp.get_bvect_cart()) ;
441 taij_down.set_etat_qcq() ;
442
443 for (int i=1; i<=3; i++) {
444 for (int j=1; j<=3; j++) {
445 taij_down.set(i,j) = 0.5 * pow(confo, 7.) * flat_dshift(i,j)
446 / lapconf ;
447 }
448 }
449
450 taij_down.std_spectral_base() ;
451 taij_down.annule_domain(0) ;
452
453 taij_quad = 0. ;
454
455 for (int i=1; i<=3; i++) {
456 for (int j=1; j<=3; j++) {
457 taij_quad += taij_down(i,j) * taij(i,j) ;
458 }
459 }
461
462 }
463
464}
465}
Vector shift_rs
Part of the shift vector from the numerical computation.
Definition blackhole.h:112
Scalar taij_quad
Part of the scalar generated by .
Definition blackhole.h:135
Scalar taij_quad_rs
Part of the scalar.
Definition blackhole.h:138
Sym_tensor taij
Trace of the extrinsic curvature.
Definition blackhole.h:127
Sym_tensor taij_rs
Part of the extrinsic curvature tensor.
Definition blackhole.h:130
Map & mp
Mapping associated with the black hole.
Definition blackhole.h:80
Scalar lapconf
A function (lapse function * conformal factor) lapconf generated by the black hole.
Definition blackhole.h:97
Vector shift
Shift vector generated by the black hole.
Definition blackhole.h:109
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the black hol...
Definition blackhole.h:143
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition blackhole.h:85
Scalar confo
Conformal factor generated by the black hole.
Definition blackhole.h:118
double mass_bh
Gravitational mass of BH.
Definition blackhole.h:88
void extr_curv_bh()
Computes taij , taij_quad from shift , lapse , confo .
Coord cosp
Definition map.h:724
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 sint
Definition map.h:721
Coord r
r coordinate centered on the grid
Definition map.h:718
Coord sinp
Definition map.h:723
Coord cost
Definition map.h:722
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Tensor field of valence 1.
Definition vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:816
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:481
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition tensor.C:926
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.