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
28char 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
65namespace Lorene {
66Vector 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
virtual double rad_ah() const
Radius of the apparent horizon.
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 ...
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.
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
Coord sint
Definition map.h:721
Coord r
r coordinate centered on the grid
Definition map.h:718
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Multi-domain grid.
Definition grilles.h:273
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Scalar & dsdt() const
Returns of *this .
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
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
const Scalar & stdsdp() const
Returns of *this .
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
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
Basic array class.
Definition tbl.h:161
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
Definition tbl.h:400
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition tbl.h:403
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 pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:481
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.