LORENE
blackhole_killing.C
1 /*
2  * Methods of class Black_hole to compute Killing vectors
3  *
4  * (see file blackhole.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 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_killing.C,v 1.5 2016/12/05 16:17:48 j_novak Exp $
32  * $Log: blackhole_killing.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/07/02 20:45:07 k_taniguchi
43  * A bug removed.
44  *
45  * Revision 1.1 2008/05/15 19:33:12 k_taniguchi
46  * *** empty log message ***
47  *
48  *
49  * $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_killing.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 "unites.h"
62 #include "utilitaires.h"
63 
64  //---------------------------------------------//
65  // Killing vectors on the AH //
66  //---------------------------------------------//
67 
68 namespace Lorene {
69 Vector Black_hole::killing_vect_bh(const Tbl& xi_i, const double& phi_i,
70  const double& theta_i, const int& nrk_phi,
71  const int& nrk_theta) const {
72 
73  using namespace Unites ;
74 
75  assert(xi_i.get_ndim() == 1) ;
76  assert(xi_i.get_dim(0) == 3) ;
77 
78  const Mg3d* mg = mp.get_mg() ;
79  int nr = mg->get_nr(1) ;
80  int nt = mg->get_nt(1) ;
81  int np = mg->get_np(1) ;
82 
83  // Vector which is returned to the main code
84  // Spherical basis, covariant
85  Vector killing(mp, COV, mp.get_bvect_spher()) ;
86 
87  if (kerrschild) {
88 
89  cout << "Not yet prepared!!!" << endl ;
90  abort() ;
91 
92  }
93  else { // Isotropic coordinates
94 
95  // Solution of the Killing vector on the equator
96  // ---------------------------------------------
97 
98  double dp = 2. * M_PI / double(np) ; // Angular step
99 
100  // Killing vector on the equator
101  // np+1 is for the check of xi(phi=0)= xi(phi=2pi)
102  Tbl xi_t(np+1) ;
103  xi_t.set_etat_qcq() ;
104  Tbl xi_p(np+1) ;
105  xi_p.set_etat_qcq() ;
106  Tbl xi_l(np+1) ;
107  xi_l.set_etat_qcq() ;
108 
109  xi_t.set(0) = xi_i(0) ;
110  xi_p.set(0) = xi_i(1) ;
111  xi_l.set(0) = xi_i(2) ;
112 
113  Tbl xi(3) ;
114  xi.set_etat_qcq() ;
115 
116  Tbl xi_ini(3) ;
117  xi_ini.set_etat_qcq() ;
118  xi_ini.set(0) = xi_i(0) ;
119  xi_ini.set(1) = xi_i(1) ;
120  xi_ini.set(2) = xi_i(2) ;
121 
122  double pp_0 = phi_i ; // azimuthal angle phi
123 
124  for (int i=1; i<np+1; i++) {
125 
126  xi = runge_kutta_phi_bh(xi_ini, pp_0, nrk_phi) ;
127 
128  xi_t.set(i) = xi(0) ;
129  xi_p.set(i) = xi(1) ;
130  xi_l.set(i) = xi(2) ;
131 
132  // New data for the next step
133  // -------------------------
134  pp_0 += dp ; // New angle
135  xi_ini = xi ;
136 
137  }
138 
139  /*
140  for (int i=0; i<np+1; i++) {
141 
142  cout << "xi_t xi_p xi_l" << endl ;
143  cout << xi_t(i) << " " << xi_p(i) << " " << xi_l(i) << endl ;
144 
145  }
146  arrete() ;
147  */
148 
149  // Normalization of the Killing vector
150  // -----------------------------------
151 
152  // Initialization of the Killing vector to the phi direction
153  Scalar xi_phi(mp) ;
154  xi_phi = 0.5 ; // If we use "1." for the initialization value,
155  // the state flag becomes "ETATUN" which does not
156  // work for "set_grid_point".
157 
158  for (int k=0; k<np; k++) {
159  xi_phi.set_grid_point(0, k, nt-1, nr-1) = xi_p(k) ;
160  xi_phi.set_grid_point(1, k, nt-1, 0) = xi_p(k) ;
161  }
162  xi_phi.std_spectral_base() ;
163  /*
164  for (int l=0; l<2; l++) {
165  for (int k=0; k<np; k++) {
166  for (int j=0; j<nt; j++) {
167  for (int i=0; i<nr; i++) {
168  cout << "(l,k,j,i)=" << l << "," << k << "," << j
169  << "," << i << ": "
170  << xi_phi.val_grid_point(l,k,j,i) << endl ;
171  }
172  arrete() ;
173  }
174  arrete() ;
175  }
176  arrete() ;
177  }
178  */
179 
180  // Normalization of the Killing vector
181  Scalar rr(mp) ;
182  rr = mp.r ;
183  rr.std_spectral_base() ;
184 
185  Scalar st(mp) ;
186  st = mp.sint ;
187  st.std_spectral_base() ;
188 
189  Scalar source_phi(mp) ;
190  source_phi = pow(confo, 2.) * rr * st / xi_phi ;
191  source_phi.std_spectral_base() ;
192 
193  double rah = rad_ah() ;
194 
195  int nn = 1000 ;
196  double hh = 2. * M_PI / double(nn) ;
197  double integ = 0. ;
198 
199  int mm ;
200  double t1, t2, t3, t4, t5 ;
201 
202  // Boole's Rule (Newton-Cotes Integral) for integration
203  // ----------------------------------------------------
204 
205  assert(nn%4 == 0) ;
206  mm = nn/4 ;
207 
208  for (int i=0; i<mm; i++) {
209 
210  t1 = hh * double(4*i) ;
211  t2 = hh * double(4*i+1) ;
212  t3 = hh * double(4*i+2) ;
213  t4 = hh * double(4*i+3) ;
214  t5 = hh * double(4*i+4) ;
215 
216  integ += (hh/45.) * (14.*source_phi.val_point(rah,M_PI/2.,t1)
217  + 64.*source_phi.val_point(rah,M_PI/2.,t2)
218  + 24.*source_phi.val_point(rah,M_PI/2.,t3)
219  + 64.*source_phi.val_point(rah,M_PI/2.,t4)
220  + 14.*source_phi.val_point(rah,M_PI/2.,t5)
221  ) ;
222 
223  }
224 
225  cout << "Black_hole:: t_f = " << integ << endl ;
226  double ratio = 0.5 * integ / M_PI ;
227 
228  cout << "Black_hole:: t_f / 2M_PI = " << ratio << endl ;
229 
230  for (int k=0; k<np; k++) {
231  xi_p.set(k) = xi_phi.val_grid_point(1, k, nt-1, 0) * ratio ;
232  }
233 
234  /*
235  for (int k=0; k<np; k++) {
236  cout << "Normalized xi_p" << "(" << k << ") :" << xi_p(k) << endl ;
237  }
238  */
239 
240  // Solution of the Killing vector to the pole angle
241  // ------------------------------------------------
242 
243  double dt = 0.5 * M_PI / double(nt-1) ; // Angular step
244 
245  // Killing vector to the polar angle
246  Tbl xi_th(nt, np) ;
247  xi_th.set_etat_qcq() ;
248  Tbl xi_ph(nt, np) ;
249  xi_ph.set_etat_qcq() ;
250  Tbl xi_ll(nt, np) ;
251  xi_ll.set_etat_qcq() ;
252 
253  // Values on the equator
254  for (int i=0; i<np; i++) {
255 
256  xi_th.set(nt-1, i) = xi_t(i) ;
257  xi_ph.set(nt-1, i) = xi_p(i) ;
258  xi_ll.set(nt-1, i) = xi_l(i) ;
259 
260  }
261 
262  for (int i=0; i<np; i++) {
263 
264  // Value at theta=pi/2, phi=phi_i
265  xi_ini.set(0) = xi_t(i) ;
266  xi_ini.set(1) = xi_p(i) ;
267  xi_ini.set(2) = xi_l(i) ;
268 
269  double pp = double(i) * dp ;
270  double tt_0 = theta_i ; // polar angle theta
271 
272  for (int j=1; j<nt; j++) {
273 
274  xi = runge_kutta_theta_bh(xi_ini, tt_0, pp, nrk_theta) ;
275 
276  xi_th.set(nt-1-j, i) = xi(0) ;
277  xi_ph.set(nt-1-j, i) = xi(1) ;
278  xi_ll.set(nt-1-j, i) = xi(2) ;
279 
280  // New data for the nxt step
281  // -------------------------
282  tt_0 -= dt ; // New angle
283  xi_ini = xi ;
284 
285  } // End of the loop to the polar direction
286 
287  } // End of the loop to the azimuhtal direction
288 
289 
290  // Construction of the Killing vector
291  // ----------------------------------
292 
293  // Definition of the vector is at the top of this code
294  killing.set_etat_qcq() ;
295  killing.set(1) = 0. ; // r component
296  killing.set(2) = 0.5 ; // initialization of the theta component
297  killing.set(3) = 0.5 ; // initialization of the phi component
298 
299  for (int l=0; l<2; l++) {
300  for (int i=0; i<nr; i++) {
301  for (int j=0; j<nt; j++) {
302  for (int k=0; k<np; k++) {
303  (killing.set(2)).set_grid_point(l, k, j, i) = xi_th(j, k) ;
304  (killing.set(3)).set_grid_point(l, k, j, i) = xi_ph(j, k) ;
305  }
306  }
307  }
308  }
309  killing.std_spectral_base() ;
310 
311  // Check the normalization
312  // -----------------------
313 
314  double check_norm = 0. ;
315  source_phi = pow(confo, 2.) * rr * st / killing(3) ;
316  source_phi.std_spectral_base() ;
317 
318  for (int i=0; i<mm; i++) {
319 
320  t1 = hh * double(4*i) ;
321  t2 = hh * double(4*i+1) ;
322  t3 = hh * double(4*i+2) ;
323  t4 = hh * double(4*i+3) ;
324  t5 = hh * double(4*i+4) ;
325 
326  check_norm += (hh/45.) *
327  ( 14.*source_phi.val_point(rah,M_PI/4.,t1)
328  + 64.*source_phi.val_point(rah,M_PI/4.,t2)
329  + 24.*source_phi.val_point(rah,M_PI/4.,t3)
330  + 64.*source_phi.val_point(rah,M_PI/4.,t4)
331  + 14.*source_phi.val_point(rah,M_PI/4.,t5) ) ;
332 
333  }
334 
335  cout << "Black_hole:: t_f for M_PI/4 = " << check_norm / M_PI
336  << " M_PI" << endl ;
337 
338  } // End of the loop for isotropic coordinates
339 
340  return killing ;
341 
342 }
343 }
Map & mp
Mapping associated with the black hole.
Definition: blackhole.h:80
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
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
Tbl runge_kutta_theta_bh(const Tbl &xi_i, const double &theta_i, const double &phi, const int &nrk) const
Compute a forth-order Runge-Kutta integration to the theta direction for the solution of the Killing ...
Tensor field of valence 1.
Definition: vector.h:188
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:643
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Scalar confo
Conformal factor generated by the black hole.
Definition: blackhole.h:118
Coord sint
Definition: map.h:739
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
Definition: tbl.h:420
virtual double rad_ah() const
Radius of the apparent horizon.
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:423
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
Definition: scalar.C:896
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Multi-domain grid.
Definition: grilles.h:279
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition: scalar.h:690
Vector killing_vect_bh(const Tbl &xi_i, const double &phi_i, const double &theta_i, const int &nrk_phi, const int &nrk_theta) const
Compute the Killing vector of a black hole normalized so that its affine length is 2 M_PI...
Basic array class.
Definition: tbl.h:164
Tbl runge_kutta_phi_bh(const Tbl &xi_i, const double &phi_i, const int &nrk) const
Compute a forth-order Runge-Kutta integration to the phi direction for the solution of the Killing ve...
Coord r
r coordinate centered on the grid
Definition: map.h:736