tslice_conf_init.C
/*
* Method of class Time_slice_conf to compute valid initial data
*
* (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_conf_init_C[] = "$Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_conf_init.C,v 1.13 2014/10/13 08:53:48 j_novak Exp $" ;
/*
* $Id: tslice_conf_init.C,v 1.13 2014/10/13 08:53:48 j_novak Exp $
* $Log: tslice_conf_init.C,v $
* Revision 1.13 2014/10/13 08:53:48 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.12 2014/10/06 15:13:22 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.11 2010/10/20 07:58:09 j_novak
* Better implementation of the explicit time-integration. Not fully-tested yet.
*
* Revision 1.10 2008/12/04 18:22:49 j_novak
* Enhancement of the dzpuis treatment + various bug fixes.
*
* Revision 1.9 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.8 2004/05/17 19:53:13 e_gourgoulhon
* Added arguments graph_device and method_poisson_vect.
*
* Revision 1.7 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.6 2004/05/10 09:12:01 e_gourgoulhon
* Added a call to del_deriv() at the end.
*
* Revision 1.5 2004/05/03 14:48:48 e_gourgoulhon
* Treatment of special cases nn_jp1.etat = ETATUN and psi_jp1.etat = ETATUN.
*
* Revision 1.4 2004/04/29 17:10:36 e_gourgoulhon
* Added argument pdt and update of depth slices at the end,
* taking into account the known time derivatives.
*
* Revision 1.3 2004/04/08 16:45:11 e_gourgoulhon
* Use of new methods set_*.
*
* Revision 1.2 2004/04/07 07:58:21 e_gourgoulhon
* Constructor as Minkowski slice: added call to std_spectral_base()
* after setting the lapse to 1.
*
* Revision 1.1 2004/04/05 21:25:37 e_gourgoulhon
* First version.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_conf_init.C,v 1.13 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 "utilitaires.h"
namespace Lorene {
void Time_slice_conf::initial_data_cts(const Sym_tensor& uu,
const Scalar& trk_in, const Scalar& trk_point,
double pdt, double precis, int method_poisson_vect,
const char* graph_device, const Scalar* p_ener_dens,
const Vector* p_mom_dens, const Scalar* p_trace_stress) {
using namespace Unites ;
// Verifications
// -------------
double tr_uu = max(maxabs(uu.trace(tgam()), "trace tgam_{ij} u^{ij}")) ;
if (tr_uu > 1.e-7) {
cerr <<
"Time_slice_conf::initial_data_cts : the trace of u^{ij} with respect\n"
<< " to the conformal metric tgam_{ij} is not zero !\n"
<< " error = " << tr_uu << endl ;
abort() ;
}
assert(trk_in.check_dzpuis(2)) ;
assert(trk_point.check_dzpuis(4)) ;
// Initialisations
// ---------------
double ttime = the_time[jtime] ;
trk_evol.update(trk_in, jtime, ttime) ;
// Reset of quantities depending on K:
k_dd_evol.downdate(jtime) ;
k_uu_evol.downdate(jtime) ;
set_hata(psi4()*psi()*psi()* uu / (2.* nn()) ) ;
const Map& map = uu.get_mp() ;
const Base_vect& triad = *(uu.get_triad()) ;
// For graphical outputs:
int ngraph0 = 10 ; // index of the first graphic device to be used
int nz = map.get_mg()->get_nzone() ;
double ray_des = 1.25 * map.val_r(nz-2, 1., 0., 0.) ; // outermost radius
// for plots
Scalar ener_dens(map) ;
if (p_ener_dens != 0x0) ener_dens = *(p_ener_dens) ;
else ener_dens.set_etat_zero() ;
Vector mom_dens(map, CON, triad) ;
if (p_mom_dens != 0x0) mom_dens = *(p_mom_dens) ;
else mom_dens.set_etat_zero() ;
Scalar trace_stress(map) ;
if (p_trace_stress != 0x0) trace_stress = *(p_trace_stress) ;
else trace_stress.set_etat_zero() ;
Scalar tmp(map) ;
Scalar source_psi(map) ;
Scalar source_nn(map) ;
Vector source_beta(map, CON, triad) ;
// Iteration
// ---------
int imax = 100 ;
for (int i=0; i<imax; i++) {
//===============================================
// Computations of sources
//===============================================
const Vector& dpsi = psi().derive_cov(ff) ; // D_i Psi
const Vector& dln_psi = ln_psi().derive_cov(ff) ; // D_i ln(Psi)
const Vector& dnn = nn().derive_cov(ff) ; // D_i N
Sym_tensor taa = aa().up_down(tgam()) ;
Scalar aa_quad = contract(taa, 0, 1, aa(), 0, 1) ;
// Source for Psi
// --------------
tmp = 0.125* psi() * tgam().ricci_scal()
- contract(hh(), 0, 1, dpsi.derive_cov(ff), 0, 1 ) ;
tmp.inc_dzpuis() ; // dzpuis : 3 -> 4
tmp -= contract(hdirac(), 0, dpsi, 0) ;
source_psi = tmp - psi()*psi4()* ( 0.5*qpig* ener_dens
+ 0.125* aa_quad
- 8.33333333333333e-2* trk()*trk() ) ;
// Source for N
// ------------
source_nn = psi4()*( nn()*( qpig* (ener_dens + trace_stress) + aa_quad
- 0.3333333333333333* trk()*trk() )
- trk_point )
- 2.* contract(dln_psi, 0, nn().derive_con(tgam()), 0)
- contract(hdirac(), 0, dnn, 0) ;
tmp = psi4()* contract(beta(), 0, trk().derive_cov(ff), 0)
- contract( hh(), 0, 1, dnn.derive_cov(ff), 0, 1 ) ;
tmp.inc_dzpuis() ; // dzpuis: 3 -> 4
source_nn += tmp ;
// Source for beta
// ---------------
source_beta = 2.* contract(aa(), 1,
dnn - 6.*nn() * dln_psi, 0) ;
source_beta += 2.* nn() * ( 2.*qpig* psi4() * mom_dens
+ 0.66666666666666666* trk().derive_con(tgam())
- contract(tgam().connect().get_delta(), 1, 2,
aa(), 0, 1) ) ;
Vector vtmp = 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)
- hdirac().derive_lie(beta())
+ uu.divergence(ff) ;
vtmp.inc_dzpuis() ; // dzpuis: 3 -> 4
source_beta -= vtmp ;
source_beta += 0.66666666666666666* beta().divergence(ff) * hdirac() ;
//=============================================
// Resolution of elliptic equations
//=============================================
// Resolution of the Poisson equation for Psi
// ------------------------------------------
Scalar psi_jp1 = source_psi.poisson() + 1. ;
if (psi_jp1.get_etat() == ETATUN) psi_jp1.std_spectral_base() ;
// Test:
maxabs(psi_jp1.laplacian() - source_psi,
"Absolute error in the resolution of the equation for Psi") ;
des_meridian(psi_jp1, 0., ray_des, "Psi", ngraph0, graph_device) ;
// Resolution of the Poisson equation for the lapse
// ------------------------------------------------
Scalar nn_jp1 = source_nn.poisson() + 1. ;
if (nn_jp1.get_etat() == ETATUN) nn_jp1.std_spectral_base() ;
// Test:
maxabs(nn_jp1.laplacian() - source_nn,
"Absolute error in the resolution of the equation for N") ;
des_meridian(nn_jp1, 0., ray_des, "N", ngraph0+1, graph_device) ;
// Resolution of the vector Poisson equation for the shift
//---------------------------------------------------------
Vector beta_jp1 = source_beta.poisson(0.3333333333333333, ff,
method_poisson_vect) ;
des_meridian(beta_jp1(1), 0., ray_des, "\\gb\\ur\\d", ngraph0+2,
graph_device) ;
des_meridian(beta_jp1(2), 0., ray_des, "\\gb\\u\\gh\\d", ngraph0+3,
graph_device) ;
des_meridian(beta_jp1(3), 0., ray_des, "\\gb\\u\\gf\\d", ngraph0+4,
graph_device) ;
// Test:
Vector test_beta = (beta_jp1.derive_con(ff)).divergence(ff)
+ 0.3333333333333333 * (beta_jp1.divergence(ff)).derive_con(ff) ;
test_beta.inc_dzpuis() ;
maxabs(test_beta - source_beta,
"Absolute error in the resolution for beta") ;
//===========================================
// Convergence control
//===========================================
double diff_psi = max( diffrel(psi(), psi_jp1) ) ;
double diff_nn = max( diffrel(nn(), nn_jp1) ) ;
double diff_beta = max( diffrel(beta(), beta_jp1) ) ;
cout << "step = " << i << " : diff_psi = " << diff_psi
<< " diff_nn = " << diff_nn
<< " diff_beta = " << diff_beta << endl ;
if ( (diff_psi < precis) && (diff_nn < precis) && (diff_beta < precis) )
break ;
//=============================================
// Updates for next step
//=============================================
set_psi_del_npsi(psi_jp1) ;
n_evol.update(nn_jp1, jtime, ttime) ;
beta_evol.update(beta_jp1, jtime, ttime) ;
// New value of A^{ij}:
Sym_tensor aa_jp1 = ( beta().ope_killing_conf(tgam()) + uu )
/ (2.* nn()) ;
set_hata( aa_jp1 / (psi4()*psi()*psi()) ) ;
}
//==================================================================
// End of iteration
//===================================================================
npsi_evol.update( n_evol[jtime]*psi_evol[jtime], jtime, ttime ) ;
A_hata() ;
B_hata() ;
// Push forward in time to enable the computation of time derivatives
// ------------------------------------------------------------------
double ttime1 = ttime ;
int jtime1 = jtime ;
for (int j=1; j < depth; j++) {
jtime1++ ;
ttime1 += pdt ;
psi_evol.update(psi_evol[jtime], jtime1, ttime1) ;
npsi_evol.update(npsi_evol[jtime], jtime1, ttime1) ;
n_evol.update(n_evol[jtime], jtime1, ttime1) ;
beta_evol.update(beta_evol[jtime], jtime1, ttime1) ;
hh_evol.update(hh_evol[jtime], jtime1, ttime1) ;
hata_evol.update(hata_evol[jtime], jtime1, ttime1) ;
A_hata_evol.update(A_hata_evol[jtime], jtime1, ttime1) ;
B_hata_evol.update(B_hata_evol[jtime], jtime1, ttime1) ;
trk_evol.update(trk_evol[jtime], jtime1, ttime1) ;
the_time.update(ttime1, jtime1, ttime1) ;
}
jtime += depth - 1 ;
// Taking into account the time derivative of h^{ij} and K :
// ---------------------------------------------------------
Sym_tensor uu0 = uu ;
uu0.dec_dzpuis(2) ; // dzpuis: 2 --> 0
for (int j=1; j < depth; j++) {
hh_evol.update(hh_evol[jtime] - j*pdt* uu0,
jtime-j, the_time[jtime-j]) ;
trk_evol.update(trk_evol[jtime] - j*pdt* trk_point,
jtime-j, the_time[jtime-j]) ;
}
// Reset of derived quantities (at the new time step jtime)
// ---------------------------
del_deriv() ;
}
}
Lorene::Metric::ricci_scal
virtual const Scalar & ricci_scal() const
Returns the Ricci scalar.
Definition metric.C:350
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::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::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::Time_slice_conf::hh_evol
Evolution_std< Sym_tensor > hh_evol
Values at successive time steps of the components .
Definition time_slice.h:530
Lorene::Time_slice_conf::hh
virtual const Sym_tensor & hh(Param *=0x0, Param *=0x0) const
Deviation of the conformal metric from the flat metric : .
Definition time_slice_conf.C:758
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::set_psi_del_npsi
virtual void set_psi_del_npsi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
Definition time_slice_conf.C:404
Lorene::Time_slice_conf::ln_psi
const Scalar & ln_psi() const
Logarithm of at the current time step (jtime ).
Definition time_slice_conf.C:719
Lorene::Time_slice_conf::hata_evol
Evolution_std< Sym_tensor > hata_evol
Values at successive time steps of the components .
Definition time_slice.h:542
Lorene::Time_slice_conf::hdirac
virtual const Vector & hdirac() const
Vector which vanishes in Dirac gauge.
Definition time_slice_conf.C:815
Lorene::Time_slice_conf::A_hata
virtual const Scalar & A_hata() const
Returns the potential A of .
Definition time_slice_conf.C:664
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::B_hata
virtual const Scalar & B_hata() const
Returns the potential of .
Definition time_slice_conf.C:678
Lorene::Time_slice_conf::del_deriv
virtual void del_deriv() const
Deletes all the derived quantities.
Definition time_slice_conf.C:347
Lorene::Time_slice_conf::npsi_evol
Evolution_std< Scalar > npsi_evol
Values at successive time steps of the factor .
Definition time_slice.h:522
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::set_hata
virtual void set_hata(const Sym_tensor &hata_in)
Sets the conformal representation of the traceless part of the extrinsic curvature: .
Definition time_slice_conf.C:534
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::psi_evol
Evolution_std< Scalar > psi_evol
Values at successive time steps of the conformal factor relating the physical metric to the conform...
Definition time_slice.h:517
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_conf::A_hata_evol
Evolution_std< Scalar > A_hata_evol
Potential A associated with the symmetric tensor .
Definition time_slice.h:547
Lorene::Time_slice_conf::ff
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition time_slice.h:507
Lorene::Time_slice_conf::B_hata_evol
Evolution_std< Scalar > B_hata_evol
Potential associated with the symmetric tensor .
Definition time_slice.h:552
Lorene::Time_slice_conf::trk
virtual const Scalar & trk() const
Trace K of the extrinsic curvature at the current time step (jtime ).
Definition time_slice_conf.C:796
Lorene::Time_slice::jtime
int jtime
Time step index of the latest slice.
Definition time_slice.h:190
Lorene::Time_slice::trk_evol
Evolution_std< Scalar > trk_evol
Values at successive time steps of the trace K of the extrinsic curvature.
Definition time_slice.h:224
Lorene::Time_slice::k_dd_evol
Evolution_std< Sym_tensor > k_dd_evol
Values at successive time steps of the covariant components of the extrinsic curvature tensor .
Definition time_slice.h:208
Lorene::Time_slice::beta
virtual const Vector & beta() const
shift vector at the current time step (jtime )
Definition time_slice_access.C:87
Lorene::Time_slice::depth
int depth
Number of stored time slices.
Definition time_slice.h:179
Lorene::Time_slice::n_evol
Evolution_std< Scalar > n_evol
Values at successive time steps of the lapse function N.
Definition time_slice.h:216
Lorene::Time_slice::the_time
Evolution_std< double > the_time
Time label of each slice.
Definition time_slice.h:193
Lorene::Time_slice::k_uu_evol
Evolution_std< Sym_tensor > k_uu_evol
Values at successive time steps of the contravariant components of the extrinsic curvature tensor .
Definition time_slice.h:213
Lorene::Time_slice::beta_evol
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Definition time_slice.h:219
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::divergence
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
Lorene::Vector::ope_killing_conf
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition vector.C:462
Lorene::Vector::derive_lie
Vector derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition vector.C:392
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
des_meridian
void des_meridian(const Scalar &uu, double r_min, double r_max, const char *nomy, int ngraph, const char *device=0x0, bool closeit=false, bool draw_bound=true)
Draws 5 profiles of a scalar field along various radial axes in two meridional planes and .
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::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