60 #include "hole_bhns.h" 62 #include "utilitaires.h" 70 const double& theta_i,
const int& nrk_phi,
71 const int& nrk_theta)
const {
80 int nt = mg->get_nt(1) ;
81 int np = mg->get_np(1) ;
89 cout <<
"Not yet prepared!!!" << endl ;
98 double dp = 2. * M_PI / double(np) ;
109 xi_t.
set(0) = xi_i(0) ;
110 xi_p.
set(0) = xi_i(1) ;
111 xi_l.
set(0) = xi_i(2) ;
118 xi_ini.
set(0) = xi_i(0) ;
119 xi_ini.
set(1) = xi_i(1) ;
120 xi_ini.
set(2) = xi_i(2) ;
122 double pp_0 = phi_i ;
124 for (
int i=1; i<np+1; i++) {
128 xi_t.
set(i) = xi(0) ;
129 xi_p.
set(i) = xi(1) ;
130 xi_l.
set(i) = xi(2) ;
139 cout <<
"Hole_bhns::killing_vect :" << endl ;
141 cout <<
" xi_p(0) : " << xi_p(0) << endl ;
142 cout <<
" xi_p(np) : " << xi_p(np) << endl ;
153 for (
int k=0; k<np; k++) {
175 double hh = 2. * M_PI / double(nn) ;
179 double t1, t2, t3, t4, t5 ;
187 for (
int i=0; i<mm; i++) {
189 t1 = hh * double(4*i) ;
190 t2 = hh * double(4*i+1) ;
191 t3 = hh * double(4*i+2) ;
192 t4 = hh * double(4*i+3) ;
193 t5 = hh * double(4*i+4) ;
195 integ += (hh/45.) * (14.*source_phi.
val_point(rah,M_PI/2.,t1)
196 + 64.*source_phi.
val_point(rah,M_PI/2.,t2)
197 + 24.*source_phi.
val_point(rah,M_PI/2.,t3)
198 + 64.*source_phi.
val_point(rah,M_PI/2.,t4)
199 + 14.*source_phi.
val_point(rah,M_PI/2.,t5)
204 cout <<
"Hole_bhns:: t_f = " << integ << endl ;
205 double ratio = 0.5 * integ / M_PI ;
207 cout <<
"Hole_bhns:: t_f / 2M_PI = " << ratio << endl ;
209 for (
int k=0; k<np; k++) {
217 double dt = 0.5 * M_PI / double(nt-1) ;
228 for (
int i=0; i<np; i++) {
230 xi_th.
set(nt-1, i) = xi_t(i) ;
231 xi_ph.
set(nt-1, i) = xi_p(i) ;
232 xi_ll.
set(nt-1, i) = xi_l(i) ;
236 for (
int i=0; i<np; i++) {
239 xi_ini.
set(0) = xi_t(i) ;
240 xi_ini.
set(1) = xi_p(i) ;
241 xi_ini.
set(2) = xi_l(i) ;
243 double pp = double(i) * dp ;
244 double tt_0 = theta_i ;
246 for (
int j=1; j<nt; j++) {
250 xi_th.
set(nt-1-j, i) = xi(0) ;
251 xi_ph.
set(nt-1-j, i) = xi(1) ;
252 xi_ll.
set(nt-1-j, i) = xi(2) ;
263 cout <<
"Hole_bhns::killing_vect :" << endl ;
265 cout <<
" xi_ph(nt-1,0) : " << xi_ph(nt-1,0) << endl ;
266 cout <<
" xi_ph(0,0) : " << xi_ph(0,0) << endl ;
267 cout <<
" xi_ph(0,np/4) : " << xi_ph(0,np/4) << endl ;
268 cout <<
" xi_ph(0,np/2) : " << xi_ph(0,np/2) << endl ;
269 cout <<
" xi_ph(0,3np/4) : " << xi_ph(0,3*np/4) << endl ;
277 killing.set(1) = 0.5 ;
279 killing.set(2) = 0.5 ;
280 killing.set(3) = 0.5 ;
282 for (
int l=0; l<2; l++) {
283 for (
int i=0; i<nr; i++) {
284 for (
int j=0; j<nt; j++) {
285 for (
int k=0; k<np; k++) {
286 (killing.set(1)).set_grid_point(l, k, j, i) = xi_ll(j, k) ;
287 (killing.set(2)).set_grid_point(l, k, j, i) = xi_th(j, k) ;
288 (killing.set(3)).set_grid_point(l, k, j, i) = xi_ph(j, k) ;
293 killing.std_spectral_base() ;
298 double check_norm1 = 0. ;
299 double check_norm2 = 0. ;
300 source_phi =
pow(
confo_tot, 2.) * rr * st / killing(3) ;
303 for (
int i=0; i<mm; i++) {
305 t1 = hh * double(4*i) ;
306 t2 = hh * double(4*i+1) ;
307 t3 = hh * double(4*i+2) ;
308 t4 = hh * double(4*i+3) ;
309 t5 = hh * double(4*i+4) ;
311 check_norm1 += (hh/45.) *
312 ( 14.*source_phi.
val_point(rah,M_PI/4.,t1)
313 + 64.*source_phi.
val_point(rah,M_PI/4.,t2)
314 + 24.*source_phi.
val_point(rah,M_PI/4.,t3)
315 + 64.*source_phi.
val_point(rah,M_PI/4.,t4)
316 + 14.*source_phi.
val_point(rah,M_PI/4.,t5) ) ;
318 check_norm2 += (hh/45.) *
319 ( 14.*source_phi.
val_point(rah,M_PI/8.,t1)
320 + 64.*source_phi.
val_point(rah,M_PI/8.,t2)
321 + 24.*source_phi.
val_point(rah,M_PI/8.,t3)
322 + 64.*source_phi.
val_point(rah,M_PI/8.,t4)
323 + 14.*source_phi.
val_point(rah,M_PI/8.,t5) ) ;
328 cout <<
"Hole_bhns:: t_f for M_PI/4 = " << check_norm1 / M_PI
330 cout <<
"Hole_bhns:: t_f for M_PI/8 = " << check_norm2 / M_PI
338 dldt = killing(1).
dsdt() ;
341 dldp = killing(1).
stdsdp() ;
344 xidl = killing(2) * dldt + killing(3) * dldp ;
347 double xidl_error1 = 0. ;
348 double xidl_error2 = 0. ;
349 double xidl_norm = 0. ;
351 for (
int j=0; j<nt; j++) {
352 for (
int k=0; k<np/2; k++) {
357 for (
int j=0; j<nt; j++) {
358 for (
int k=np/2; k<np; k++) {
363 for (
int j=0; j<nt; j++) {
364 for (
int k=0; k<np; k++) {
370 cout <<
"Relative error in the 1st condition : " 371 << xidl_error1 / xidl_norm / double(nt) / double(np) * 2.
373 << xidl_error2 / xidl_norm / double(nt) / double(np) * 2.
375 << (xidl_error1+xidl_error2)/xidl_norm/
double(nt)/double(np)
384 dxitstdt = xitst.
dsdt() ;
394 dxi = dxitstdt + st * dxipdp ;
397 double dxi_error1 = 0. ;
398 double dxi_error2 = 0. ;
399 double dxi_norm = 0. ;
401 for (
int j=0; j<nt; j++) {
402 for (
int k=0; k<np/2; k++) {
407 for (
int j=0; j<nt; j++) {
408 for (
int k=np/2; k<np; k++) {
413 for (
int j=0; j<nt; j++) {
414 for (
int k=0; k<np; k++) {
419 cout <<
"Relative error in the 2nd condition : " 420 << dxi_error1 / dxi_norm / double(nt) / double(np) * 2.
422 << dxi_error2 / dxi_norm / double(nt) / double(np) * 2.
424 << (dxi_error1+dxi_error2)/dxi_norm/
double(nt)/double(np)
433 dxipstdt = xipst.
dsdt() ;
443 dxi2l = dxipstdt - st * dxitdp
447 double dxi2l_error1 = 0. ;
448 double dxi2l_error2 = 0. ;
449 double dxi2l_norm = 0. ;
451 for (
int j=0; j<nt; j++) {
452 for (
int k=0; k<np/2; k++) {
457 for (
int j=0; j<nt; j++) {
458 for (
int k=np/2; k<np; k++) {
463 for (
int j=0; j<nt; j++) {
464 for (
int k=0; k<np; k++) {
469 cout <<
"Relative error in the 3rd condition : " 470 << dxi2l_error1 / dxi2l_norm / double(nt) / double(np) * 2.
472 << dxi2l_error2 / dxi2l_norm / double(nt) / double(np) * 2.
474 << (dxi2l_error1+dxi2l_error2)/dxi2l_norm/
double(nt)/double(np)
Map & mp
Mapping associated with the black hole.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Standard units of space, time and mass.
const Scalar & dsdt() const
Returns of *this .
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Tensor field of valence 0 (or component of a tensorial field).
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
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. ...
Tensor field of valence 1.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
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 .
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Cmp pow(const Cmp &, int)
Power .
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.
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.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Coord r
r coordinate centered on the grid