Excised_slice/isol_hole_compute_metric.C
/*
* Method of class Isol_Hole to compute metric data associated to a quasistationary single
* black hole spacetime slice.
*
* (see file isol_hole.h for documentation).
*
*/
/*
* Copyright (c) 2009 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 Lorene
#include "param_elliptic.h"
#include "proto.h"
#include "excised_slice.h"
#include "unites.h"
namespace Lorene {
void Excised_slice::compute_stat_metric(double precis, double Omega, bool NorKappa,
Scalar boundNoK, bool isCF,double relax, int mer_max,
int mer_max2, bool isvoid) {
// Fundamental constants and units
// -------------------------------
using namespace Unites ;
// Verifying that we are in the right case
assert(type_hor==1);
// Noticing the user that the iteration has started
cout << "================================================================" << endl;
cout << "STARTING THE MAIN ITERATION FOR COMPUTING METRIC FIELDS" << endl;
cout << " iteration parameters are the following: " << endl;
cout << " convergence precision required:" << precis << endl;
cout << " max number of global steps :" << mer_max << endl;
cout << " relaxation parameter :" << relax << endl;
cout << "================================================================" << endl;
// Construction of a multi-grid (Mg3d) and an affine mapping from the class mapping
// --------------------------------------------------------------------------------
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 int nz = (*mgrid).get_nzone(); // Number of domains
int nt = (*mgrid).get_nt(1); // Number of collocation points in theta in each domain
const int np = (*mgrid).get_np(1);
const Coord& rr = (*map).r;
Scalar rrr (*map) ;
rrr = rr ;
rrr.std_spectral_base();
// 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.
assert((rrr.val_grid_point(1,0,0,0) - 1.) <= 1.e-9);
// Angular mapping defined as well
//--------------------------------
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); //2D mapping; check if this is useful.
const Metric_flat& mets = (*map).flat_met_spher() ;
//----------------
// Initializations
// ---------------
Scalar logn (*map) ; logn = log(lapse) ;
logn.std_spectral_base();
Scalar logpsi(*map) ; logpsi = log(conf_fact) ;
logpsi.std_spectral_base();
Scalar psi4 (*map) ; psi4 = conf_fact*conf_fact*conf_fact*conf_fact ;
Scalar npsi (*map) ; npsi =conf_fact*lapse ;
Scalar conf_fact_new(*map) ; conf_fact_new.annule_hard(); conf_fact_new.std_spectral_base();
Scalar npsi_new(*map); npsi_new.annule_hard(); npsi_new.std_spectral_base();
Vector shift_new (*map, CON, (*map).get_bvect_spher());
for(int i=1; i<=3; i++){
shift_new.set(i)=0;
}
shift_new.std_spectral_base();
// Non conformally flat variables
//-------------------------------
Sym_tensor gamtuu = mets.con() + hij;
Metric gamt(gamtuu);
Metric gam(gamt.cov()*psi4) ;
Sym_tensor gamma = gam.cov();
// Extrinsic curvature variables
//------------------------------
Sym_tensor aa(*map, CON, (*map).get_bvect_spher());
for (int iii= 1; iii<=3; iii++){
for(int j=1; j<=3; j++){
aa.set(iii,j)= 0;
}
}
aa.std_spectral_base();
Scalar aa_quad_scal(*map) ; aa_quad_scal = 0. ;
Sym_tensor aa_hat(*map, CON, (*map).get_bvect_spher());
for (int iii= 1; iii<=3; iii++){
for(int j=1; j<=3; j++){
aa_hat.set(iii,j)= 0;
}
}
Sym_tensor kuu = aa/psi4 ;
Sym_tensor kuu2 = aa_hat/(psi4*psi4*sqrt(psi4));
Sym_tensor kdd = contract (gamma, 0, contract(gamma, 1, kuu, 0),1);
// (2,1)-rank delta tensor: difference between ricci rotation coefficients.
//-------------------------------------------------------------------------
Tensor delta = -0.5*contract( hij, 1, gamt.cov().derive_cov(mets), 2);
Scalar tmp(*map);
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
for (int k=1; k<=3; k++) {
tmp = 0.;
tmp = -0.5 *(gamt.cov().derive_con(mets))(i,j,k);
for (int l=1; l<=3; l++) {
tmp += -0.5*( gamt.cov()(i,l)*(hij.derive_cov(mets))(k,l,j)
+ gamt.cov()(l,j)*(hij.derive_cov(mets))(k,l,i));
}
delta.set(k,i,j) += tmp ;
}
}
}
// Conformal Rstar scalar(eq 61, Bonazzola et al. 2003)
Scalar Rstar =
0.25 * contract(gamt.con(), 0, 1,
contract(hij.derive_cov(mets), 0, 1, gamt.cov().derive_cov(mets), 0, 1), 0, 1 )
- 0.5 * contract(gamt.con(), 0, 1,
contract(hij.derive_cov(mets), 0, 1, gamt.cov().derive_cov(mets), 0, 2), 0, 1 ) ;
Scalar norm(*map);
Scalar norm3(*map);
if (isvoid == false){
cout <<"FAIL: case of non-void spacetime not treated yet" << endl;
}
else {
// Parameters for the iteration
//-----------------------------
double diff_ent = 1 ; // initialization of difference marker between two iterations;
int util = 0; // Tool used to stop tensorial iteration at any wished step "util"
for(int mer=0 ;(diff_ent > precis) && (mer<mer_max) ; mer++) {
//Global relaxation coefficient
// Scalar variables linked to the norm of normal vector to horizon.
norm = sqrt(1. + hij(1,1)); norm.std_spectral_base();
norm3 = sqrt(1. + hij(3,3)); norm3.std_spectral_base();
// Solving for (Psi-1) //
// Setting of the boundary
//------------------------
double diric= 0.;
double neum = 1.;
Vector ssalt = rrr.derive_cov(gam);
Vector ssaltcon = ssalt.up_down(gam);
Scalar ssnormalt = sqrt(contract (ssalt,0, ssaltcon, 0));
ssnormalt.std_spectral_base();
ssalt.annule_domain(nz-1);
ssalt.annule_domain(0);
ssaltcon.annule_domain(nz-1);
ssaltcon.annule_domain(0);
ssalt = ssalt/ssnormalt;
ssaltcon = ssaltcon/ssnormalt;
// \tilde{s} in the notations of Gourgoulhon and Jaramillo, 2006
Vector ssconalt = ssaltcon*conf_fact*conf_fact;
ssconalt.std_spectral_base();
ssconalt.annule_domain(nz-1);
Scalar bound3bis = -((1./conf_fact)*contract((contract(kdd,1,ssconalt,0)),0, ssconalt,0));
bound3bis.annule_domain(nz-1);
bound3bis += -conf_fact*ssconalt.divergence(gamt);
bound3bis.annule_domain(nz-1);
bound3bis = 0.25*bound3bis;
bound3bis += -contract(conf_fact.derive_cov(gamt),0,ssconalt,0) + conf_fact.dsdr();
bound3bis.annule_domain(nz-1);
bound3bis.std_spectral_base();
bound3bis.set_spectral_va().ylm();
Mtbl_cf *boundd3bis = bound3bis.set_spectral_va().c_cf;
// Computing the source
//---------------------
Scalar source_conf_fact(*map) ; source_conf_fact=3. ; // Pour le fun...
source_conf_fact.std_spectral_base();
Scalar d2logpsi = contract(conf_fact.derive_cov(mets).derive_cov(mets), 0, 1, hij, 0,1);
d2logpsi.inc_dzpuis(1);
source_conf_fact = -(0.125* aa_quad_scal )/(psi4*conf_fact*conf_fact*conf_fact)
+ conf_fact* 0.125* Rstar - d2logpsi;
source_conf_fact.std_spectral_base();
if (source_conf_fact.get_etat() == ETATZERO) {
source_conf_fact.annule_hard() ;
source_conf_fact.set_dzpuis(4) ;
source_conf_fact.std_spectral_base() ;
}
source_conf_fact.set_spectral_va().ylm();
// System inversion
//-----------------
Param_elliptic source11(source_conf_fact);
// Resolution has been done for quantity Q-1,
// because our solver gives a vanishing solution at infinity!
conf_fact_new =
source_conf_fact.sol_elliptic_boundary(source11, *boundd3bis, diric , neum) + 1 ;
// tests for resolution
//---------------------
Scalar baba2 = (conf_fact_new-1).laplacian();
// cout << "psi+1-resolution" << endl;
// maxabs (baba2 - source_conf_fact);
Scalar psinewbis = conf_fact_new -1. ; psinewbis.annule_domain(nz -1);
psinewbis.std_spectral_base();
psinewbis = psinewbis.dsdr();
Scalar psinewfin2 (map_2) ;
psinewfin2.allocate_all();
psinewfin2.set_etat_qcq();
psinewfin2.std_spectral_base();
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
psinewfin2.set_grid_point(0, k, j, 0) =
psinewbis.val_grid_point(1, k,j,0) - bound3bis.val_grid_point(1, k, j, 0);
}
// maxabs (psinewfin2);
// Update during the loop
//-----------------------
conf_fact = conf_fact_new* (1-relax) + conf_fact* relax ;
psi4 = conf_fact*conf_fact*conf_fact*conf_fact;
logpsi = log(conf_fact) ;
logpsi.std_spectral_base();
// Solving for (N*Psi -1)/
// Setting of the boundary
//------------------------
assert (NorKappa == false) ;
Scalar bound(*map);
bound = (boundNoK)*conf_fact -1;
bound.annule_domain(nz -1);
bound.std_spectral_base();
bound.set_spectral_va().ylm();
Mtbl_cf *boundd = bound.get_spectral_va().c_cf;
diric =1;
neum = 0 ;
// Computing the source ...
//-------------------------
Scalar d2lognpsi = contract(npsi.derive_cov(mets).derive_cov(mets), 0, 1, hij, 0,1);
d2lognpsi.inc_dzpuis(1); // dzpuis correction.
Scalar source_npsi = npsi*(aa_quad_scal*(7./8.)/(psi4*psi4) + Rstar/8.) - d2lognpsi;
source_npsi.std_spectral_base();
if (source_npsi.get_etat() == ETATZERO) {
source_npsi.annule_hard() ;
source_npsi.set_dzpuis(4) ;
source_npsi.std_spectral_base() ;
}
// Inversion of the operator
//--------------------------
Param_elliptic source1 (source_npsi);
npsi_new = source_npsi.sol_elliptic_boundary(source1, *boundd, diric, neum) ;
npsi_new = npsi_new +1;
// Resolution tests in npsi
//-------------------------
Scalar baba = npsi_new.laplacian();
// cout << "resolution_npsi" << endl;
// maxabs (baba - source_npsi);
// cout << "bound_npsi" << endl;
Scalar npsibound2 (map_2) ;
npsibound2.allocate_all();
npsibound2.set_etat_qcq();
npsibound2.std_spectral_base();
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
npsibound2.set_grid_point(0, k, j, 0) =
npsi_new.val_grid_point(1, k,j,0) - bound.val_grid_point(1, k, j, 0) -1.;
}
// maxabs (npsibound2);
// Update during the loop
//-----------------------
npsi = npsi_new*(1-relax) + npsi* relax;
lapse = npsi/conf_fact;
logn = log(lapse);
logn.std_spectral_base();
//Resolution in Beta //
// Setting of the boundary
//------------------------
bound = (boundNoK)/(conf_fact*conf_fact) ;
bound.annule_domain(nz -1);
// Rotation parameter for spacetime
Scalar hor_rot(*map); hor_rot.annule_hard(); hor_rot = Omega;
hor_rot.std_spectral_base() ; hor_rot.mult_rsint();
hor_rot.annule_domain(nz -1);
Vector limit = shift_new;
Vector ephi(*map, CON, (*map).get_bvect_spher());
ephi.set(1).annule_hard();
ephi.set(2).annule_hard();
ephi.set(3) = 1;
ephi.std_spectral_base();
ephi.annule_domain(nz -1);
limit = bound*ssconalt + hor_rot*ephi;
// Boundary is fixed by value of 3 components of a vector (rather than value of potentials)
limit.std_spectral_base();
Scalar Vrb = limit(1); Vrb.set_spectral_va().ylm();
Scalar mmub = limit.mu(); mmub.set_spectral_va().ylm();
Scalar etab = limit.eta(); etab.set_spectral_va().ylm();
// Computing the source
//---------------------
Vector deltaA = - 2*lapse*contract(delta, 1,2, aa, 0,1);
Vector hijddb = - contract (shift.derive_cov(mets).derive_cov(mets), 1,2, hij, 0,1) ;
Vector hijddivb =
- 0.3333333333333* contract (shift.divergence(mets).derive_cov(mets),0, hij,1);
hijddb.inc_dzpuis(); // dzpuis fixing patch...
hijddivb.inc_dzpuis();
Vector sourcevect2(*map,CON, (*map).get_bvect_spher());
sourcevect2 = 2.* contract(aa, 1, lapse.derive_cov(mets),0)
- 12*lapse*contract(aa, 1, logpsi.derive_cov(mets), 0)
+ deltaA + hijddb + hijddivb ;
sourcevect2.set(1).set_dzpuis(4);
sourcevect2.set(2).set_dzpuis(4);
sourcevect2.set(3).set_dzpuis(4);
sourcevect2.std_spectral_base();
if(sourcevect2.eta().get_etat() == ETATZERO)
{ sourcevect2.set(2).annule_hard();}
double lam = (1./3.);
// System inversion
//-----------------
sourcevect2.poisson_boundary2(lam, shift_new, Vrb, etab, mmub, 1., 0., 1. ,0. ,1. ,0.) ;
// resolution tests
//-----------------
Vector source2 = contract(shift_new.derive_con(mets).derive_cov(mets), 1,2)
+ lam* contract(shift_new.derive_cov(mets), 0,1).derive_con(mets);
source2.inc_dzpuis(1);
// maxabs (source2 - sourcevect2);
Scalar mufin = shift_new.mu();
mufin.set_spectral_va().coef();
Scalar mufin2 (map_2) ;
mufin2.allocate_all();
mufin2.set_etat_qcq();
mufin2.std_spectral_base();
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
mufin2.set_grid_point(0, k, j, 0) =
mufin.val_grid_point(1, k,j,0) - mmub.val_grid_point(1, k, j, 0);
}
// maxabs (mufin2);
Scalar brfin = shift_new(1);
brfin.set_spectral_va().coef();
Scalar brfin2 (map_2) ;
brfin2.allocate_all();
brfin2.set_etat_qcq();
brfin2.std_spectral_base();
for (int k=0; k<np; k++)
for (int j=0; j<nt; j++) {
brfin2.set_grid_point(0, k, j, 0) =
brfin.val_grid_point(1, k,j,0) - Vrb.val_grid_point(1, k, j, 0);
}
// maxabs (brfin2);
// Update during the loop
//-----------------------
for (int ii=1; ii <=3; ii++){
shift.set(ii) = shift_new(ii)*(1-relax) + shift(ii)* relax;
}
diff_ent = max(maxabs(npsi_new - npsi )); // Convergence parameter (discutable relevance...)
// Tensor hij resolution ///
if (isCF == false){
if (diff_ent <=5.e-3) { // No resolution until we are close to the result.
//WARNING; parameter maybe to be changed according to convergence parameters
//in Poisson-Hole/Kerr1.9/pplncp.C
util = util+1; // Loop marker for NCF equation.
// Local convergence can be asked for hij equation.
// Allows to satisfy integrability conditions.
// Here we ask for a local convergence; this can be improved later.
// WARNING; parameter maybe to be changed according to convergence parameters
// in Poisson-Hole/Kerr1.9/pplncp.C
double precis2 = 1.e5*precis ;
double diff_ent2 = 1 ; // Local convergence marker
// Local relaxation parameter. If not 1, the determinant condition won't be satisfied
// on a particular iteration.
double relax2 = 0.5;
Sym_tensor sourcehij = hij; // Random initialization...
// cuts off high spherical harmonics with threshold being the last argument.
// coupe_l_tous (hij, aa, nn, ppsi, bb, nt, 6);
for(int mer2=0 ;(diff_ent2 > precis2) && (mer2<mer_max2) ; mer2++) {
// Calculation of the source
//--------------------------
// The double Lie derivative term is taken care of in the subroutine.
secmembre_kerr(sourcehij);
//System inversion (note that no boundary condition is imposed)
//-------------------------------------------------------------
Sym_tensor hij_new = hij;
hij_new = boundfree_tensBC (sourcehij, shift , conf_fact, lapse, hij, precis2);
cout << "maximum of convergence marker for the subiteration" << endl;
diff_ent2 = max(maxabs(hij - hij_new));
hij = relax2*hij_new + (1 - relax2)*hij;
cout << "mer2, diffent2" << endl;
cout << mer2 << endl;
cout << diff_ent2 << endl;
}
// Resolution tests
//-----------------
Sym_tensor gammatilde = mets.con() + hij;
Metric gammatilde2(gammatilde); Scalar detgam = gammatilde2.determinant();
// cout << "determinant of result" << endl;
// maxabs (detgam-1.);
// cout << "comment l'equation en hij est elle verifiee?" << endl;
Sym_tensor test =contract (hij.derive_cov(mets).derive_con(mets), 2,3);
test.annule(nz-1, nz-1);
test = test
- hij.derive_lie(shift).derive_lie(shift)
/ ((lapse/(conf_fact*conf_fact))*(lapse/(conf_fact*conf_fact)));
test.annule(nz-1, nz-1);
Sym_tensor youps = test
- sourcehij
/ ((lapse/(conf_fact*conf_fact))*(lapse/(conf_fact*conf_fact)));
// maxabs (youps);
// maxabs((youps).trace(mets));
// cout << " AAABBB" << endl;
// maxabs((youps).compute_A());
// maxabs((youps).compute_tilde_B());
}
}
// Global variable update after an entire loop ////
gamtuu = mets.con() + hij;
gamt = gamtuu; // Metric
gam = gamt.cov()*psi4;
gamma = gam.cov();
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
tmp = 0;
tmp = ((shift.derive_con(mets))(i,j) + (shift.derive_con(mets))(j,i)
- (2./3.)*shift.divergence(mets) * mets.con()(i,j))*(1./(2.*lapse));
aa.set(i,j) = tmp ;
}
}
//Non conformally flat correction; we assume here dhij/dt = 0.
aa = aa - (hij.derive_lie(shift) + (2./3.)*shift.divergence(mets)*hij)*(1./(2.*lapse));
aa_hat = aa*psi4*sqrt(psi4); // Rescaling of traceless exrinsic curvature.
aa_hat.std_spectral_base();
hatA = aa_hat; // Probably obsolete very soon... replace aa_hat by hatA.
Sym_tensor aaud = aa.up_down(gamt);
Sym_tensor aaud_hat = aa_hat.up_down(gamt);
aa_quad_scal = contract(contract (aa_hat, 0, aaud_hat, 0), 0,1);
delta = -0.5*contract( hij, 1, gamt.cov().derive_cov(mets), 2);
for (int i=1; i<=3; i++) {
for (int j=1; j<=3; j++) {
for (int k=1; k<=3; k++) {
tmp = 0.;
tmp = -0.5 *(gamt.cov().derive_con(mets))(i,j,k);
for (int l=1; l<=3; l++) {
tmp += -0.5*( gamt.cov()(i,l)*(hij.derive_cov(mets))(k,l,j)
+ gamt.cov()(l,j)*(hij.derive_cov(mets))(k,l,i));
}
delta.set(k,i,j) += tmp ;
}
}
}
Rstar =
0.25 * contract(gamt.con(), 0, 1,
contract(gamt.con().derive_cov(mets), 0, 1,
gamt.cov().derive_cov(mets), 0, 1), 0, 1 )
- 0.5 * contract(gamt.con(), 0, 1,
contract(gamt.con().derive_cov(mets), 0, 1,
gamt.cov().derive_cov(mets), 0, 2), 0, 1 ) ;
kuu = aa/(psi4);
kdd = contract (gamma, 0, contract(gamma, 1, kuu, 0),1);
// Convergence markers
//--------------------
cout << "diffent" << endl;
cout<< diff_ent << endl;
cout <<"mer" << mer << endl;
//------------------------------------------------
// Check of Einstein equations (3+1 form)
//------------------------------------------------
// Lapse
//------
Scalar lapse2(*map) ;
lapse2 = lapse ;
lapse2.std_spectral_base();
// 3-metric
//---------
// const Metric gam(mets.cov()*psi4) ;
Sym_tensor gamt2(*map, COV, (*map).get_bvect_spher());
for (int i=1; i<=3; i++)
for (int j=1; j<=3; j++)
{ gamt2.set(i,j)=gam.cov()(i,j);
}
//Shift
//-----
Vector beta = shift ;
// Extrinsic curvature
//--------------------
Scalar TrK3(*map);
Sym_tensor k_uu = aa/(psi4) ;
k_uu.dec_dzpuis(k_uu(1,1).get_dzpuis());
// Another way of computing the same thing, just to be sure...
Sym_tensor k_dd = k_uu.up_down(gam);
TrK3 = k_uu.trace(gam);
// TrK3.spectral_display("TraceKvraie", 1.e-10);
// Hamiltonian constraint
//-----------------------
Scalar ham_constr = gam.ricci_scal() ;
ham_constr.dec_dzpuis(3) ;
ham_constr += TrK3*TrK3 - contract(k_uu, 0, 1, k_dd, 0, 1) ;
// maxabs(ham_constr, "Hamiltonian constraint: ") ;
ham_constr.set_spectral_va().ylm();
// ham_constr.spectral_display("ham_constr", 1.e-9);
// Momentum constraint
//-------------------
Vector mom_constr = k_uu.divergence(gam) - TrK3.derive_con(gam) ;
mom_constr.dec_dzpuis(2) ;
// maxabs(mom_constr, "Momentum constraint: ") ;
// mom_constr(1).spectral_display("mom1", 1.e-9) ;
// mom_constr(2).spectral_display("mom2", 1.e-9) ;
// mom_constr(3).spectral_display("mom3", 1.e-9) ;
// Evolution equations
//--------------------
Sym_tensor evol_eq = lapse2*gam.ricci() - lapse2.derive_cov(gam).derive_cov(gam);
evol_eq.dec_dzpuis() ;
evol_eq += k_dd.derive_lie(beta) ;
evol_eq.dec_dzpuis(2) ;
evol_eq += lapse2*(TrK3*k_dd - 2*contract(k_dd, 1, k_dd.up(0, gam), 0) ) ;
// maxabs(evol_eq, "Evolution equations: ") ;
// evol_eq.trace(gam).spectral_display("evoltrace", 1.e-10);
// maxabs (evol_eq.trace(gam));
// evol_eq.spectral_display("evol", 1.e-10);
}
cout << "================================================================" << endl;
cout << " THE ITERATION HAS NOW CONVERGED" << endl;
cout << "================================================================" << endl;
}
return;
}
}
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::compute_stat_metric
void compute_stat_metric(double precis, double Omega_i, bool NorKappa_i, Scalar NoK_i, bool isCF_i=true, double relax=0.2, int mer_max=2000, int mer_max2=200, bool isvoid=true)
If hor_type=1, computes a quasi-stationary 3-slice from the chosen parameters.
Definition Excised_slice/isol_hole_compute_metric.C:38
Lorene::Excised_slice::type_hor
int type_hor
Chosen horizon type.
Definition excised_slice.h:52
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::mp
const Map & mp
Mapping associated with the slice.
Definition excised_slice.h:46
Lorene::Excised_slice::shift
Vector shift
Shift vector.
Definition excised_slice.h:64
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::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
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::determinant
virtual const Scalar & determinant() const
Returns the determinant.
Definition metric.C:392
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::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Lorene::Param_elliptic
This class contains the parameters needed to call the general elliptic solver.
Definition param_elliptic.h:135
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::laplacian
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition scalar_deriv.C:436
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
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::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::dsdr
const Scalar & dsdr() const
Returns of *this .
Definition scalar_deriv.C:113
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:604
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version).
Definition scalar.h:601
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
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::allocate_all
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition scalar.C:365
Lorene::Scalar::mult_rsint
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
Definition scalar_r_manip.C:281
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Lorene::Scalar::sol_elliptic_boundary
Scalar sol_elliptic_boundary(Param_elliptic &params, const Mtbl_cf &bound, double fact_dir, double fact_neu) const
Resolution of a general elliptic equation, putting zero at infinity and with inner boundary condition...
Definition scalar_pde.C:256
Lorene::Scalar::set_grid_point
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition scalar.h:684
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::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::Tenseur::derive_con
const Tenseur & derive_con(const Metrique &) const
Returns the contravariant derivative of *this , with respect to met .
Definition tenseur.C:1568
Lorene::Tensor
Tensor handling.
Definition tensor.h:288
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition valeur_coef.C:148
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::mu
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written:
Definition vector_etamu.C:98
Lorene::Vector::eta
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written:
Definition vector_etamu.C:66
Lorene::Vector::poisson_boundary2
void poisson_boundary2(double lam, Vector &resu, Scalar boundvr, Scalar boundeta, Scalar boundmu, double dir_vr, double neum_vr, double dir_eta, double neum_eta, double dir_mu, double neum_mu) const
Alternative to previous poisson_boundary method for vectors ; this uses method 6 for vectorial solvin...
Definition vector_poisson_boundary2.C:80
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::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
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::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::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::annule
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
Definition tensor.C:671
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::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::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::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::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
Unites
Standard units of space, time and mass.