blackhole_killing.C
/*
* Methods of class Black_hole to compute Killing vectors
*
* (see file blackhole.h for documentation).
*
*/
/*
* Copyright (c) 2007 Keisuke Taniguchi
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License version 2
* as published by the Free Software Foundation.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
char blackhole_killing_C[] = "$Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_killing.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $" ;
/*
* $Id: blackhole_killing.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $
* $Log: blackhole_killing.C,v $
* Revision 1.4 2014/10/13 08:52:46 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.3 2014/10/06 15:13:02 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.2 2008/07/02 20:45:07 k_taniguchi
* A bug removed.
*
* Revision 1.1 2008/05/15 19:33:12 k_taniguchi
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_killing.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $
*
*/
// C++ headers
//#include <>
// C headers
#include <cmath>
// Lorene headers
#include "blackhole.h"
#include "unites.h"
#include "utilitaires.h"
//---------------------------------------------//
// Killing vectors on the AH //
//---------------------------------------------//
namespace Lorene {
Vector Black_hole::killing_vect_bh(const Tbl& xi_i, const double& phi_i,
const double& theta_i, const int& nrk_phi,
const int& nrk_theta) const {
using namespace Unites ;
assert(xi_i.get_ndim() == 1) ;
assert(xi_i.get_dim(0) == 3) ;
const Mg3d* mg = mp.get_mg() ;
int nr = mg->get_nr(1) ;
int nt = mg->get_nt(1) ;
int np = mg->get_np(1) ;
// Vector which is returned to the main code
// Spherical basis, covariant
Vector killing(mp, COV, mp.get_bvect_spher()) ;
if (kerrschild) {
cout << "Not yet prepared!!!" << endl ;
abort() ;
}
else { // Isotropic coordinates
// Solution of the Killing vector on the equator
// ---------------------------------------------
double dp = 2. * M_PI / double(np) ; // Angular step
// Killing vector on the equator
// np+1 is for the check of xi(phi=0)= xi(phi=2pi)
Tbl xi_t(np+1) ;
xi_t.set_etat_qcq() ;
Tbl xi_p(np+1) ;
xi_p.set_etat_qcq() ;
Tbl xi_l(np+1) ;
xi_l.set_etat_qcq() ;
xi_t.set(0) = xi_i(0) ;
xi_p.set(0) = xi_i(1) ;
xi_l.set(0) = xi_i(2) ;
Tbl xi(3) ;
xi.set_etat_qcq() ;
Tbl xi_ini(3) ;
xi_ini.set_etat_qcq() ;
xi_ini.set(0) = xi_i(0) ;
xi_ini.set(1) = xi_i(1) ;
xi_ini.set(2) = xi_i(2) ;
double pp_0 = phi_i ; // azimuthal angle phi
for (int i=1; i<np+1; i++) {
xi = runge_kutta_phi_bh(xi_ini, pp_0, nrk_phi) ;
xi_t.set(i) = xi(0) ;
xi_p.set(i) = xi(1) ;
xi_l.set(i) = xi(2) ;
// New data for the next step
// -------------------------
pp_0 += dp ; // New angle
xi_ini = xi ;
}
/*
for (int i=0; i<np+1; i++) {
cout << "xi_t xi_p xi_l" << endl ;
cout << xi_t(i) << " " << xi_p(i) << " " << xi_l(i) << endl ;
}
arrete() ;
*/
// Normalization of the Killing vector
// -----------------------------------
// Initialization of the Killing vector to the phi direction
Scalar xi_phi(mp) ;
xi_phi = 0.5 ; // If we use "1." for the initialization value,
// the state flag becomes "ETATUN" which does not
// work for "set_grid_point".
for (int k=0; k<np; k++) {
xi_phi.set_grid_point(0, k, nt-1, nr-1) = xi_p(k) ;
xi_phi.set_grid_point(1, k, nt-1, 0) = xi_p(k) ;
}
xi_phi.std_spectral_base() ;
/*
for (int l=0; l<2; l++) {
for (int k=0; k<np; k++) {
for (int j=0; j<nt; j++) {
for (int i=0; i<nr; i++) {
cout << "(l,k,j,i)=" << l << "," << k << "," << j
<< "," << i << ": "
<< xi_phi.val_grid_point(l,k,j,i) << endl ;
}
arrete() ;
}
arrete() ;
}
arrete() ;
}
*/
// Normalization of the Killing vector
Scalar rr(mp) ;
rr = mp.r ;
rr.std_spectral_base() ;
Scalar st(mp) ;
st = mp.sint ;
st.std_spectral_base() ;
Scalar source_phi(mp) ;
source_phi = pow(confo, 2.) * rr * st / xi_phi ;
source_phi.std_spectral_base() ;
double rah = rad_ah() ;
int nn = 1000 ;
double hh = 2. * M_PI / double(nn) ;
double integ = 0. ;
int mm ;
double t1, t2, t3, t4, t5 ;
// Boole's Rule (Newton-Cotes Integral) for integration
// ----------------------------------------------------
assert(nn%4 == 0) ;
mm = nn/4 ;
for (int i=0; i<mm; i++) {
t1 = hh * double(4*i) ;
t2 = hh * double(4*i+1) ;
t3 = hh * double(4*i+2) ;
t4 = hh * double(4*i+3) ;
t5 = hh * double(4*i+4) ;
integ += (hh/45.) * (14.*source_phi.val_point(rah,M_PI/2.,t1)
+ 64.*source_phi.val_point(rah,M_PI/2.,t2)
+ 24.*source_phi.val_point(rah,M_PI/2.,t3)
+ 64.*source_phi.val_point(rah,M_PI/2.,t4)
+ 14.*source_phi.val_point(rah,M_PI/2.,t5)
) ;
}
cout << "Black_hole:: t_f = " << integ << endl ;
double ratio = 0.5 * integ / M_PI ;
cout << "Black_hole:: t_f / 2M_PI = " << ratio << endl ;
for (int k=0; k<np; k++) {
xi_p.set(k) = xi_phi.val_grid_point(1, k, nt-1, 0) * ratio ;
}
/*
for (int k=0; k<np; k++) {
cout << "Normalized xi_p" << "(" << k << ") :" << xi_p(k) << endl ;
}
*/
// Solution of the Killing vector to the pole angle
// ------------------------------------------------
double dt = 0.5 * M_PI / double(nt-1) ; // Angular step
// Killing vector to the polar angle
Tbl xi_th(nt, np) ;
xi_th.set_etat_qcq() ;
Tbl xi_ph(nt, np) ;
xi_ph.set_etat_qcq() ;
Tbl xi_ll(nt, np) ;
xi_ll.set_etat_qcq() ;
// Values on the equator
for (int i=0; i<np; i++) {
xi_th.set(nt-1, i) = xi_t(i) ;
xi_ph.set(nt-1, i) = xi_p(i) ;
xi_ll.set(nt-1, i) = xi_l(i) ;
}
for (int i=0; i<np; i++) {
// Value at theta=pi/2, phi=phi_i
xi_ini.set(0) = xi_t(i) ;
xi_ini.set(1) = xi_p(i) ;
xi_ini.set(2) = xi_l(i) ;
double pp = double(i) * dp ;
double tt_0 = theta_i ; // polar angle theta
for (int j=1; j<nt; j++) {
xi = runge_kutta_theta_bh(xi_ini, tt_0, pp, nrk_theta) ;
xi_th.set(nt-1-j, i) = xi(0) ;
xi_ph.set(nt-1-j, i) = xi(1) ;
xi_ll.set(nt-1-j, i) = xi(2) ;
// New data for the nxt step
// -------------------------
tt_0 -= dt ; // New angle
xi_ini = xi ;
} // End of the loop to the polar direction
} // End of the loop to the azimuhtal direction
// Construction of the Killing vector
// ----------------------------------
// Definition of the vector is at the top of this code
killing.set_etat_qcq() ;
killing.set(1) = 0. ; // r component
killing.set(2) = 0.5 ; // initialization of the theta component
killing.set(3) = 0.5 ; // initialization of the phi component
for (int l=0; l<2; l++) {
for (int i=0; i<nr; i++) {
for (int j=0; j<nt; j++) {
for (int k=0; k<np; k++) {
(killing.set(2)).set_grid_point(l, k, j, i) = xi_th(j, k) ;
(killing.set(3)).set_grid_point(l, k, j, i) = xi_ph(j, k) ;
}
}
}
}
killing.std_spectral_base() ;
// Check the normalization
// -----------------------
double check_norm = 0. ;
source_phi = pow(confo, 2.) * rr * st / killing(3) ;
source_phi.std_spectral_base() ;
for (int i=0; i<mm; i++) {
t1 = hh * double(4*i) ;
t2 = hh * double(4*i+1) ;
t3 = hh * double(4*i+2) ;
t4 = hh * double(4*i+3) ;
t5 = hh * double(4*i+4) ;
check_norm += (hh/45.) *
( 14.*source_phi.val_point(rah,M_PI/4.,t1)
+ 64.*source_phi.val_point(rah,M_PI/4.,t2)
+ 24.*source_phi.val_point(rah,M_PI/4.,t3)
+ 64.*source_phi.val_point(rah,M_PI/4.,t4)
+ 14.*source_phi.val_point(rah,M_PI/4.,t5) ) ;
}
cout << "Black_hole:: t_f for M_PI/4 = " << check_norm / M_PI
<< " M_PI" << endl ;
} // End of the loop for isotropic coordinates
return killing ;
}
}
Lorene::Black_hole::runge_kutta_phi_bh
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...
Definition blackhole_rk_phi.C:67
Lorene::Black_hole::mp
Map & mp
Mapping associated with the black hole.
Definition blackhole.h:80
Lorene::Black_hole::rad_ah
virtual double rad_ah() const
Radius of the apparent horizon.
Definition blackhole_global.C:505
Lorene::Black_hole::runge_kutta_theta_bh
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 ...
Definition blackhole_rk_theta.C:67
Lorene::Black_hole::killing_vect_bh
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.
Definition blackhole_killing.C:66
Lorene::Black_hole::kerrschild
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition blackhole.h:85
Lorene::Black_hole::confo
Scalar confo
Conformal factor generated by the black hole.
Definition blackhole.h:118
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
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::Scalar::val_grid_point
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
Lorene::Scalar::set_grid_point
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
Lorene::Scalar::val_point
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
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::get_ndim
int get_ndim() const
Gives the number of dimensions (ie dim.ndim).
Definition tbl.h:400
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Tbl::get_dim
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i]).
Definition tbl.h:403
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
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
Lorene
Lorene prototypes.
Definition app_hor.h:64
Unites
Standard units of space, time and mass.