excised_slice.C
/*
* Methods of class Excised_slice
*
* (see file excised_slice.h for documentation)
* Disclaimer: the class Isol_hole() is redundant with this class under a set of parameters;
* therefore it has to go at some point.
*/
/*
* Copyright (c) 2010 Nicolas Vasset
* 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 as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* 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
*
*/
// Headers C
#include "math.h"
// Headers Lorene
#include "excised_slice.h"
#include "spheroid.h"
#include "excision_surf.h"
#include "excision_hor.h"
#include "utilitaires.h"
#include "param.h"
#include "unites.h"
#include "proto.h"
namespace Lorene {
// Fundamental constants and units
// -------------------------------
using namespace Unites ;
//--------------//
// Constructors //
//--------------//
// Standard constructor
// --------------------
Excised_slice::Excised_slice (const Map& mpi, int hor_type, int scheme_type)
: mp(mpi),
type_hor(hor_type),
field_set(scheme_type),
lapse(mpi),
conf_fact(mpi),
shift(mpi, CON, mpi.get_bvect_spher()),
hij(mpi, CON, mpi.get_bvect_spher()),
hatA(mpi, CON, mpi.get_bvect_spher()),
Xx(mpi, CON, mpi.get_bvect_spher()){
// Pointers of derived quantities initialized to zero :
set_der_0x0() ;
// Initializing primary quantities.
lapse = 1. ;
conf_fact = 1. ;
shift.set_etat_zero() ;
hij.set_etat_zero();
hatA.set_etat_zero();
Xx.set_etat_zero();
}
// Copy constructor
// ----------------
Excised_slice::Excised_slice(const Excised_slice& ih)
: mp(ih.mp),
type_hor(ih.type_hor),
field_set(ih.field_set),
lapse(ih.lapse),
conf_fact(ih.conf_fact),
shift(ih.shift),
hij(ih.hij),
hatA(ih.hatA),
Xx(ih.Xx){
set_der_0x0() ;
}
// Constructor from a file
// -----------------------
Excised_slice::Excised_slice(const Map& mpi, int hor_type, int scheme_type, FILE* fich)
: mp(mpi),
type_hor(hor_type),
field_set(scheme_type),
lapse(mpi,*(mpi.get_mg()), fich),
conf_fact(mpi,*(mpi.get_mg()), fich),
shift(mpi, mpi.get_bvect_spher(), fich),
hij(mpi, mpi.get_bvect_spher(), fich),
hatA(mpi, mpi.get_bvect_spher(), fich),
Xx(mpi, mpi.get_bvect_spher(), fich){
// Pointers of derived quantities initialized to zero
// --------------------------------------------------
set_der_0x0() ;
}
//------------//
// Destructor //
//------------//
Excised_slice::~Excised_slice(){
del_deriv() ;
}
//----------------------------------//
// Management of derived quantities //
//----------------------------------//
void Excised_slice::del_deriv() const {
if (p_hor != 0x0) delete p_hor ;
if (p_adm_mass != 0x0) delete p_adm_mass ;
if (p_komar_angmom != 0x0) delete p_komar_angmom ;
if (p_virial_residue != 0x0) delete p_virial_residue ;
Excised_slice::set_der_0x0() ;
}
void Excised_slice::set_der_0x0() const {
p_hor = 0x0;
p_adm_mass = 0x0 ;
p_komar_angmom = 0x0 ;
p_virial_residue = 0x0 ;
}
//--------------//
// Assignment //
//--------------//
// Assignment to another Excised_slice
// ----------------------------
void Excised_slice::operator=(const Excised_slice& ih) {
assert( &(ih.mp) == &mp ) ; // Same mapping
type_hor = ih.type_hor ;
field_set = ih.field_set ;
lapse = ih.lapse ;
conf_fact = ih.conf_fact ;
shift = ih.shift ;
hij = ih.hij ;
hatA = ih.hatA ;
Xx = ih.Xx;
del_deriv() ; // Deletes all derived quantities
}
//--------------//
// Outputs //
//--------------//
// Save in a file
// --------------
void Excised_slice::sauve(FILE* fich) const {
lapse.sauve(fich) ;
conf_fact.sauve(fich) ;
shift.sauve(fich);
hij.sauve(fich);
hatA.sauve(fich);
Xx.sauve(fich);}
// Prints out maximal errors in Einstein equations for the obtained metric fields
void Excised_slice::Einstein_errors() {
const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
const Metric_flat& mets = (*map).flat_met_spher() ;
Sym_tensor gamtcon = mets.con() + hij;
Metric gamt(gamtcon);
Sym_tensor gamcon = gamtcon/(conf_fact*conf_fact*conf_fact*conf_fact);
gamcon.std_spectral_base();
Metric gam(gamcon);
Sym_tensor k_uu = hatA/(pow(conf_fact,10)) ;
k_uu.std_spectral_base();
k_uu.dec_dzpuis(k_uu(1,1).get_dzpuis()); //WTF?
Sym_tensor k_dd = k_uu.up_down(gam);
Scalar TrK3 = k_uu.trace(gam);
// Hamiltonian constraint
//-----------------------
Scalar ham_constr = gam.ricci_scal() ;
ham_constr.dec_dzpuis(3) ; // To check
ham_constr += TrK3*TrK3 - contract(k_uu, 0, 1, k_dd, 0, 1) ;
maxabs(ham_constr, "Hamiltonian constraint: ") ;
// Momentum constraint
//-------------------
Vector mom_constr = k_uu.divergence(gam) - TrK3.derive_con(gam) ;
mom_constr.dec_dzpuis(2) ; // To check
maxabs(mom_constr, "Momentum constraint: ") ;
// Evolution equations
//--------------------
Sym_tensor evol_eq = lapse*gam.ricci()
- lapse.derive_cov(gam).derive_cov(gam);
evol_eq.dec_dzpuis() ;
evol_eq += k_dd.derive_lie(shift) ;
evol_eq.dec_dzpuis(2) ; // To check
evol_eq += lapse*(TrK3*k_dd - 2*contract(k_dd, 1, k_dd.up(0, gam), 0) ) ;
maxabs(evol_eq, "Evolution equations: ") ;
return;
}
//----------------------------//
// Accessors/ Derived data //
//----------------------------//
// Computation of the Spheroid corresponding to the black hole excision surface.
Spheroid Excised_slice::hor() {
const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
const Mg3d* mgrid = (*map).get_mg();
// Construct angular grid for h(theta,phi)
const Mg3d* g_angu = (*mgrid).get_angu_1dom() ;
const Coord& rr = (*map).r;
Scalar rrr (*map) ;
rrr = rr ;
rrr.std_spectral_base();
assert((rrr.val_grid_point(1,0,0,0) - 1.) <= 1.e-9); // For now the code handles only horizons at r=1, corresponding to the first shell inner boundary. This test assures this is the case with our mapping.
// Angular mapping defined for one domain (argument of spheroid Class)
//--------------------------------------------------------------------
double r_limits2[] = {rrr.val_grid_point(1,0,0,0), rrr.val_grid_point(2,0,0,0)} ;
const Map_af map_2(*g_angu, r_limits2);
//Full 3-metric and extrrinsic curvature
const Metric_flat& mets = (*map).flat_met_spher() ;
Sym_tensor gamtcon = mets.con() + hij;
Metric gamt(gamtcon);
Sym_tensor gamcon = gamtcon/(conf_fact*conf_fact*conf_fact*conf_fact);
Metric gam(gamcon);
Sym_tensor kuu = hatA/pow(conf_fact,10) ;
kuu.std_spectral_base();
Sym_tensor kdd = kuu.up_down(gam);
//---------------------------------------------------------
// Construction of the spheroid associated with those data
//--------------------------------------------------------
double hor_posd = rrr.val_grid_point(1,0,0,0);
Scalar hor_pos(map_2); hor_pos = hor_posd; hor_pos.std_spectral_base();
Spheroid hor_loc(hor_pos, gam, kdd);
return hor_loc;
}
// Computation of the ADM mass of the BH spacetime
double Excised_slice::adm_mass() {
const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
const Metric_flat& mets = (*map).flat_met_spher() ;
Sym_tensor gamtcon = mets.con() + hij;
Metric gamt(gamtcon);
Sym_tensor gamcon = gamtcon/(conf_fact*conf_fact*conf_fact*conf_fact);
Metric gam(gamcon);
Scalar detgam = sqrt((gam.cov())(2,2)*(gam.cov())(3,3) - (gam.cov())(2,3)*(gam.cov())(3,2));
detgam.std_spectral_base();
Vector corr_adm = - (0.125*contract(gamt.cov().derive_con(mets),1,2));
Scalar admintegr = conf_fact.dsdr() + corr_adm(1);
double M_ADM = - (1/(2.*3.1415927*ggrav))*(*map).integrale_surface_infini(admintegr*detgam);
return M_ADM;
}
// Computation of the Komar angular momentum w.r.t. assumed rotational symmetry
double Excised_slice:: komar_angmom() {
const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
const Metric_flat& mets = (*map).flat_met_spher() ;
Sym_tensor gamtcon = mets.con() + hij;
Sym_tensor gamcon = gamtcon/(conf_fact*conf_fact*conf_fact*conf_fact);
gamcon.std_spectral_base();
Metric gam(gamcon);
Sym_tensor k_uu = hatA/(pow(conf_fact,10));
k_uu.std_spectral_base();
k_uu.dec_dzpuis(k_uu(1,1).get_dzpuis()); //WTF?
Sym_tensor k_dd = k_uu.up_down(gam);
Scalar detgam = sqrt((gam.cov())(2,2)*(gam.cov())(3,3) - (gam.cov())(2,3)*(gam.cov())(3,2));
detgam.std_spectral_base();
Scalar contraction = k_dd(1,3); contraction.mult_r_dzpuis(2); contraction.mult_sint();
double angu_komar = - (1/(8.*3.1415927*ggrav))*(*map).integrale_surface_infini(detgam*contraction);
return angu_komar;
}
// Computation of the Virial residual, as rescaled difference at infinity
// between the ADM mass and the Komar integral associated to the mass
double Excised_slice::virial_residue() {
const Mg3d* mgrid = mp.get_mg();
const int nz = (*mgrid).get_nzone(); // Number of domains
Valeur** devel_psi (conf_fact.asymptot(1)) ;
Valeur** devel_n (lapse.asymptot(1)) ;
double erreur = (2*(*devel_psi[1])(nz-1, 0, 0, 0)
+ (*devel_n[1])(nz-1, 0, 0, 0))/fabs ((*devel_n[1])(nz-1, 0, 0, 0)) ;
return erreur;
}
}
Lorene::Coord
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Lorene::Excised_slice::hatA
Sym_tensor hatA
Rescaled tracefree extrinsic curvature tensor: rescaled same way as Cordero et al.
Definition excised_slice.h:74
Lorene::Excised_slice::p_adm_mass
double * p_adm_mass
Computation of the ADM mass of the BH spacetime.
Definition excised_slice.h:91
Lorene::Excised_slice::komar_angmom
double komar_angmom()
Computation of the Komar angular momentum w.r.t.
Definition excised_slice.C:308
Lorene::Excised_slice::type_hor
int type_hor
Chosen horizon type.
Definition excised_slice.h:52
Lorene::Excised_slice::Einstein_errors
void Einstein_errors()
Prints out errors in Einstein equations for the data obtained.
Definition excised_slice.C:193
Lorene::Excised_slice::p_komar_angmom
double * p_komar_angmom
Computation of the Komar angular momentum w.r.t.
Definition excised_slice.h:96
Lorene::Excised_slice::del_deriv
virtual void del_deriv() const
Deletes all the derived quantities.
Definition excised_slice.C:135
Lorene::Excised_slice::p_hor
Spheroid * p_hor
Computation of the spheroid associated with the black hole horizon.
Definition excised_slice.h:89
Lorene::Excised_slice::virial_residue
double virial_residue()
Computation of the Virial residual, as difference at infinity between the ADM mass and the Komar inte...
Definition excised_slice.C:333
Lorene::Excised_slice::hij
Sym_tensor hij
Deviation tensor( non-flat part of the conformal 3-metric on the slice; see Bonazzola et al.
Definition excised_slice.h:69
Lorene::Excised_slice::set_der_0x0
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition excised_slice.C:144
Lorene::Excised_slice::lapse
Scalar lapse
Lapse function.
Definition excised_slice.h:58
Lorene::Excised_slice::p_virial_residue
double * p_virial_residue
Computation of the Virial residual, as difference at infinity between the ADM mass and the Komar inte...
Definition excised_slice.h:101
Lorene::Excised_slice::operator=
void operator=(const Excised_slice &)
Assignment to another Excised_slice.
Definition excised_slice.C:159
Lorene::Excised_slice::~Excised_slice
virtual ~Excised_slice()
Destructor.
Definition excised_slice.C:124
Lorene::Excised_slice::mp
const Map & mp
Mapping associated with the slice.
Definition excised_slice.h:46
Lorene::Excised_slice::field_set
int field_set
Chose field set type.
Definition excised_slice.h:55
Lorene::Excised_slice::shift
Vector shift
Shift vector.
Definition excised_slice.h:64
Lorene::Excised_slice::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition excised_slice.C:181
Lorene::Excised_slice::Excised_slice
Excised_slice(const Map &mp_i, int hor_type, int scheme_type)
Standard constructor.
Definition excised_slice.C:58
Lorene::Excised_slice::adm_mass
double adm_mass()
Computation of the ADM mass of the BH spacetime.
Definition excised_slice.C:288
Lorene::Excised_slice::Xx
Vector Xx
Longitudinal part of the extrinsic curvature.
Definition excised_slice.h:79
Lorene::Excised_slice::hor
Spheroid hor()
Spheroid at the excised boundary associated with the black hole MOTS on the slice.
Definition excised_slice.C:246
Lorene::Map_af
Affine radial mapping.
Definition map.h:2027
Lorene::Metric_flat
Flat metric for tensor calculation.
Definition metric.h:261
Lorene::Metric_flat::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric_flat.C:153
Lorene::Metric
Metric for tensor calculation.
Definition metric.h:90
Lorene::Metric::cov
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
Lorene::Metric::ricci
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect ).
Definition metric.C:338
Lorene::Metric::ricci_scal
virtual const Scalar & ricci_scal() const
Returns the Ricci scalar.
Definition metric.C:350
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_angu_1dom
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
Definition mg3d.C:494
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::dsdr
const Scalar & dsdr() const
Returns of *this .
Definition scalar_deriv.C:113
Lorene::Scalar::mult_sint
void mult_sint()
Multiplication by .
Definition scalar_th_manip.C:84
Lorene::Scalar::mult_r_dzpuis
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r .
Definition scalar_r_manip.C:221
Lorene::Scalar::derive_con
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition scalar_deriv.C:402
Lorene::Scalar::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Definition scalar_r_manip.C:418
Lorene::Spheroid
Spheroidal 2-surfaces embedded in a time-slice of the 3+1 formalism.
Definition spheroid.h:84
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Lorene::Sym_tensor::derive_lie
Sym_tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition sym_tensor.C:360
Lorene::Sym_tensor::divergence
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition sym_tensor.C:349
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
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::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:808
Lorene::Tensor::up_down
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
Definition tensor_calculus.C:305
Lorene::Tensor_sym::derive_con
const Tensor_sym & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition tensor_sym_calculus.C:204
Lorene::Tensor::up
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
Definition tensor_calculus.C:225
Lorene::Tensor::trace
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Definition tensor_calculus.C:94
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::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene::maxabs
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Definition tensor_calculus_ext.C:606
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::get_bvect_spher
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
Lorene::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:139