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