LORENE
hole_bhns_killing.C
1 /*
2  * Methods of class Hole_bhns to compute Killing vectors
3  *
4  * (see file hole_bhns.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2006-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 hole_bhns_killing_C[] = "$Header: /cvsroot/Lorene/C++/Source/Hole_bhns/hole_bhns_killing.C,v 1.4 2014/10/13 08:53:00 j_novak Exp $" ;
29 
30 /*
31  * $Id: hole_bhns_killing.C,v 1.4 2014/10/13 08:53:00 j_novak Exp $
32  * $Log: hole_bhns_killing.C,v $
33  * Revision 1.4 2014/10/13 08:53:00 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.3 2014/10/06 15:13:10 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.2 2008/07/02 21:01:11 k_taniguchi
40  * A bug removed.
41  *
42  * Revision 1.1 2008/05/15 19:09:53 k_taniguchi
43  * *** empty log message ***
44  *
45  *
46  * $Header: /cvsroot/Lorene/C++/Source/Hole_bhns/hole_bhns_killing.C,v 1.4 2014/10/13 08:53:00 j_novak Exp $
47  *
48  */
49 
50 // C++ headers
51 //#include <>
52 
53 // C headers
54 #include <cmath>
55 
56 // Lorene headers
57 #include "hole_bhns.h"
58 #include "unites.h"
59 #include "utilitaires.h"
60 
61  //---------------------------------------------//
62  // Killing vectors on the AH //
63  //---------------------------------------------//
64 
65 namespace Lorene {
66 Vector Hole_bhns::killing_vect(const Tbl& xi_i, const double& phi_i,
67  const double& theta_i, const int& nrk_phi,
68  const int& nrk_theta) const {
69 
70  using namespace Unites ;
71 
72  assert(xi_i.get_ndim() == 1) ;
73  assert(xi_i.get_dim(0) == 3) ;
74 
75  const Mg3d* mg = mp.get_mg() ;
76  int nr = mg->get_nr(1) ;
77  int nt = mg->get_nt(1) ;
78  int np = mg->get_np(1) ;
79 
80  // Vector which is returned to the main code
81  // Spherical basis, covariant
82  Vector killing(mp, COV, mp.get_bvect_spher()) ;
83 
84  if (kerrschild) {
85 
86  cout << "Not yet prepared!!!" << endl ;
87  abort() ;
88 
89  }
90  else { // Isotropic coordinates
91 
92  // Solution of the Killing vector on the equator
93  // ---------------------------------------------
94 
95  double dp = 2. * M_PI / double(np) ; // Angular step
96 
97  // Killing vector on the equator
98  // np+1 is for the check of xi(phi=0)= xi(phi=2pi)
99  Tbl xi_t(np+1) ;
100  xi_t.set_etat_qcq() ;
101  Tbl xi_p(np+1) ;
102  xi_p.set_etat_qcq() ;
103  Tbl xi_l(np+1) ;
104  xi_l.set_etat_qcq() ;
105 
106  xi_t.set(0) = xi_i(0) ;
107  xi_p.set(0) = xi_i(1) ;
108  xi_l.set(0) = xi_i(2) ;
109 
110  Tbl xi(3) ;
111  xi.set_etat_qcq() ;
112 
113  Tbl xi_ini(3) ;
114  xi_ini.set_etat_qcq() ;
115  xi_ini.set(0) = xi_i(0) ;
116  xi_ini.set(1) = xi_i(1) ;
117  xi_ini.set(2) = xi_i(2) ;
118 
119  double pp_0 = phi_i ; // azimuthal angle phi
120 
121  for (int i=1; i<np+1; i++) {
122 
123  xi = runge_kutta_phi(xi_ini, pp_0, nrk_phi) ;
124 
125  xi_t.set(i) = xi(0) ;
126  xi_p.set(i) = xi(1) ;
127  xi_l.set(i) = xi(2) ;
128 
129  // New data for the next step
130  // -------------------------
131  pp_0 += dp ; // New angle
132  xi_ini = xi ;
133 
134  }
135 
136  cout << "Hole_bhns::killing_vect :" << endl ;
137  cout.precision(16) ;
138  cout << " xi_p(0) : " << xi_p(0) << endl ;
139  cout << " xi_p(np) : " << xi_p(np) << endl ;
140 
141  // Normalization of the Killing vector
142  // -----------------------------------
143 
144  // Initialization of the Killing vector to the phi direction
145  Scalar xi_phi(mp) ;
146  xi_phi = 0.5 ; // If we use "1." for the initialization value,
147  // the state flag becomes "ETATUN" which does not
148  // work for "set_grid_point".
149 
150  for (int k=0; k<np; k++) {
151  xi_phi.set_grid_point(0, k, nt-1, nr-1) = xi_p(k) ;
152  xi_phi.set_grid_point(1, k, nt-1, 0) = xi_p(k) ;
153  }
154  xi_phi.std_spectral_base() ;
155 
156  // Normalization of the Killing vector
157  Scalar rr(mp) ;
158  rr = mp.r ;
159  rr.std_spectral_base() ;
160 
161  Scalar st(mp) ;
162  st = mp.sint ;
163  st.std_spectral_base() ;
164 
165  Scalar source_phi(mp) ;
166  source_phi = pow(confo_tot, 2.) * rr * st / xi_phi ;
167  source_phi.std_spectral_base() ;
168 
169  double rah = rad_ah() ;
170 
171  int nn = 1000 ;
172  double hh = 2. * M_PI / double(nn) ;
173  double integ = 0. ;
174 
175  int mm ;
176  double t1, t2, t3, t4, t5 ;
177 
178  // Boole's Rule (Newton-Cotes Integral) for integration
179  // ----------------------------------------------------
180 
181  assert(nn%4 == 0) ;
182  mm = nn/4 ;
183 
184  for (int i=0; i<mm; i++) {
185 
186  t1 = hh * double(4*i) ;
187  t2 = hh * double(4*i+1) ;
188  t3 = hh * double(4*i+2) ;
189  t4 = hh * double(4*i+3) ;
190  t5 = hh * double(4*i+4) ;
191 
192  integ += (hh/45.) * (14.*source_phi.val_point(rah,M_PI/2.,t1)
193  + 64.*source_phi.val_point(rah,M_PI/2.,t2)
194  + 24.*source_phi.val_point(rah,M_PI/2.,t3)
195  + 64.*source_phi.val_point(rah,M_PI/2.,t4)
196  + 14.*source_phi.val_point(rah,M_PI/2.,t5)
197  ) ;
198 
199  }
200 
201  cout << "Hole_bhns:: t_f = " << integ << endl ;
202  double ratio = 0.5 * integ / M_PI ;
203 
204  cout << "Hole_bhns:: t_f / 2M_PI = " << ratio << endl ;
205 
206  for (int k=0; k<np; k++) {
207  xi_p.set(k) = xi_phi.val_grid_point(1, k, nt-1, 0) * ratio ;
208  }
209 
210 
211  // Solution of the Killing vector to the pole angle
212  // ------------------------------------------------
213 
214  double dt = 0.5 * M_PI / double(nt-1) ; // Angular step
215 
216  // Killing vector to the polar angle
217  Tbl xi_th(nt, np) ;
218  xi_th.set_etat_qcq() ;
219  Tbl xi_ph(nt, np) ;
220  xi_ph.set_etat_qcq() ;
221  Tbl xi_ll(nt, np) ;
222  xi_ll.set_etat_qcq() ;
223 
224  // Values on the equator
225  for (int i=0; i<np; i++) {
226 
227  xi_th.set(nt-1, i) = xi_t(i) ;
228  xi_ph.set(nt-1, i) = xi_p(i) ;
229  xi_ll.set(nt-1, i) = xi_l(i) ;
230 
231  }
232 
233  for (int i=0; i<np; i++) {
234 
235  // Value at theta=pi/2, phi=phi_i
236  xi_ini.set(0) = xi_t(i) ;
237  xi_ini.set(1) = xi_p(i) ;
238  xi_ini.set(2) = xi_l(i) ;
239 
240  double pp = double(i) * dp ;
241  double tt_0 = theta_i ; // polar angle theta
242 
243  for (int j=1; j<nt; j++) {
244 
245  xi = runge_kutta_theta(xi_ini, tt_0, pp, nrk_theta) ;
246 
247  xi_th.set(nt-1-j, i) = xi(0) ;
248  xi_ph.set(nt-1-j, i) = xi(1) ;
249  xi_ll.set(nt-1-j, i) = xi(2) ;
250 
251  // New data for the nxt step
252  // -------------------------
253  tt_0 -= dt ; // New angle
254  xi_ini = xi ;
255 
256  } // End of the loop to the polar direction
257 
258  } // End of the loop to the azimuhtal direction
259 
260  cout << "Hole_bhns::killing_vect :" << endl ;
261  cout.precision(16) ;
262  cout << " xi_ph(nt-1,0) : " << xi_ph(nt-1,0) << endl ;
263  cout << " xi_ph(0,0) : " << xi_ph(0,0) << endl ;
264  cout << " xi_ph(0,np/4) : " << xi_ph(0,np/4) << endl ;
265  cout << " xi_ph(0,np/2) : " << xi_ph(0,np/2) << endl ;
266  cout << " xi_ph(0,3np/4) : " << xi_ph(0,3*np/4) << endl ;
267 
268 
269  // Construction of the Killing vector
270  // ----------------------------------
271 
272  // Definition of the vector is at the top of this code
273  killing.set_etat_qcq() ;
274  killing.set(1) = 0.5 ; // initialization of L instead of
275  // the r component
276  killing.set(2) = 0.5 ; // initialization of the theta component
277  killing.set(3) = 0.5 ; // initialization of the phi component
278 
279  for (int l=0; l<2; l++) {
280  for (int i=0; i<nr; i++) {
281  for (int j=0; j<nt; j++) {
282  for (int k=0; k<np; k++) {
283  (killing.set(1)).set_grid_point(l, k, j, i) = xi_ll(j, k) ;
284  (killing.set(2)).set_grid_point(l, k, j, i) = xi_th(j, k) ;
285  (killing.set(3)).set_grid_point(l, k, j, i) = xi_ph(j, k) ;
286  }
287  }
288  }
289  }
290  killing.std_spectral_base() ;
291 
292  // Check the normalization
293  // -----------------------
294 
295  double check_norm1 = 0. ;
296  double check_norm2 = 0. ;
297  source_phi = pow(confo_tot, 2.) * rr * st / killing(3) ;
298  source_phi.std_spectral_base() ;
299 
300  for (int i=0; i<mm; i++) {
301 
302  t1 = hh * double(4*i) ;
303  t2 = hh * double(4*i+1) ;
304  t3 = hh * double(4*i+2) ;
305  t4 = hh * double(4*i+3) ;
306  t5 = hh * double(4*i+4) ;
307 
308  check_norm1 += (hh/45.) *
309  ( 14.*source_phi.val_point(rah,M_PI/4.,t1)
310  + 64.*source_phi.val_point(rah,M_PI/4.,t2)
311  + 24.*source_phi.val_point(rah,M_PI/4.,t3)
312  + 64.*source_phi.val_point(rah,M_PI/4.,t4)
313  + 14.*source_phi.val_point(rah,M_PI/4.,t5) ) ;
314 
315  check_norm2 += (hh/45.) *
316  ( 14.*source_phi.val_point(rah,M_PI/8.,t1)
317  + 64.*source_phi.val_point(rah,M_PI/8.,t2)
318  + 24.*source_phi.val_point(rah,M_PI/8.,t3)
319  + 64.*source_phi.val_point(rah,M_PI/8.,t4)
320  + 14.*source_phi.val_point(rah,M_PI/8.,t5) ) ;
321 
322  }
323 
324  cout.precision(16) ;
325  cout << "Hole_bhns:: t_f for M_PI/4 = " << check_norm1 / M_PI
326  << " M_PI" << endl ;
327  cout << "Hole_bhns:: t_f for M_PI/8 = " << check_norm2 / M_PI
328  << " M_PI" << endl ;
329 
330  // Checking the accuracy of the Killing vector
331  // -------------------------------------------
332 
333  // xi^i D_i L = 0
334  Scalar dldt(mp) ;
335  dldt = killing(1).dsdt() ;
336 
337  Scalar dldp(mp) ;
338  dldp = killing(1).stdsdp() ;
339 
340  Scalar xidl(mp) ;
341  xidl = killing(2) * dldt + killing(3) * dldp ;
342  xidl.std_spectral_base() ;
343 
344  double xidl_error1 = 0. ;
345  double xidl_error2 = 0. ;
346  double xidl_norm = 0. ;
347 
348  for (int j=0; j<nt; j++) {
349  for (int k=0; k<np/2; k++) {
350  xidl_error1 += xidl.val_grid_point(1, k, j, 0) ;
351  }
352  }
353 
354  for (int j=0; j<nt; j++) {
355  for (int k=np/2; k<np; k++) {
356  xidl_error2 += xidl.val_grid_point(1, k, j, 0) ;
357  }
358  }
359 
360  for (int j=0; j<nt; j++) {
361  for (int k=0; k<np; k++) {
362  xidl_norm += fabs(xidl.val_grid_point(1, k, j, 0)) ;
363  }
364  }
365 
366  cout.precision(6) ;
367  cout << "Relative error in the 1st condition : "
368  << xidl_error1 / xidl_norm / double(nt) / double(np) * 2.
369  << " "
370  << xidl_error2 / xidl_norm / double(nt) / double(np) * 2.
371  << " "
372  << (xidl_error1+xidl_error2)/xidl_norm/double(nt)/double(np)
373  << endl ;
374 
375  // D_i xi^i = 0
376  Scalar xitst(mp) ;
377  xitst = pow(confo_tot, 2.) * rr * st * killing(2) ;
378  xitst.std_spectral_base() ;
379 
380  Scalar dxitstdt(mp) ;
381  dxitstdt = xitst.dsdt() ;
382 
383  Scalar xip(mp) ;
384  xip = pow(confo_tot, 2.) * rr * killing(3) ;
385  xip.std_spectral_base() ;
386 
387  Scalar dxipdp(mp) ;
388  dxipdp = xip.stdsdp() ;
389 
390  Scalar dxi(mp) ;
391  dxi = dxitstdt + st * dxipdp ;
392  dxi.std_spectral_base() ;
393 
394  double dxi_error1 = 0. ;
395  double dxi_error2 = 0. ;
396  double dxi_norm = 0. ;
397 
398  for (int j=0; j<nt; j++) {
399  for (int k=0; k<np/2; k++) {
400  dxi_error1 += dxi.val_grid_point(1, k, j, 0) ;
401  }
402  }
403 
404  for (int j=0; j<nt; j++) {
405  for (int k=np/2; k<np; k++) {
406  dxi_error2 += dxi.val_grid_point(1, k, j, 0) ;
407  }
408  }
409 
410  for (int j=0; j<nt; j++) {
411  for (int k=0; k<np; k++) {
412  dxi_norm += fabs(dxi.val_grid_point(1, k, j, 0)) ;
413  }
414  }
415 
416  cout << "Relative error in the 2nd condition : "
417  << dxi_error1 / dxi_norm / double(nt) / double(np) * 2.
418  << " "
419  << dxi_error2 / dxi_norm / double(nt) / double(np) * 2.
420  << " "
421  << (dxi_error1+dxi_error2)/dxi_norm/double(nt)/double(np)
422  << endl ;
423 
424  // e^{ij} D_i \xi_j = 2L
425  Scalar xipst(mp) ;
426  xipst = pow(confo_tot, 2.) * rr * st * killing(3) ;
427  xipst.std_spectral_base() ;
428 
429  Scalar dxipstdt(mp) ;
430  dxipstdt = xipst.dsdt() ;
431 
432  Scalar xit(mp) ;
433  xit = pow(confo_tot, 2.) * rr * killing(2) ;
434  xit.std_spectral_base() ;
435 
436  Scalar dxitdp(mp) ;
437  dxitdp = xit.stdsdp() ;
438 
439  Scalar dxi2l(mp) ;
440  dxi2l = dxipstdt - st * dxitdp
441  - 2. * killing(1) * pow(confo_tot, 4.) * rr * rr * st ;
442  dxi2l.std_spectral_base() ;
443 
444  double dxi2l_error1 = 0. ;
445  double dxi2l_error2 = 0. ;
446  double dxi2l_norm = 0. ;
447 
448  for (int j=0; j<nt; j++) {
449  for (int k=0; k<np/2; k++) {
450  dxi2l_error1 += dxi2l.val_grid_point(1, k, j, 0) ;
451  }
452  }
453 
454  for (int j=0; j<nt; j++) {
455  for (int k=np/2; k<np; k++) {
456  dxi2l_error2 += dxi2l.val_grid_point(1, k, j, 0) ;
457  }
458  }
459 
460  for (int j=0; j<nt; j++) {
461  for (int k=0; k<np; k++) {
462  dxi2l_norm += fabs(dxi2l.val_grid_point(1, k, j, 0)) ;
463  }
464  }
465 
466  cout << "Relative error in the 3rd condition : "
467  << dxi2l_error1 / dxi2l_norm / double(nt) / double(np) * 2.
468  << " "
469  << dxi2l_error2 / dxi2l_norm / double(nt) / double(np) * 2.
470  << " "
471  << (dxi2l_error1+dxi2l_error2)/dxi2l_norm/double(nt)/double(np)
472  << endl ;
473 
474  cout.precision(4) ;
475 
476  } // End of the loop for isotropic coordinates
477 
478  return killing ;
479 
480 }
481 }
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:783
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.
const Scalar & dsdt() const
Returns of *this .
Definition: scalar_deriv.C:208
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition: blackhole.h:85
Tbl runge_kutta_theta(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 ...
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
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:637
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Coord sint
Definition: map.h:721
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
Definition: tbl.h:400
Vector killing_vect(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...
virtual double rad_ah() const
Radius of the apparent horizon.
const Scalar & stdsdp() const
Returns of *this .
Definition: scalar_deriv.C:238
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:403
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
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:890
Tbl runge_kutta_phi(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...
Scalar confo_tot
Total conformal factor.
Definition: hole_bhns.h:169
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Multi-domain grid.
Definition: grilles.h:273
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:684
Basic array class.
Definition: tbl.h:161
Coord r
r coordinate centered on the grid
Definition: map.h:718