Excised_slice/secmembre_hij_stat.C
// Header Lorene:
#include "param_elliptic.h"
#include "excised_slice.h"
#include "unites.h"
// Computes the rhs of hyperbolic equation for conformal metric assuming statioarity;
// WARNING; up to now, we are only able to handle void spacetimes.
namespace Lorene {
void Excised_slice::secmembre_kerr(Sym_tensor& source_hh){
// Getting: hij; hatA; lapse; conf_fact; shift;
// hij; aa; nn; ppsi; bb;
using namespace Unites;
const int nz = (*(hij.get_mp().get_mg())).get_nzone();
const Vector& beta = shift;
const Sym_tensor& hh = hij;
const Scalar& psi4 = conf_fact*conf_fact*conf_fact*conf_fact;
Scalar ln_psi = log(conf_fact); ln_psi.std_spectral_base();
const Scalar qq = lapse*conf_fact*conf_fact;
Sym_tensor aa = hatA/(psi4*sqrt(psi4));
aa.std_spectral_base(); //(check...)
const Metric_flat& ff = (hij.get_mp()).flat_met_spher() ;
const Sym_tensor& tgam_uu = ff.con() + hh;
const Metric tgam(tgam_uu);
const Base_vect_spher& otriad = hij.get_mp().get_bvect_spher();
// On met a zero les quantites supposees etre de "matiere"
Sym_tensor strain_tens = hij;
for (int ii=1; ii<=3; ii++)
for(int jj=1; jj<=3; jj++)
{ strain_tens.set(ii,jj).annule_hard();
}
strain_tens.std_spectral_base();
// On met a zero les quantites derivee temporelle
Vector beta_point = shift;
for (int ii=1; ii<=3; ii++)
{ beta_point.set(ii).annule_hard();
}
beta_point.annule_domain(nz-1) ; // Pour faire passer un assert, je ne comprends pas pourquoi...
beta_point.std_spectral_base();
Scalar lapse_point(hij.get_mp());
lapse_point.annule_hard();
lapse_point.annule_domain(nz -1);
lapse_point.std_spectral_base();
Sym_tensor hh_point = hij;
for (int ii=1; ii<=3; ii++)
for(int jj=1; jj<=3; jj++)
{ hh_point.set(ii,jj).annule_hard();
}
hh_point.annule_domain(nz-1);
hh_point.std_spectral_base();
// Note: Il sera probablement necessaire de ne pas mettre a zero hh point;
//Sym_tensor Rrij(map, CON, map.get_bvect_spher());
// Rrij = 0.5*[ contract (( h_iju, 0, h_iju.derive_cov(mets).derive_cov(mets), 3), 0,3) - contract((contract(h_iju.derive_cov(mets),1, h_iju.derive_cov(mets),2)), 1,3)] ;
//==================================
// Source for hij
//==================================
const Sym_tensor& tgam_dd = tgam.cov() ; // {\tilde \gamma}_{ij}
// const Sym_tensor& tgam_uu = tgam().con() ; // {\tilde \gamma}^{ij}
const Tensor_sym& dtgam = tgam_dd.derive_cov(ff) ;// D_k {\tilde \gamma}_{ij} // ff etant la metrique plate
const Tensor_sym& dhh = hh.derive_cov(ff) ; // D_k h^{ij}
const Vector& dln_psi = ln_psi.derive_cov(ff) ; // D_i ln(Psi)
const Vector& tdln_psi_u = ln_psi.derive_con(tgam) ; // tD^i ln(Psi)
const Vector& tdnn_u = lapse.derive_con(tgam) ; // tD^i N
const Vector& dqq = qq.derive_cov(ff) ; // D_i Q
const Scalar& div_beta = beta.divergence(ff) ; // D_k beta^k
Scalar tmp(hij.get_mp()) ;
Sym_tensor sym_tmp(hij.get_mp(), CON, otriad) ;
// Quadratic part of the Ricci tensor of gam_tilde
// ------------------------------------------------
Sym_tensor ricci_star(hij.get_mp(), CON, otriad) ;
ricci_star = contract(hh, 0, 1, dhh.derive_cov(ff), 2, 3) ;
ricci_star.inc_dzpuis() ; // dzpuis : 3 --> 4
// des_profile (ricci_star(1,1), 1, 8, 1,1, "riccistar"); // A enlever
for (int i=1; i<=3; i++) {
for (int j=1; j<=i; j++) {
tmp = 0 ;
for (int k=1; k<=3; k++) {
for (int l=1; l<=3; l++) {
tmp += dhh(i,k,l) * dhh(j,l,k) ;
}
}
sym_tmp.set(i,j) = tmp ;
}
}
ricci_star -= sym_tmp ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=i; j++) {
tmp = 0 ;
for (int k=1; k<=3; k++) {
for (int l=1; l<=3; l++) {
for (int m=1; m<=3; m++) {
for (int n=1; n<=3; n++) {
tmp += 0.5 * tgam_uu(i,k)* tgam_uu(j,l)
* dhh(m,n,k) * dtgam(m,n,l)
+ tgam_dd(n,l) * dhh(m,n,k)
* (tgam_uu(i,k) * dhh(j,l,m) + tgam_uu(j,k) * dhh(i,l,m) )
- tgam_dd(k,l) *tgam_uu(m,n) * dhh(i,k,m) * dhh(j,l,n) ;
}
}
}
}
sym_tmp.set(i,j) = tmp ;
}
}
ricci_star += sym_tmp ;
ricci_star = 0.5 * ricci_star ;
// Curvature scalar of conformal metric :
// -------------------------------------
Scalar tricci_scal =
0.25 * contract(tgam_uu, 0, 1,
contract(dhh, 0, 1, dtgam, 0, 1), 0, 1 )
- 0.5 * contract(tgam_uu, 0, 1,
contract(dhh, 0, 1, dtgam, 0, 2), 0, 1 ) ;
// Full quadratic part of source for h : S^{ij}
// --------------------------------------------
Sym_tensor ss(hij.get_mp(), CON, otriad) ; // Source secondaire
sym_tmp = lapse * (ricci_star + 8.* tdln_psi_u * tdln_psi_u)
+ 4.*( tdln_psi_u * tdnn_u + tdnn_u * tdln_psi_u )
- 0.3333333333333333 *
( lapse * (tricci_scal + 8.* contract(dln_psi, 0, tdln_psi_u, 0) )
+ 8.* contract(dln_psi, 0, tdnn_u, 0) ) *tgam_uu ;
ss = sym_tmp / psi4 ;
sym_tmp = contract(tgam_uu, 1,
contract(tgam_uu, 1, dqq.derive_cov(ff), 0), 1) ;
sym_tmp.inc_dzpuis() ; // dzpuis : 3 --> 4
// des_profile (sym_tmp(1,1), 1, 8, 1,1, "sym_tmp"); // A enlever
for (int i=1; i<=3; i++) {
for (int j=1; j<=i; j++) {
tmp = 0 ;
for (int k=1; k<=3; k++) {
for (int l=1; l<=3; l++) {
tmp += ( hh(i,k)*dhh(l,j,k) + hh(k,j)*dhh(i,l,k)
- hh(k,l)*dhh(i,j,k) ) * dqq(l) ;
}
}
sym_tmp.set(i,j) += 0.5 * tmp ;
}
}
tmp = qq.derive_con(tgam).divergence(tgam) ;
tmp.inc_dzpuis() ; // dzpuis : 3 --> 4 // reverifier pourquoi
// des_profile (tmp, 1, 8, 1,1, "tmp"); // A enlever
sym_tmp -= 0.3333333333333333 * tmp *tgam_uu ;
ss -= sym_tmp / (psi4*conf_fact*conf_fact) ; // Voir dans quel sens sont construits psi et psi4 (eviter les multiplications d'erreurs)
for (int i=1; i<=3; i++) {
for (int j=1; j<=i; j++) {
tmp = 0 ;
for (int k=1; k<=3; k++) {
for (int l=1; l<=3; l++) {
tmp += tgam_dd(k,l) * aa(i,k) * aa(j,l) ;
}
}
sym_tmp.set(i,j) = tmp ;
}
}
tmp = psi4 * strain_tens.trace(tgam) ; // S = S_i^i
ss += (2.*lapse) * ( sym_tmp);// - qpig*( psi4* strain_tens
// - 0.3333333333333333 * tmp * tgam_uu );
Sym_tensor ss2 =2.*lapse*( qpig*(psi4*strain_tens - 0.33333333333333 * tmp * tgam_uu));
ss2.inc_dzpuis(4); // A retirer!
// des_profile (ss2(1,1), 1, 8, 1,1, "ss2"); // A enlever
// cout << zone << endl;
// ss2.annule_domain(nz-1);
ss += -ss2; // ATTENTION!!!! A RETABLIR!!!!
// maxabs(ss, "ss tot") ;
// Source for h^{ij}
// -----------------
Sym_tensor lbh = hh.derive_lie(beta) ;
source_hh =// (lapse*lapse/psi4 - 1.)
// * hh.derive_con(ff).divergence(ff)
+ 2.* hh_point.derive_lie(beta); // - lbh.derive_lie(beta) ; // La double derivee de
// Lie en Beta est retiree (prise en charge dans tensorelliptic.C)
source_hh.inc_dzpuis() ;
// des_profile (source_hh(1,1), 1, 8, 1,1, "sourcehh"); // A enlever
source_hh += 2.* lapse * ss ;
//## Provisory: waiting for the Lie derivative to allow
// derivation with respect to a vector with dzpuis != 0
Vector vtmp = beta_point ;
// vtmp.dec_dzpuis(2) ; // A remettre si jamais beta point est non nul;
sym_tmp = hh.derive_lie(vtmp) ;
sym_tmp.inc_dzpuis(2) ;
// des_profile (sym_tmp(1,1), 1, 8, 1,1, "sym_tmp"); // A enlever
source_hh += sym_tmp
+ 1.3333333333333333 * div_beta* (hh_point - lbh)
+ 2. * (lapse_point - lapse.derive_lie(beta)) * aa ;
for (int i=1; i<=3; i++) {
for (int j=1; j<=i; j++) {
tmp = 0. ;
for (int k=1; k<=3; k++) {
tmp += ( hh.derive_con(ff)(k,j,i)
+ hh.derive_con(ff)(i,k,j)
- hh.derive_con(ff)(i,j,k) ) * dqq(k) ;
}
sym_tmp.set(i,j) = tmp ;
}
}
source_hh -= lapse / (psi4*conf_fact*conf_fact) * sym_tmp ;
tmp = beta_point.divergence(ff) - div_beta.derive_lie(beta) ;
tmp.inc_dzpuis() ;
// des_profile (tmp, 1, 8, 1,1, "tmp"); // A enlever
source_hh += 0.6666666666666666*
( tmp - 0.6666666666666666* div_beta * div_beta ) * hh ;
// Term (d/dt - Lie_beta) (L beta)^{ij}--> sym_tmp
// ---------------------------------------
Sym_tensor l_beta = beta.ope_killing_conf(ff) ; // Attention aux headers a inclure
sym_tmp = beta_point.ope_killing_conf(ff) - l_beta.derive_lie(beta) ;
sym_tmp.inc_dzpuis() ;
// Final source:
// ------------
source_hh += 0.6666666666666666* div_beta * l_beta - sym_tmp ;
// Invert it (because the right source is for a Laplace operator, not a Dalembertain operator as it is initially:
source_hh = -source_hh;
// Annulation de la source
// for (int ii=1; ii<=3; ii++)
// for (int jj=1; jj<=3; jj++){
// source_hh.set(ii,jj).annule_hard();
// }
// source_hh.std_spectral_base();
return;
}
}
Lorene::Base_vect_spher
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
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::secmembre_kerr
void secmembre_kerr(Sym_tensor &source_hh)
If hor_type=1, computes the rhs of hyperbolic equation for conformal metric assuming stationarity; WA...
Definition Excised_slice/secmembre_hij_stat.C:12
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::lapse
Scalar lapse
Lapse function.
Definition excised_slice.h:58
Lorene::Excised_slice::shift
Vector shift
Shift vector.
Definition excised_slice.h:64
Lorene::Metric_flat
Flat metric for tensor calculation.
Definition metric.h:261
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::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::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::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
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::derive_lie
Scalar derive_lie(const Vector &v) const
Computes the derivative of this along a vector field v.
Definition scalar_deriv.C:419
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::Tensor_sym
Symmetric tensors (with respect to two of their arguments).
Definition tensor.h:1037
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::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::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::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
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::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_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::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
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::trace
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Definition tensor_calculus.C:94
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::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Unites
Standard units of space, time and mass.