blackhole_extr_curv.C
/*
* Method of class Black_hole to compute the extrinsic curvature tensor
*
* (see file blackhole.h for documentation).
*
*/
/*
* Copyright (c) 2005-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_extr_curv_C[] = "$Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_extr_curv.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $" ;
/*
* $Id: blackhole_extr_curv.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $
* $Log: blackhole_extr_curv.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:27:14 k_taniguchi
* Change of some parameters.
*
* Revision 1.1 2007/06/22 01:19:32 k_taniguchi
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_extr_curv.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 "utilitaires.h"
#include "unites.h"
namespace Lorene {
void Black_hole::extr_curv_bh() {
// Fundamental constants and units
// -------------------------------
using namespace Unites ;
if (kerrschild) {
double mass = ggrav * mass_bh ;
Scalar rr(mp) ;
rr = mp.r ;
rr.std_spectral_base() ;
Scalar st(mp) ;
st = mp.sint ;
st.std_spectral_base() ;
Scalar ct(mp) ;
ct = mp.cost ;
ct.std_spectral_base() ;
Scalar sp(mp) ;
sp = mp.sinp ;
sp.std_spectral_base() ;
Scalar cp(mp) ;
cp = mp.cosp ;
cp.std_spectral_base() ;
Vector ll(mp, CON, mp.get_bvect_cart()) ;
ll.set_etat_qcq() ;
ll.set(1) = st * cp ;
ll.set(2) = st * sp ;
ll.set(3) = ct ;
ll.std_spectral_base() ;
// Computation of \tilde{A}^{ij}
// -----------------------------
Scalar divshift(mp) ;
divshift = shift_rs(1).dsdx() + shift_rs(2).dsdy()
+ shift_rs(3).dsdz() ;
divshift.std_spectral_base() ;
Sym_tensor flat_taij(mp, CON, mp.get_bvect_cart()) ;
flat_taij.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
flat_taij.set(i,j) = shift_rs(j).deriv(i)
+ shift_rs(i).deriv(j)
- 2. * divshift * flat.con()(i,j) / 3. ;
}
}
flat_taij.std_spectral_base() ;
Scalar lapse_bh(mp) ;
lapse_bh = 1. / sqrt(1. + 2. * mass / rr) ;
lapse_bh.std_spectral_base() ;
lapse_bh.annule_domain(0) ;
lapse_bh.raccord(1) ;
Sym_tensor curv_taij(mp, CON, mp.get_bvect_cart()) ;
curv_taij.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
curv_taij.set(i,j) = -2. * lapse_bh * lapse_bh * mass
* (ll(i)*(shift_rs(j).dsdr()) + ll(j)*(shift_rs(i).dsdr())
- 2. * ll(i) * ll(j) * divshift / 3.) / rr ;
}
}
curv_taij.std_spectral_base() ;
Sym_tensor resi_taij(mp, CON, mp.get_bvect_cart()) ;
resi_taij.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
resi_taij.set(i,j) = 2. * lapse_bh * lapse_bh * mass
* ( ll(i) * shift_rs(j) + ll(j) * shift_rs(i)
+ ( flat.con()(i,j)
- lapse_bh*lapse_bh*(9.+14.*mass/rr)*ll(i)*ll(j) )
* ( ll(1)*shift_rs(1)+ll(2)*shift_rs(2)
+ll(3)*shift_rs(3) )/ 3. )
/ rr / rr ;
}
}
resi_taij.std_spectral_base() ;
resi_taij.inc_dzpuis(2) ;
taij_rs = 0.5 * pow(confo, 7.)
* (flat_taij + curv_taij + resi_taij) / lapconf ;
taij_rs.std_spectral_base() ;
taij_rs.annule_domain(0) ;
Sym_tensor taij_bh(mp, CON, mp.get_bvect_cart()) ;
taij_bh.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
taij_bh.set(i,j) = 2.*pow(lapse_bh,6.)*mass*(2.+3.*mass/rr)
*( (1.+2.*mass/rr) * flat.con()(i,j)
- (3.+2.*mass/rr) * ll(i) * ll(j) )
*pow(confo, 7.)/lapconf/3./rr/rr ;
}
}
taij_bh.std_spectral_base() ;
taij_bh.inc_dzpuis(2) ;
taij_bh.annule_domain(0) ;
taij = taij_rs + taij_bh ;
taij.std_spectral_base() ;
taij.annule_domain(0) ;
/*
Sym_tensor taij_ks_con(mp, CON, mp.get_bvect_cart()) ;
taij_ks_con.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
taij_ks_con.set(i,j) = 2.*pow(lap_bh2,2.5)*mass
* (2.+3.*mass/rr)
* ( (1.+2.*mass/rr)*flat.con()(i,j)
- (3.+2.*mass/rr)*ll(i)*ll(j) ) / 3. / rr / rr ;
}
}
taij_ks_con.std_spectral_base() ;
taij_ks_con.annule_domain(0) ;
taij_ks_con.inc_dzpuis(2) ;
cout << "taij(1,1) - taij_ks_con(1,1) : " << endl ;
cout << taij(1,1) - taij_ks_con(1,1) << endl ;
arrete() ;
cout << "taij_ks_con(1,1) : " << endl ;
cout << taij_ks_con(1,1) << endl ;
arrete() ;
cout << "taij(1,1) : " << endl ;
cout << taij(1,1) << endl ;
arrete() ;
*/
// Computation of \tilde{A}^{ij} \tilde{A}_{ij}
// --------------------------------------------
Sym_tensor flat_dshift(mp, COV, mp.get_bvect_cart()) ;
flat_dshift.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
flat_dshift.set(i,j) = flat.cov()(j,1)*(shift_rs(1).deriv(i))
+ flat.cov()(j,2)*(shift_rs(2).deriv(i))
+ flat.cov()(j,3)*(shift_rs(3).deriv(i))
+ flat.cov()(i,1)*(shift_rs(1).deriv(j))
+ flat.cov()(i,2)*(shift_rs(2).deriv(j))
+ flat.cov()(i,3)*(shift_rs(3).deriv(j))
- 2. * divshift * flat.cov()(i,j) / 3. ;
}
}
flat_dshift.std_spectral_base() ;
Sym_tensor curv_dshift(mp, COV, mp.get_bvect_cart()) ;
curv_dshift.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
curv_dshift.set(i,j) = 2. * mass
*( ll(j) *( ll(1)*(shift_rs(1).deriv(i))
+ ll(2)*(shift_rs(2).deriv(i))
+ ll(3)*(shift_rs(3).deriv(i)) )
+ ll(i) *( ll(1)*(shift_rs(1).deriv(j))
+ ll(2)*(shift_rs(2).deriv(j))
+ ll(3)*(shift_rs(3).deriv(j)) )
- 2. * divshift * ll(i) * ll(j) / 3. ) / rr ;
}
}
curv_dshift.std_spectral_base() ;
Sym_tensor tmp1(mp, COV, mp.get_bvect_cart()) ;
tmp1.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
tmp1.set(i,j) = 2. * mass
*(ll(j)*( (flat.cov()(i,1)+2.*mass*ll(i)*ll(1)/rr)
*shift_rs(1)
+ (flat.cov()(i,2)+2.*mass*ll(i)*ll(2)/rr)
*shift_rs(2)
+ (flat.cov()(i,3)+2.*mass*ll(i)*ll(3)/rr)
*shift_rs(3)
)
+ ll(i)*( (flat.cov()(j,1)+2.*mass*ll(j)*ll(1)/rr)
*shift_rs(1)
+ (flat.cov()(j,2)+2.*mass*ll(j)*ll(2)/rr)
*shift_rs(2)
+ (flat.cov()(j,3)+2.*mass*ll(j)*ll(3)/rr)
*shift_rs(3) )
) / rr / rr ;
}
}
tmp1.std_spectral_base() ;
tmp1.inc_dzpuis(2) ;
Sym_tensor tmp2(mp, COV, mp.get_bvect_cart()) ;
tmp2.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
tmp2.set(i,j) = 2. * mass * lapse_bh * lapse_bh
* (ll(1)*shift_rs(1)+ll(2)*shift_rs(2)+ll(3)*shift_rs(3))
* (flat.cov()(i,j)
- (9.+28.*mass/rr+24.*mass*mass/rr/rr)*ll(i)*ll(j))
/ 3. / rr / rr ;
}
}
tmp2.std_spectral_base() ;
tmp2.inc_dzpuis(2) ;
Sym_tensor taij_rs_down(mp, COV, mp.get_bvect_cart()) ;
taij_rs_down.set_etat_qcq() ;
taij_rs_down = 0.5 * pow(confo, 7.)
* (flat_dshift + curv_dshift + tmp1 + tmp2) / lapconf ;
taij_rs_down.std_spectral_base() ;
taij_rs_down.annule_domain(0) ;
Sym_tensor taij_bh_down(mp, COV, mp.get_bvect_cart()) ;
taij_bh_down.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
taij_bh_down.set(i,j) = 2.*pow(lapse_bh,4.)*mass*(2.+3.*mass/rr)
*pow(confo,7.)*(flat.cov()(i,j)-(3.+4.*mass/rr)*ll(i)*ll(j))
/lapconf/3./rr/rr ;
}
}
taij_bh_down.std_spectral_base() ;
taij_bh_down.inc_dzpuis(2) ;
taij_bh_down.annule_domain(0) ;
/*
Sym_tensor taij_ks_cov(mp, COV, mp.get_bvect_cart()) ;
taij_ks_cov.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
taij_ks_cov.set(i,j) = 2.*pow(lap_bh2,1.5)*mass * (2.+3.*mass/rr)
* ( flat.cov()(i,j)
- (3.+4.*mass/rr)*ll(i)*ll(j) ) / 3. / rr / rr ;
}
}
taij_ks_cov.std_spectral_base() ;
taij_ks_cov.annule_domain(0) ;
taij_ks_cov.inc_dzpuis(2) ;
cout << "taij_down(1,1) - taij_ks_cov(1,1) : " << endl ;
cout << taij_down(1,1) - taij_ks_cov(1,1) << endl ;
arrete() ;
cout << "taij_ks_cov(1,1) : " << endl ;
cout << taij_ks_cov(1,1) << endl ;
arrete() ;
cout << "taij_down(1,1) : " << endl ;
cout << taij_down(1,1) << endl ;
arrete() ;
*/
Scalar taij_quad_rsrs(mp) ;
taij_quad_rsrs = 0. ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
taij_quad_rsrs += taij_rs_down(i,j) * taij_rs(i,j) ;
}
}
taij_quad_rsrs.std_spectral_base() ;
Scalar taij_quad_rsbh1(mp) ;
taij_quad_rsbh1 = 0. ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
taij_quad_rsbh1 += taij_rs_down(i,j) * taij_bh(i,j) ;
}
}
taij_quad_rsbh1.std_spectral_base() ;
Scalar taij_quad_rsbh2(mp) ;
taij_quad_rsbh2 = 0. ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
taij_quad_rsbh2 += taij_bh_down(i,j) * taij_rs(i,j) ;
}
}
taij_quad_rsbh2.std_spectral_base() ;
taij_quad_rs = taij_quad_rsrs + taij_quad_rsbh1 + taij_quad_rsbh2 ;
taij_quad_rs.std_spectral_base() ;
Scalar taij_quad_bh(mp) ;
taij_quad_bh = 8.*pow(lapse_bh,10.)*mass*mass*(2.+3.*mass/rr)
*(2.+3.*mass/rr)*pow(confo,12.)/3./pow(rr,4.)/lapconf/lapconf ;
taij_quad_bh.std_spectral_base() ;
taij_quad_bh.inc_dzpuis(4) ;
taij_quad = taij_quad_rs + taij_quad_bh ;
taij_quad.std_spectral_base() ;
/*
Scalar taij_quad_ks(mp) ;
taij_quad_ks = 8. * pow(lap_bh2,3.) * mass * mass * (2.+3.*mass/rr)
* (2.+3.*mass/rr) / 3. / pow(rr, 4.) ;
taij_quad_ks.std_spectral_base() ;
taij_quad_ks.annule_domain(0) ;
taij_quad_ks.inc_dzpuis(4) ;
cout << "taij_quad - taij_quad_ks : " << endl ;
cout << taij_quad - taij_quad_ks << endl ;
arrete() ;
*/
}
else { // Isotropic coordinates with the maximal slicing
// Computation of \tilde{A}^{ij}
// -----------------------------
Scalar divshift(mp) ;
divshift = shift(1).dsdx() + shift(2).dsdy() + shift(3).dsdz() ;
divshift.std_spectral_base() ;
Sym_tensor flat_taij(mp, CON, mp.get_bvect_cart()) ;
flat_taij.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
flat_taij.set(i,j) = shift(j).deriv(i) + shift(i).deriv(j)
- 2. * divshift * flat.con()(i,j) / 3. ;
}
}
flat_taij.std_spectral_base() ;
taij = 0.5 * pow(confo, 7.) * flat_taij / lapconf ;
taij.std_spectral_base() ;
taij.annule_domain(0) ;
// Computation of \tilde{A}^{ij} \tilde{A}_{ij}
// --------------------------------------------
Sym_tensor flat_dshift(mp, COV, mp.get_bvect_cart()) ;
flat_dshift.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
flat_dshift.set(i,j) = flat.cov()(j,1)*(shift(1).deriv(i))
+ flat.cov()(j,2)*(shift(2).deriv(i))
+ flat.cov()(j,3)*(shift(3).deriv(i))
+ flat.cov()(i,1)*(shift(1).deriv(j))
+ flat.cov()(i,2)*(shift(2).deriv(j))
+ flat.cov()(i,3)*(shift(3).deriv(j))
- 2. * divshift * flat.cov()(i,j) / 3. ;
}
}
Sym_tensor taij_down(mp, COV, mp.get_bvect_cart()) ;
taij_down.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
taij_down.set(i,j) = 0.5 * pow(confo, 7.) * flat_dshift(i,j)
/ lapconf ;
}
}
taij_down.std_spectral_base() ;
taij_down.annule_domain(0) ;
taij_quad = 0. ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
taij_quad += taij_down(i,j) * taij(i,j) ;
}
}
taij_quad.std_spectral_base() ;
}
}
}
Lorene::Black_hole::shift_rs
Vector shift_rs
Part of the shift vector from the numerical computation.
Definition blackhole.h:112
Lorene::Black_hole::taij_quad
Scalar taij_quad
Part of the scalar generated by .
Definition blackhole.h:135
Lorene::Black_hole::taij_quad_rs
Scalar taij_quad_rs
Part of the scalar.
Definition blackhole.h:138
Lorene::Black_hole::taij
Sym_tensor taij
Trace of the extrinsic curvature.
Definition blackhole.h:127
Lorene::Black_hole::taij_rs
Sym_tensor taij_rs
Part of the extrinsic curvature tensor.
Definition blackhole.h:130
Lorene::Black_hole::mp
Map & mp
Mapping associated with the black hole.
Definition blackhole.h:80
Lorene::Black_hole::lapconf
Scalar lapconf
A function (lapse function * conformal factor) lapconf generated by the black hole.
Definition blackhole.h:97
Lorene::Black_hole::shift
Vector shift
Shift vector generated by the black hole.
Definition blackhole.h:109
Lorene::Black_hole::flat
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the black hol...
Definition blackhole.h:143
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::Black_hole::mass_bh
double mass_bh
Gravitational mass of BH.
Definition blackhole.h:88
Lorene::Black_hole::extr_curv_bh
void extr_curv_bh()
Computes taij , taij_quad from shift , lapse , confo .
Definition blackhole_extr_curv.C:62
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::inc_dzpuis
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Definition scalar_r_manip.C:470
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
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::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::Tensor::inc_dzpuis
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:816
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
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::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
Lorene::Tensor::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition tensor.C:926
Lorene
Lorene prototypes.
Definition app_hor.h:64
Unites
Standard units of space, time and mass.