compobj.C
/*
* Methods of the class Compobj
*
* (see file compobj.h for documentation).
*
*/
/*
* Copyright (c) 2012 Claire Some, Eric Gourgoulhon
*
* 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 compobj_C[] = "$Header: /cvsroot/Lorene/C++/Source/Compobj/compobj.C,v 1.9 2014/10/13 08:52:49 j_novak Exp $" ;
/*
* $Id: compobj.C,v 1.9 2014/10/13 08:52:49 j_novak Exp $
* $Log: compobj.C,v $
* Revision 1.9 2014/10/13 08:52:49 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.8 2014/05/16 11:55:18 o_straub
* fixed: GYOTO output from compobj & compobj_QI
*
* Revision 1.7 2014/04/11 17:22:07 o_straub
* Risco and Rms output for GYOTO
*
* Revision 1.6 2013/07/25 19:44:11 o_straub
* calculation of the marginally bound radius
*
* Revision 1.5 2013/04/03 12:10:13 e_gourgoulhon
* Added member kk to Compobj; suppressed tkij
*
* Revision 1.4 2012/12/03 15:27:30 c_some
* Small changes
*
* Revision 1.3 2012/11/22 16:04:51 c_some
* Minor modifications
*
* Revision 1.2 2012/11/20 16:24:09 c_some
* Added computation of ADM mass (method mass_q())
*
* Revision 1.1 2012/11/15 16:20:51 c_some
* New class Compobj
*
*
* $Header: /cvsroot/Lorene/C++/Source/Compobj/compobj.C,v 1.9 2014/10/13 08:52:49 j_novak Exp $
*
*/
// C headers
#include <cassert>
#include <cmath>
// Lorene headers
#include "compobj.h"
#include "nbr_spx.h"
#include "utilitaires.h"
//--------------//
// Constructors //
//--------------//
// Standard constructor
// --------------------
namespace Lorene {
Compobj::Compobj(Map& map_i) :
mp(map_i) ,
nn(map_i) ,
beta(map_i, CON, map_i.get_bvect_spher()) ,
gamma(map_i.flat_met_spher()) ,
ener_euler(map_i) ,
mom_euler(map_i, CON, map_i.get_bvect_spher()) ,
stress_euler(map_i, COV, map_i.get_bvect_spher()) ,
kk(map_i, COV, map_i.get_bvect_spher())
{
// Pointers of derived quantities initialized to zero :
set_der_0x0() ;
// Some initialisations:
nn = 1 ;
nn.std_spectral_base() ;
beta.set_etat_zero() ;
ener_euler = 0 ;
mom_euler.set_etat_zero() ;
stress_euler.set_etat_zero() ;
kk.set_etat_zero() ;
}
// Copy constructor
// --------------------
Compobj::Compobj(const Compobj& co) :
mp(co.mp) ,
nn(co.nn) ,
beta(co.beta) ,
gamma(co.gamma) ,
ener_euler(co.ener_euler) ,
mom_euler(co.mom_euler) ,
stress_euler(co.stress_euler) ,
kk(co.kk)
{
// Pointers of derived quantities initialized to zero :
set_der_0x0() ;
}
// Constructor from a file
// -----------------------
Compobj::Compobj(Map& map_i, FILE* fich) :
mp(map_i) ,
nn(map_i, *(map_i.get_mg()), fich) ,
beta(map_i, map_i.get_bvect_spher(), fich) ,
gamma(map_i, fich) ,
ener_euler(map_i, *(map_i.get_mg()), fich) ,
mom_euler(map_i, map_i.get_bvect_spher(), fich) ,
stress_euler(map_i, map_i.get_bvect_spher(), fich) ,
kk(map_i, COV, map_i.get_bvect_spher())
{
// Pointers of derived quantities initialized to zero :
set_der_0x0() ;
}
//------------//
// Destructor //
//------------//
Compobj::~Compobj(){
del_deriv() ;
}
//----------------------------------//
// Management of derived quantities //
//----------------------------------//
void Compobj::del_deriv() const {
if (p_adm_mass != 0x0) delete p_adm_mass ;
Compobj::set_der_0x0() ;
}
void Compobj::set_der_0x0() const {
p_adm_mass = 0x0 ;
}
//--------------//
// Assignment //
//--------------//
// Assignment to another Compobj
// ----------------------------
void Compobj::operator=(const Compobj& co) {
assert( &(co.mp) == &mp ) ; // Same mapping
nn = co.nn ;
beta = co.beta ;
gamma = co.gamma ;
ener_euler = co.ener_euler ;
mom_euler = co.mom_euler ;
stress_euler = co.stress_euler ;
kk = co.kk ;
del_deriv() ; // Deletes all derived quantities
}
//--------------//
// Outputs //
//--------------//
// Save in a file
// --------------
void Compobj::sauve(FILE* fich) const {
nn.sauve(fich) ;
beta.sauve(fich) ;
gamma.sauve(fich) ;
ener_euler.sauve(fich) ;
mom_euler.sauve(fich) ;
stress_euler.sauve(fich) ;
}
// Save in a file for GYOTO input
// ------------------------------
void Compobj::gyoto_data(const char* file_name) const {
FILE* file_out = fopen(file_name, "w") ;
double total_time = 0. ; // for compatibility
fwrite_be(&total_time, sizeof(double), 1, file_out) ;
mp.get_mg()->sauve(file_out) ;
mp.sauve(file_out) ;
nn.sauve(file_out) ;
beta.sauve(file_out) ;
gamma.cov().sauve(file_out) ;
gamma.con().sauve(file_out) ;
kk.sauve(file_out) ;
fclose(file_out) ;
cout << "WRITING TO GYOTO FILE - OK: " << endl ;
}
// Printing
// --------
ostream& operator<<(ostream& ost, const Compobj& co) {
co >> ost ;
return ost ;
}
ostream& Compobj::operator>>(ostream& ost) const {
ost << endl << "Compact object (class Compobj) " << endl ;
ost << "Mapping : " << mp << endl ;
ost << "Central values of various fields : " << endl ;
ost << "-------------------------------- " << endl ;
ost << " lapse function : N_c = " << nn.val_grid_point(0,0,0,0) << endl ;
ost << " metric components gamma_{ij} : " << endl
<< " ( " << gamma.cov()(1,1).val_grid_point(0,0,0,0) << " "
<< gamma.cov()(1,2).val_grid_point(0,0,0,0) << " "
<< gamma.cov()(1,3).val_grid_point(0,0,0,0) << " )" << endl
<< " ( " << gamma.cov()(2,1).val_grid_point(0,0,0,0) << " "
<< gamma.cov()(2,2).val_grid_point(0,0,0,0) << " "
<< gamma.cov()(2,3).val_grid_point(0,0,0,0) << " )" << endl
<< " ( " << gamma.cov()(3,1).val_grid_point(0,0,0,0) << " "
<< gamma.cov()(3,2).val_grid_point(0,0,0,0) << " "
<< gamma.cov()(3,3).val_grid_point(0,0,0,0) << " )" << endl ;
ost << " components of the extrinsic curvature K_{ij} : " << endl
<< " ( " << kk(1,1).val_grid_point(0,0,0,0) << " "
<< kk(1,2).val_grid_point(0,0,0,0) << " "
<< kk(1,3).val_grid_point(0,0,0,0) << " )" << endl
<< " ( " << kk(2,1).val_grid_point(0,0,0,0) << " "
<< kk(2,2).val_grid_point(0,0,0,0) << " "
<< kk(2,3).val_grid_point(0,0,0,0) << " )" << endl
<< " ( " << kk(3,1).val_grid_point(0,0,0,0) << " "
<< kk(3,2).val_grid_point(0,0,0,0) << " "
<< kk(3,3).val_grid_point(0,0,0,0) << " )" << endl ;
ost << " energy density / Eulerian observer : E_c = " << ener_euler.val_grid_point(0,0,0,0) << endl ;
ost << " components of the stress tensor S_{ij} / Eulerian observer : " << endl
<< " ( " << stress_euler(1,1).val_grid_point(0,0,0,0) << " "
<< stress_euler(1,2).val_grid_point(0,0,0,0) << " "
<< stress_euler(1,3).val_grid_point(0,0,0,0) << " )" << endl
<< " ( " << stress_euler(2,1).val_grid_point(0,0,0,0) << " "
<< stress_euler(2,2).val_grid_point(0,0,0,0) << " "
<< stress_euler(2,3).val_grid_point(0,0,0,0) << " )" << endl
<< " ( " << stress_euler(3,1).val_grid_point(0,0,0,0) << " "
<< stress_euler(3,2).val_grid_point(0,0,0,0) << " "
<< stress_euler(3,3).val_grid_point(0,0,0,0) << " )" << endl ;
//## ost << endl << "ADM mass : " << adm_mass() << endl ;
return ost ;
}
//-------------------------//
// Computational methods //
//-------------------------//
// Extrinsic curvature
void Compobj::extrinsic_curvature() {
cout << "WARNING: Compobj::extrinsic_curvature() NOT TESTED !" << endl ;
// Gradient of the shift D_j beta_i
Vector cobeta = beta.down(0, gamma) ;
Tensor dn = cobeta.derive_cov(gamma) ;
kk.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
for (int j=i; j<=3; j++) {
kk.set(i, j) = (dn(i, j) + dn(j, i))/(2*nn) ;
}
}
}
// Gravitational mass
double Compobj::adm_mass() const {
if (p_adm_mass == 0x0) { // a new computation is required
const Sym_tensor& gam_dd = gamma.cov() ; // components \gamma_{ij} of the 3-metric
Metric_flat ff(mp, *(gam_dd.get_triad())) ;
Vector ww = gam_dd.derive_con(ff).trace(1,2).up(0,ff)
- gam_dd.trace(ff).derive_con(ff) ;
p_adm_mass = new double( ww.flux(__infinity, ff) / (16.* M_PI) ) ;
}
return *p_adm_mass ;
}
}
Lorene::Compobj::adm_mass
virtual double adm_mass() const
ADM mass (computed as a surface integral at spatial infinity).
Definition compobj.C:310
Lorene::Compobj::kk
Sym_tensor kk
Extrinsic curvature tensor .
Definition compobj.h:153
Lorene::Compobj::mom_euler
Vector mom_euler
Total 3-momentum density in the Eulerian frame.
Definition compobj.h:147
Lorene::Compobj::operator<<
friend ostream & operator<<(ostream &, const Compobj &)
Display.
Definition compobj.C:233
Lorene::Compobj::del_deriv
virtual void del_deriv() const
Deletes all the derived quantities.
Definition compobj.C:155
Lorene::Compobj::stress_euler
Sym_tensor stress_euler
Stress tensor with respect to the Eulerian observer.
Definition compobj.h:150
Lorene::Compobj::ener_euler
Scalar ener_euler
Total energy density E in the Eulerian frame.
Definition compobj.h:144
Lorene::Compobj::set_der_0x0
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition compobj.C:163
Lorene::Compobj::operator=
void operator=(const Compobj &)
Assignment to another Compobj.
Definition compobj.C:175
Lorene::Compobj::operator>>
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition compobj.C:239
Lorene::Compobj::gamma
Metric gamma
3-metric
Definition compobj.h:141
Lorene::Compobj::Compobj
Compobj(Map &map_i)
Standard constructor.
Definition compobj.C:82
Lorene::Compobj::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition compobj.C:196
Lorene::Compobj::nn
Scalar nn
Lapse function N .
Definition compobj.h:135
Lorene::Compobj::beta
Vector beta
Shift vector .
Definition compobj.h:138
Lorene::Compobj::gyoto_data
void gyoto_data(const char *file_name) const
Save in a file for GYOTO.
Definition compobj.C:210
Lorene::Compobj::p_adm_mass
double * p_adm_mass
ADM mass.
Definition compobj.h:158
Lorene::Compobj::extrinsic_curvature
virtual void extrinsic_curvature()
Computation of the extrinsic curvature.
Definition compobj.C:290
Lorene::Compobj::~Compobj
virtual ~Compobj()
Destructor.
Definition compobj.C:144
Lorene::Compobj::mp
Map & mp
Mapping describing the coordinate system (r,theta,phi).
Definition compobj.h:132
Lorene::Metric_flat
Flat metric for tensor calculation.
Definition metric.h:261
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Lorene::Tensor
Tensor handling.
Definition tensor.h:288
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::flux
double flux(double radius, const Metric &met) const
Computes the flux of the vector accross a sphere r = const.
Definition vector.C:807
Lorene::fwrite_be
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition fwrite_be.C:70
Lorene::Tensor::derive_con
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition tensor.C:1014
Lorene::Tensor::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:866
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::derive_cov
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor.C:1002
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
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
Lorene::flat_met_spher
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition map.C:321