tslice_dirac_max_solve.C
/*
* Methods of class Tslice_dirac_max for solving Einstein equations
*
* (see file time_slice.h for documentation).
*
*/
/*
* Copyright (c) 2004 Eric Gourgoulhon & Jerome Novak
*
* 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 tslice_dirac_max_solve_C[] = "$Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_dirac_max_solve.C,v 1.18 2014/10/13 08:53:48 j_novak Exp $" ;
/*
* $Id: tslice_dirac_max_solve.C,v 1.18 2014/10/13 08:53:48 j_novak Exp $
* $Log: tslice_dirac_max_solve.C,v $
* Revision 1.18 2014/10/13 08:53:48 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.17 2014/10/06 15:13:22 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.16 2010/10/20 07:58:10 j_novak
* Better implementation of the explicit time-integration. Not fully-tested yet.
*
* Revision 1.15 2008/12/02 15:02:22 j_novak
* Implementation of the new constrained formalism, following Cordero et al. 2009
* paper. The evolution eqs. are solved as a first-order system. Not tested yet!
*
* Revision 1.14 2004/07/08 12:29:01 j_novak
* use of new method Tensor::annule_extern_cn
*
* Revision 1.13 2004/06/30 08:02:40 j_novak
* Added filtering in l of khi_new and mu_new. ki_source is forced to go to
* zero at least as r^2.
*
* Revision 1.12 2004/06/17 07:07:11 e_gourgoulhon
* Method solve_hij:
* -- replaced the attenuation of khi_source with tempo by a call to
* the new method Tensor::annule_extern_c2
* -- suppressed filtre_r on khi_source and khi_new
* -- added graphs of W^i and LW^{ij} (transverse decomp of S^ij).
*
* Revision 1.11 2004/06/15 09:43:36 j_novak
* Attenuation of the source for khi in the last shell (temporary?).
*
* Revision 1.10 2004/06/14 20:47:31 e_gourgoulhon
* Added argument method_poisson to method solve_hij.
*
* Revision 1.9 2004/06/03 10:02:45 j_novak
* Some filtering is done on source_khi and khi_new.
*
* Revision 1.8 2004/05/24 21:00:44 e_gourgoulhon
* Method solve_hij: khi and mu.smooth_decay(2,2) --> smooth_decay(2,1) ;
* added exponential_decay(khi) and exponential_decay(mu) after the
* call to smooth_decay. Method exponential_decay is provisory defined
* in this file.
*
* Revision 1.7 2004/05/17 19:56:25 e_gourgoulhon
* -- Method solve_beta: added argument method
* -- Method solve_hij: added argument graph_device
*
* Revision 1.6 2004/05/12 15:24:20 e_gourgoulhon
* Reorganized the #include 's, taking into account that
* time_slice.h contains now an #include "metric.h".
*
* Revision 1.5 2004/05/05 14:47:05 e_gourgoulhon
* Modified text and graphical outputs.
*
* Revision 1.4 2004/05/03 15:06:27 e_gourgoulhon
* Added matter source in solve_hij.
*
* Revision 1.3 2004/05/03 14:50:00 e_gourgoulhon
* Finished the implementation of method solve_hij.
*
* Revision 1.2 2004/04/30 14:36:15 j_novak
* Added the method Tslice_dirac_max::solve_hij(...)
* NOT READY YET!!!
*
* Revision 1.1 2004/04/30 10:52:14 e_gourgoulhon
* First version.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_dirac_max_solve.C,v 1.18 2014/10/13 08:53:48 j_novak Exp $
*
*/
// C headers
#include <cassert>
// Lorene headers
#include "time_slice.h"
#include "unites.h"
#include "graphique.h"
#include "proto.h"
//----------------------------//
// Equation for N\Psi //
//----------------------------//
namespace Lorene {
Scalar Tslice_dirac_max::solve_npsi(const Scalar* p_ener_dens,
const Scalar* p_trace_stress) const {
using namespace Unites ;
const Map& map = npsi().get_mp() ;
Scalar ener_dens(map) ;
if (p_ener_dens != 0x0) ener_dens = *(p_ener_dens) ;
else ener_dens.set_etat_zero() ;
Scalar trace_stress(map) ;
if (p_trace_stress != 0x0) trace_stress = *(p_trace_stress) ;
else trace_stress.set_etat_zero() ;
// Source for N\Psi
// ----------------
const Vector& dnpsi = npsi().derive_cov(ff) ; // D_i N\Psi
const Tensor_sym& dhh = hh().derive_cov(ff) ; // D_k h^{ij}
const Tensor_sym& dtgam = tgam().cov().derive_cov(ff) ;
// D_k {\tilde \gamma}_{ij}
Sym_tensor taa = aa().up_down(tgam()) ;
Scalar aa_quad = contract(taa, 0, 1, aa(), 0, 1) ;
Scalar tildeR = 0.25 * contract( dhh, 0, 1, dtgam, 0, 1 ).trace(tgam())
- 0.5 * contract( dhh, 0, 1, dtgam, 0, 2 ).trace(tgam()) ;
Scalar source_npsi = - contract( hh(), 0, 1, dnpsi.derive_cov(ff), 0, 1 ) ;
source_npsi.inc_dzpuis() ;
source_npsi += npsi()*( 0.5*qpig*(ener_dens + 2.*trace_stress )/( psi()*psi() )
+ 0.875*aa_quad/( psi4()*psi4() ) + 0.125*tildeR ) ;
// Resolution of the Poisson equation for N\Psi
// --------------------------------------------
Scalar npsi_new = source_npsi.poisson() + 1. ;
if (npsi_new.get_etat() == ETATUN) npsi_new.std_spectral_base() ;
#ifndef NDEBUG
// Test:
maxabs(npsi_new.laplacian() - source_npsi,
"Absolute error in the resolution of the equation for N") ;
#endif
return npsi_new ;
}
//--------------------------//
// Equation for \Psi //
//--------------------------//
Scalar Tslice_dirac_max::solve_psi(const Scalar* p_ener_dens) const {
using namespace Unites ;
const Map& map = psi().get_mp() ;
Scalar ener_dens(map) ;
if (p_ener_dens != 0x0) ener_dens = *(p_ener_dens) ;
else ener_dens.set_etat_zero() ;
// Source for \Psi
// ---------------
const Vector& dpsi = psi().derive_cov(ff) ; // D_i Psi
const Tensor_sym& dhh = hh().derive_cov(ff) ; // D_k h^{ij}
const Tensor_sym& dtgam = tgam().cov().derive_cov(ff) ;
// D_k {\tilde \gamma}_{ij}
Sym_tensor taa = hata().up_down(tgam()) ;
Scalar aa_quad = contract(taa, 0, 1, hata(), 0, 1) ;
Scalar tildeR = 0.25 * contract( dhh, 0, 1, dtgam, 0, 1 ).trace(tgam())
- 0.5 * contract( dhh, 0, 1, dtgam, 0, 2 ).trace(tgam()) ;
Scalar source_psi = -contract( hh(), 0, 1, dpsi.derive_cov(ff), 0, 1 ) ;
source_psi.inc_dzpuis() ;
source_psi -= 0.5*qpig*ener_dens/psi()
+ 0.125*( aa_quad*pow(psi(), -7) - tildeR*psi() ) ;
// Resolution of the Poisson equation for Psi
// ------------------------------------------
Scalar psi_new = source_psi.poisson() + 1. ;
if (psi_new.get_etat() == ETATUN) psi_new.std_spectral_base() ;
#ifndef NDEBUG
// Test:
maxabs(psi_new.laplacian() - source_psi,
"Absolute error in the resolution of the equation for Psi") ;
#endif
return psi_new ;
}
//--------------------------//
// Equation for beta //
//--------------------------//
Vector Tslice_dirac_max::solve_beta(int method)
const {
// Source for beta
// ---------------
Vector source_beta =
- contract(hh(), 0, 1, beta().derive_cov(ff).derive_cov(ff), 1, 2)
- 0.3333333333333333*contract(hh(), 1, beta().divergence(ff).derive_cov(ff), 0) ;
Sym_tensor sym_tmp = 2*nn()*aa() ;
source_beta += sym_tmp.divergence(ff) ;
source_beta.inc_dzpuis() ; // dzpuis: 3 -> 4
// Resolution of the vector Poisson equation
//------------------------------------------
Vector beta_new = source_beta.poisson(0.3333333333333333, ff, method) ;
#ifndef NDEBUG
// Test:
Vector test_beta = (beta_new.derive_con(ff)).divergence(ff)
+ 0.3333333333333333 * (beta_new.divergence(ff)).derive_con(ff) ;
test_beta.inc_dzpuis() ;
maxabs(test_beta - source_beta,
"Absolute error in the resolution of the equation for beta") ;
#endif
return beta_new ;
}
}
Lorene::Metric::cov
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::derive_cov
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this.
Definition scalar_deriv.C:390
Lorene::Scalar::poisson
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition scalar_pde.C:136
Lorene::Scalar::laplacian
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition scalar_deriv.C:436
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::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
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::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
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::Tensor_sym
Symmetric tensors (with respect to two of their arguments).
Definition tensor.h:1037
Lorene::Time_slice_conf::npsi
virtual const Scalar & npsi() const
Factor at the current time step (jtime ).
Definition time_slice_conf.C:732
Lorene::Time_slice_conf::nn
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime ).
Definition time_slice_conf.C:591
Lorene::Time_slice_conf::hata
virtual const Sym_tensor & hata() const
Conformal representation of the traceless part of the extrinsic curvature: .
Definition time_slice_conf.C:772
Lorene::Time_slice_conf::psi
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
Definition time_slice_conf.C:693
Lorene::Time_slice_conf::aa
virtual Sym_tensor aa() const
Conformal representation of the traceless part of the extrinsic curvature: .
Definition time_slice_conf.C:765
Lorene::Time_slice_conf::psi4
const Scalar & psi4() const
Factor at the current time step (jtime ).
Definition time_slice_conf.C:707
Lorene::Time_slice_conf::tgam
virtual const Metric & tgam() const
Conformal metric Returns the value at the current time step (jtime ).
Definition time_slice_conf.C:747
Lorene::Time_slice::beta
virtual const Vector & beta() const
shift vector at the current time step (jtime )
Definition time_slice_access.C:87
Lorene::Tslice_dirac_max::hh
virtual const Sym_tensor & hh(Param *par_bc=0x0, Param *par_mat=0x0) const
Deviation of the conformal metric from the flat metric : .
Definition tslice_dirac_max.C:474
Lorene::Tslice_dirac_max::solve_psi
virtual Scalar solve_psi(const Scalar *ener_dens=0x0) const
Solves the elliptic equation for the conformal factor $\Psi$ (Hamiltonian constraint).
Definition tslice_dirac_max_solve.C:172
Lorene::Tslice_dirac_max::solve_npsi
virtual Scalar solve_npsi(const Scalar *ener_dens=0x0, const Scalar *trace_stress=0x0) const
Solves the elliptic equation for (maximal slicing condition + Hamiltonian constraint).
Definition tslice_dirac_max_solve.C:118
Lorene::Tslice_dirac_max::solve_beta
virtual Vector solve_beta(int method=6) const
Solves the elliptic equation for the shift vector from (Eq.
Definition tslice_dirac_max_solve.C:224
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::poisson
Vector poisson(double lambda, int method=6) const
Solves the vector Poisson equation with *this as a source.
Definition vector_poisson.C:518
Lorene::Vector::divergence
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene::Tensor_sym::derive_cov
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor_sym_calculus.C:192
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::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::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::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::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::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:139
Unites
Standard units of space, time and mass.