blackhole_r_coord.C
/*
* Method of class Black_hole to express the radial coordinate
* in Kerr-Schild coordinates by that in spatially isotropic coordinates
*
* (see file blackhole.h for documentation).
*
*/
/*
* Copyright (c) 2006-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_r_coord_C[] = "$Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_r_coord.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $" ;
/*
* $Id: blackhole_r_coord.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $
* $Log: blackhole_r_coord.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/05/15 19:29:58 k_taniguchi
* Change of some parameters.
*
* Revision 1.1 2007/06/22 01:20:33 k_taniguchi
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_r_coord.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"
// Local function
namespace Lorene {
double gg(double, const double) ;
const Scalar Black_hole::r_coord(bool neumann, bool first) const {
// Fundamental constants and units
// -------------------------------
using namespace Unites ;
const Mg3d* mg = mp.get_mg() ;
int nz = mg->get_nzone() ; // total number of domains
int nr = mg->get_nr(0) ;
int nt = mg->get_nt(0) ;
int np = mg->get_np(0) ;
double mass = ggrav * mass_bh ;
Scalar r_iso(mp) ;
r_iso = mp.r ;
r_iso.std_spectral_base() ;
Scalar r_are(mp) ;
r_are = r_iso ; // Initialization
r_are.std_spectral_base() ;
// Sets C/M^2 for each case of the lapse boundary condition
// --------------------------------------------------------
double cc ;
if (neumann) { // Neumann boundary condition
if (first) { // First condition
// d(\alpha \psi)/dr = 0
// ---------------------
cc = 2. * (sqrt(13.) - 1.) / 3. ;
}
else { // Second condition
// d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
// -----------------------------------------
cc = 4. / 3. ;
}
}
else { // Dirichlet boundary condition
if (first) { // First condition
// (\alpha \psi) = 1/2
// -------------------
cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
abort() ;
}
else { // Second condition
// (\alpha \psi) = 1/sqrt(2.) \psi_KS
// ----------------------------------
cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
abort() ;
}
}
int ll ;
double diff ;
double ratio ;
double precis = 1.e-15 ;
double dp ;
double tmp ;
double tr ;
int nn = 1000 ;
assert(nn%4 == 0) ;
int mm = nn/4 ;
double x1, x2, x3, x4, x5 ;
double hh, integ ;
// Boole's Rule (Newton-Cotes Integral) for integration
// ----------------------------------------------------
for (int l=1; l<nz; l++) {
for (int i=0; i<nr; i++) {
ratio = 1. ;
dp = 10. ;
tr = r_iso.val_grid_point(l,0,0,i) ;
while ( dp > precis ) {
diff = 1. ; // Initialization
ll = 0 ;
dp = 0.1 * dp ;
while ( diff > precis ) {
ll++ ;
tmp = ratio + ll * dp ;
double r_max = 2.*mass/tmp/tr ;
hh = r_max / double(nn) ;
integ = 0. ;
for (int n=0; n<mm; n++) {
x1 = hh * double(4*n) ;
x2 = hh * double(4*n+1) ;
x3 = hh * double(4*n+2) ;
x4 = hh * double(4*n+3) ;
x5 = hh * double(4*n+4) ;
integ += (hh/45.) * (14.*gg(x1,cc) + 64.*gg(x2,cc)
+ 24.*gg(x3,cc) + 64.*gg(x4,cc)
+ 14.*gg(x5,cc)) ;
}
diff = -log( tmp ) - integ ;
// cout << "diff: " << diff << " x: " << tmp << endl ;
}
ratio += (ll - 1) * dp ;
}
for (int j=0; j<nt; j++) {
for (int k=0; k<np; k++) {
r_are.set_grid_point(l,k,j,i) = ratio ;
}
}
// arrete() ;
}
}
r_are.std_spectral_base() ;
r_are.annule_domain(0) ;
r_are.raccord(1) ;
/*
cout << "r_are:" << endl ;
for (int l=0; l<nz; l++) {
cout << r_are.val_grid_point(l,0,0,0) << " "
<< r_are.val_grid_point(l,0,0,nr-1) << endl ;
}
*/
return r_are ;
}
//*****************************************************************
double gg(double xx, const double cc) {
double tcc2 = cc*cc/16. ;
double tmp = sqrt(1. - xx + tcc2*pow(xx, 4.)) ;
double resu = (-1. + tcc2 * pow(xx, 3.)) / tmp / (1. + tmp) ;
return resu ;
}
}
Lorene::Black_hole::r_coord
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates.
Definition blackhole_r_coord.C:66
Lorene::Black_hole::mp
Map & mp
Mapping associated with the black hole.
Definition blackhole.h:80
Lorene::Black_hole::mass_bh
double mass_bh
Gravitational mass of BH.
Definition blackhole.h:88
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_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
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::raccord
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Definition scalar_raccord.C:62
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::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
Lorene
Lorene prototypes.
Definition app_hor.h:64
Unites
Standard units of space, time and mass.