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 char 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 
61 namespace 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 ;
74  rr.std_spectral_base() ;
75  Scalar st(mp) ;
76  st = mp.sint ;
77  st.std_spectral_base() ;
78  Scalar ct(mp) ;
79  ct = mp.cost ;
80  ct.std_spectral_base() ;
81  Scalar sp(mp) ;
82  sp = mp.sinp ;
83  sp.std_spectral_base() ;
84  Scalar cp(mp) ;
85  cp = mp.cosp ;
86  cp.std_spectral_base() ;
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 ;
93  ll.std_spectral_base() ;
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 ;
178  taij.annule_domain(0) ;
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 
419  taij.annule_domain(0) ;
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 }
Lorene::Scalar::raccord
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Definition: scalar_raccord.C:62
Lorene::Metric_flat::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:153
Lorene::Metric_flat::cov
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric_flat.C:134
Lorene::Black_hole::extr_curv_bh
void extr_curv_bh()
Computes taij , taij_quad from shift , lapse , confo .
Definition: blackhole_extr_curv.C:62
Lorene::Scalar::inc_dzpuis
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...
Definition: scalar_r_manip.C:470
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
Lorene::Map::cost
Coord cost
Definition: map.h:722
Lorene::Black_hole::shift_rs
Vector shift_rs
Part of the shift vector from the numerical computation.
Definition: blackhole.h:112
Lorene::Black_hole::flat
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the black hol...
Definition: blackhole.h:143
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Black_hole::taij_rs
Sym_tensor taij_rs
Part of the extrinsic curvature tensor.
Definition: blackhole.h:130
Lorene::Black_hole::lapconf
Scalar lapconf
A function (lapse function * conformal factor) lapconf generated by the black hole.
Definition: blackhole.h:97
Lorene::Black_hole::mass_bh
double mass_bh
Gravitational mass of BH.
Definition: blackhole.h:88
Lorene::Map::sint
Coord sint
Definition: map.h:721
Lorene::Map::r
Coord r
r coordinate centered on the grid
Definition: map.h:718
Lorene::Black_hole::mp
Map & mp
Mapping associated with the black hole.
Definition: blackhole.h:80
Lorene::Black_hole::taij
Sym_tensor taij
Trace of the extrinsic curvature.
Definition: blackhole.h:127
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
Lorene::Tensor::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:481
Unites
Standard units of space, time and mass.
Lorene::Map::sinp
Coord sinp
Definition: map.h:723
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
Lorene::Tensor::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: tensor.C:926
Lorene::Black_hole::taij_quad
Scalar taij_quad
Part of the scalar generated by .
Definition: blackhole.h:135
Lorene::Map::get_bvect_cart
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
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Lorene::Vector
Tensor field of valence 1.
Definition: vector.h:188
Lorene::Tensor::inc_dzpuis
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
Lorene::Scalar::std_spectral_base
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
Lorene::Black_hole::confo
Scalar confo
Conformal factor generated by the black hole.
Definition: blackhole.h:118
Lorene::Black_hole::taij_quad_rs
Scalar taij_quad_rs
Part of the scalar.
Definition: blackhole.h:138
Lorene::Map::cosp
Coord cosp
Definition: map.h:724
Lorene::Black_hole::shift
Vector shift
Shift vector generated by the black hole.
Definition: blackhole.h:109
Lorene::Black_hole::kerrschild
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition: blackhole.h:85
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:316