hole_bhns_rk_theta.C
/*
* Methods of class Hole_bhns to compute a forth-order Runge-Kutta
* integration to the theta direction for the solution of the Killing vectors
*
* (see file hole_bhns.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 hole_bhns_rk_theta_C[] = "$Header: /cvsroot/Lorene/C++/Source/Hole_bhns/hole_bhns_rk_theta.C,v 1.4 2014/10/13 08:53:00 j_novak Exp $" ;
/*
* $Id: hole_bhns_rk_theta.C,v 1.4 2014/10/13 08:53:00 j_novak Exp $
* $Log: hole_bhns_rk_theta.C,v $
* Revision 1.4 2014/10/13 08:53:00 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.3 2014/10/06 15:13:10 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.2 2008/07/02 20:49:05 k_taniguchi
* Typos removed.
*
* Revision 1.1 2008/05/15 19:11:01 k_taniguchi
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Hole_bhns/hole_bhns_rk_theta.C,v 1.4 2014/10/13 08:53:00 j_novak Exp $
*
*/
// C++ headers
//#include <>
// C headers
#include <cmath>
// Lorene headers
#include "hole_bhns.h"
#include "unites.h"
#include "utilitaires.h"
//---------------------------------------------------//
// Forth-order Runge-Kutta to the polar angle //
//---------------------------------------------------//
namespace Lorene {
Tbl Hole_bhns::runge_kutta_theta(const Tbl& xi_i, const double& theta_i,
const double& phi,
const int& nrk_theta) const {
using namespace Unites ;
const Mg3d* mg = mp.get_mg() ;
int nt = mg->get_nt(1) ;
Tbl xi_f(3) ; // xi_f(0)=xi_hat{theta}, xi_f(1)=xi_hat{phi}, xi_f(2)=L
xi_f.set_etat_qcq() ;
if (kerrschild) {
cout << "Not yet prepared!!!" << endl ;
abort() ;
}
else { // Isotropic coordinates
// Initial data at phi on the equator
double xi_t0 = xi_i(0) ; // xi_hat{theta}
double xi_p0 = xi_i(1) ; // xi_hat{phi}
double xi_l0 = xi_i(2) ; // L
double theta0 = theta_i ;
double dt = - 0.5 * M_PI / double(nt-1) / double(nrk_theta) ;
// Compute from M_PI/2 to 0
double rah = rad_ah() ;
Scalar dlnconfo(mp) ;
dlnconfo = confo_tot.stdsdp() / confo_tot ;
dlnconfo.std_spectral_base() ;
Scalar laplnconfo(mp) ;
laplnconfo = confo_tot.lapang() / confo_tot ;
laplnconfo.std_spectral_base() ;
Scalar confo2(mp) ;
confo2 = confo_tot * confo_tot ;
confo2.std_spectral_base() ;
double xi_t1, xi_t2, xi_t3, xi_t4, xi_tf ;
double xi_p1, xi_p2, xi_p3, xi_p4, xi_pf ;
double xi_l1, xi_l2, xi_l3, xi_l4, xi_lf ;
double f1, f2, f3, f4 ;
double g1, g2, g3, g4 ;
double h1, h2, h3, h4 ;
// Forth-order Runge-Kutta
// (nrk_theta times steps between two collocation points)
// ------------------------------------------------------
for (int i=0; i<nrk_theta; i++) {
// First
f1 = -2. * xi_p0 * dlnconfo.val_point(rah, theta0, phi) ;
g1 = xi_l0 * rah * confo2.val_point(rah, theta0, phi)
+ 2. * xi_t0 * dlnconfo.val_point(rah, theta0, phi) ;
h1 = - (1. - 2.*laplnconfo.val_point(rah, theta0, phi)) * xi_p0
/ rah / confo2.val_point(rah, theta0, phi) ;
xi_t1 = dt * f1 ;
xi_p1 = dt * g1 ;
xi_l1 = dt * h1 ;
// Second
f2 = -2. * (xi_p0+0.5*xi_p1)
* dlnconfo.val_point(rah, theta0+0.5*dt, phi) ;
g2 = (xi_l0+0.5*xi_l1) * rah
* confo2.val_point(rah, theta0+0.5*dt, phi)
+ 2. * (xi_t0+0.5*xi_t1)
* dlnconfo.val_point(rah, theta0+0.5*dt, phi) ;
h2 = - (1. - 2.*laplnconfo.val_point(rah, theta0+0.5*dt, phi))
* (xi_p0+0.5*xi_p1) / rah
/ confo2.val_point(rah, theta0+0.5*dt, phi) ;
xi_t2 = dt * f2 ;
xi_p2 = dt * g2 ;
xi_l2 = dt * h2 ;
// Third
f3 = -2. * (xi_p0+0.5*xi_p2)
* dlnconfo.val_point(rah, theta0+0.5*dt, phi) ;
g3 = (xi_l0+0.5*xi_l2) * rah
* confo2.val_point(rah, theta0+0.5*dt, phi)
+ 2. * (xi_t0+0.5*xi_t2)
* dlnconfo.val_point(rah, theta0+0.5*dt, phi) ;
h3 = - (1. - 2.*laplnconfo.val_point(rah, theta0+0.5*dt, phi))
* (xi_p0+0.5*xi_p2) / rah
/ confo2.val_point(rah, theta0+0.5*dt, phi) ;
xi_t3 = dt * f3 ;
xi_p3 = dt * g3 ;
xi_l3 = dt * h3 ;
// Forth
f4 = -2. * (xi_p0+xi_p3) * dlnconfo.val_point(rah, theta0+dt, phi) ;
g4 = (xi_l0+xi_l3) * rah * confo2.val_point(rah, theta0+dt, phi)
+ 2. * (xi_t0+xi_t3) * dlnconfo.val_point(rah, theta0+dt, phi) ;
h4 = - (1. - 2.*laplnconfo.val_point(rah, theta0+dt, phi))
* (xi_p0+xi_p3) / rah / confo2.val_point(rah, theta0+dt, phi) ;
xi_t4 = dt * f4 ;
xi_p4 = dt * g4 ;
xi_l4 = dt * h4 ;
// Final results
// -------------
xi_tf = xi_t0 + (xi_t1 + 2.*xi_t2 + 2.*xi_t3 + xi_t4) / 6. ;
xi_pf = xi_p0 + (xi_p1 + 2.*xi_p2 + 2.*xi_p3 + xi_p4) / 6. ;
xi_lf = xi_l0 + (xi_l1 + 2.*xi_l2 + 2.*xi_l3 + xi_l4) / 6. ;
// Final results are put into the initial data
// in order for the next step
// -------------------------------------------
xi_t0 = xi_tf ;
xi_p0 = xi_pf ;
xi_l0 = xi_lf ;
} // End of the loop
xi_f.set(0) = xi_tf ;
xi_f.set(1) = xi_pf ;
xi_f.set(2) = xi_lf ;
}
return xi_f ;
}
}
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::kerrschild
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition blackhole.h:85
Lorene::Hole_bhns::runge_kutta_theta
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 ...
Definition hole_bhns_rk_theta.C:67
Lorene::Hole_bhns::confo_tot
Scalar confo_tot
Total conformal factor.
Definition hole_bhns.h:169
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
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::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_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::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
Lorene prototypes.
Definition app_hor.h:64
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:720
Unites
Standard units of space, time and mass.