star_bin_equilibrium.C
/*
* Method of class Star_bin to compute an equilibrium configuration
*
* (see file star.h for documentation).
*/
/*
* Copyright (c) 2004 Francois Limousin
*
* 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 star_bin_equilibrium_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_bin_equilibrium.C,v 1.29 2014/10/13 08:53:38 j_novak Exp $" ;
/*
* $Id: star_bin_equilibrium.C,v 1.29 2014/10/13 08:53:38 j_novak Exp $
* $Log: star_bin_equilibrium.C,v $
* Revision 1.29 2014/10/13 08:53:38 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.28 2014/10/06 15:13:16 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.27 2006/08/01 14:26:01 f_limousin
* Display
*
* Revision 1.26 2006/05/31 09:26:04 f_limousin
* Modif. of the size of the different domains
*
* Revision 1.25 2006/04/11 14:24:44 f_limousin
* New version of the code : improvement of the computation of some
* critical sources, estimation of the dirac gauge, helical symmetry...
*
* Revision 1.24 2005/11/03 13:27:09 f_limousin
* Final version for the letter.
*
* Revision 1.23 2005/09/14 12:48:02 f_limousin
* Comment graphical outputs.
*
* Revision 1.22 2005/09/14 12:30:52 f_limousin
* Saving of fields lnq and logn in class Star.
*
* Revision 1.21 2005/09/13 19:38:31 f_limousin
* Reintroduction of the resolution of the equations in cartesian coordinates.
*
* Revision 1.20 2005/04/08 12:36:44 f_limousin
* Just to avoid warnings...
*
* Revision 1.19 2005/02/24 16:27:21 f_limousin
* Add mermax_poisson and relax_poisson in the parameters of the function.
*
* Revision 1.18 2005/02/24 16:04:13 f_limousin
* Change the name of some variables (for instance dcov_logn --> dlogn).
* Improve the resolution of the tensorial poisson equation for hh.
*
* Revision 1.17 2005/02/18 13:14:18 j_novak
* Changing of malloc/free to new/delete + suppression of some unused variables
* (trying to avoid compilation warnings).
*
* Revision 1.16 2005/02/17 17:32:53 f_limousin
* Change the name of some quantities to be consistent with other classes
* (for instance nnn is changed to nn, shift to beta, beta to lnq...)
*
* Revision 1.15 2005/02/11 18:13:47 f_limousin
* Important modification : all the poisson equations for the metric
* quantities are now solved on an affine mapping.
*
* Revision 1.14 2004/12/17 16:23:19 f_limousin
* Modif. comments.
*
* Revision 1.13 2004/06/22 12:49:12 f_limousin
* Change qq, qq_auto and qq_comp to beta, beta_auto and beta_comp.
*
* Revision 1.12 2004/05/27 12:41:00 p_grandclement
* correction of some shadowed variables
*
* Revision 1.11 2004/05/25 14:18:00 f_limousin
* Include filters
*
* Revision 1.10 2004/05/10 10:26:22 f_limousin
* Minor changes to avoid warnings in the compilation of Lorene
*
* Revision 1.9 2004/04/08 16:32:48 f_limousin
* The new variable is ln(Q) instead of Q=psi^2*N. It improves the
* convergence of the code.
*
* Revision 1.8 2004/03/25 10:29:26 j_novak
* All LORENE's units are now defined in the namespace Unites (in file unites.h).
*
* Revision 1.7 2004/03/23 09:56:09 f_limousin
* Many minor changes
*
* Revision 1.6 2004/02/27 21:16:32 e_gourgoulhon
* Function contract_desal replaced by function contract with
* argument desaliasing set to true.
*
* Revision 1.5 2004/02/27 09:51:51 f_limousin
* Many minor changes.
*
* Revision 1.4 2004/02/21 17:05:13 e_gourgoulhon
* Method Scalar::point renamed Scalar::val_grid_point.
* Method Scalar::set_point renamed Scalar::set_grid_point.
*
* Revision 1.3 2004/01/20 15:17:48 f_limousin
* First version
*
*
* $Header: /cvsroot/Lorene/C++/Source/Star/star_bin_equilibrium.C,v 1.29 2014/10/13 08:53:38 j_novak Exp $
*
*/
// C headers
#include <cmath>
// Lorene headers
#include "cmp.h"
#include "tenseur.h"
#include "metrique.h"
#include "star.h"
#include "param.h"
#include "graphique.h"
#include "utilitaires.h"
#include "tensor.h"
#include "nbr_spx.h"
#include "unites.h"
namespace Lorene {
void Star_bin::equilibrium(double ent_c, int mermax, int mermax_potvit,
int mermax_poisson, double relax_poisson,
double relax_potvit, double thres_adapt,
Tbl& diff, double om) {
// Fundamental constants and units
// -------------------------------
using namespace Unites ;
// Initializations
// ---------------
const Mg3d* mg = mp.get_mg() ;
int nz = mg->get_nzone() ; // total number of domains
// The following is required to initialize mp_prev as a Map_et:
Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
// Domain and radial indices of points at the surface of the star:
int l_b = nzet - 1 ;
int i_b = mg->get_nr(l_b) - 1 ;
int k_b ;
int j_b ;
// Value of the enthalpy defining the surface of the star
double ent_b = 0 ;
// Error indicators
// ----------------
double& diff_ent = diff.set(0) ;
double& diff_vel_pot = diff.set(1) ;
double& diff_logn = diff.set(2) ;
double& diff_lnq = diff.set(3) ;
double& diff_beta_x = diff.set(4) ;
double& diff_beta_y = diff.set(5) ;
double& diff_beta_z = diff.set(6) ;
double& diff_h11 = diff.set(7) ;
double& diff_h21 = diff.set(8) ;
double& diff_h31 = diff.set(9) ;
double& diff_h22 = diff.set(10) ;
double& diff_h32 = diff.set(11) ;
double& diff_h33 = diff.set(12) ;
// Parameters for te function Map_et::adapt
// -----------------------------------------
Param par_adapt ;
int nitermax = 100 ;
int niter ;
int adapt_flag = 1 ; // 1 = performs the full computation,
// 0 = performs only the rescaling by
// the factor alpha_r
//## int nz_search = nzet + 1 ; // Number of domains for searching the
// enthalpy
int nz_search = nzet ; // Number of domains for searching the enthalpy
// isosurfaces
double precis_secant = 1.e-14 ;
double alpha_r ;
double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
Tbl ent_limit(nz) ;
par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
// locate zeros by the secant method
par_adapt.add_int(nzet, 1) ; // number of domains where the adjustment
// to the isosurfaces of ent is to be
// performed
par_adapt.add_int(nz_search, 2) ; // number of domains to search
// the enthalpy isosurface
par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
// 0 = performs only the rescaling by
// the factor alpha_r
par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
// (theta_*, phi_*)
par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
// (theta_*, phi_*)
par_adapt.add_int_mod(niter, 0) ; // number of iterations actually used in
// the secant method
par_adapt.add_double(precis_secant, 0) ; // required absolute precision in
// the determination of zeros by
// the secant method
par_adapt.add_double(reg_map, 1) ; // 1. = regular mapping,
// 0 = contracting mapping
par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
// distances will be multiplied
par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
// to define the isosurfaces.
Cmp ssjm1logn (ssjm1_logn) ;
Cmp ssjm1lnq (ssjm1_lnq) ;
Cmp ssjm1h11 (ssjm1_h11) ;
Cmp ssjm1h21 (ssjm1_h21) ;
Cmp ssjm1h31 (ssjm1_h31) ;
Cmp ssjm1h22 (ssjm1_h22) ;
Cmp ssjm1h32 (ssjm1_h32) ;
Cmp ssjm1h33 (ssjm1_h33) ;
ssjm1logn.set_etat_qcq() ;
ssjm1lnq.set_etat_qcq() ;
ssjm1h11.set_etat_qcq() ;
ssjm1h21.set_etat_qcq() ;
ssjm1h31.set_etat_qcq() ;
ssjm1h22.set_etat_qcq() ;
ssjm1h32.set_etat_qcq() ;
ssjm1h33.set_etat_qcq() ;
double precis_poisson = 1.e-16 ;
// Parameters for the function Scalar::poisson for logn_auto
// ---------------------------------------------------------------
Param par_logn ;
par_logn.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_logn.add_double(relax_poisson, 0) ; // relaxation parameter
par_logn.add_double(precis_poisson, 1) ; // required precision
par_logn.add_int_mod(niter, 0) ; // number of iterations actually used
par_logn.add_cmp_mod( ssjm1logn ) ;
// Parameters for the function Scalar::poisson for lnq_auto
// ---------------------------------------------------------------
Param par_lnq ;
par_lnq.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_lnq.add_double(relax_poisson, 0) ; // relaxation parameter
par_lnq.add_double(precis_poisson, 1) ; // required precision
par_lnq.add_int_mod(niter, 0) ; // number of iterations actually used -
par_lnq.add_cmp_mod( ssjm1lnq ) ;
// Parameters for the function Vector::poisson for beta method 2
// ---------------------------------------------------------------
Param par_beta2 ;
par_beta2.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_beta2.add_double(relax_poisson, 0) ; // relaxation parameter
par_beta2.add_double(precis_poisson, 1) ; // required precision
par_beta2.add_int_mod(niter, 0) ; // number of iterations actually used
Cmp ssjm1khi (ssjm1_khi) ;
Tenseur ssjm1wbeta(mp, 1, CON, mp.get_bvect_cart()) ;
ssjm1wbeta.set_etat_qcq() ;
for (int i=0; i<3; i++) {
ssjm1wbeta.set(i) = Cmp(ssjm1_wbeta(i+1)) ;
}
par_beta2.add_cmp_mod(ssjm1khi) ;
par_beta2.add_tenseur_mod(ssjm1wbeta) ;
// Parameters for the function Scalar::poisson for h11_auto
// -------------------------------------------------------------
Param par_h11 ;
par_h11.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_h11.add_double(relax_poisson, 0) ; // relaxation parameter
par_h11.add_double(precis_poisson, 1) ; // required precision
par_h11.add_int_mod(niter, 0) ; // number of iterations actually used
par_h11.add_cmp_mod( ssjm1h11 ) ;
// Parameters for the function Scalar::poisson for h21_auto
// -------------------------------------------------------------
Param par_h21 ;
par_h21.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_h21.add_double(relax_poisson, 0) ; // relaxation parameter
par_h21.add_double(precis_poisson, 1) ; // required precision
par_h21.add_int_mod(niter, 0) ; // number of iterations actually used
par_h21.add_cmp_mod( ssjm1h21 ) ;
// Parameters for the function Scalar::poisson for h31_auto
// -------------------------------------------------------------
Param par_h31 ;
par_h31.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_h31.add_double(relax_poisson, 0) ; // relaxation parameter
par_h31.add_double(precis_poisson, 1) ; // required precision
par_h31.add_int_mod(niter, 0) ; // number of iterations actually used
par_h31.add_cmp_mod( ssjm1h31 ) ;
// Parameters for the function Scalar::poisson for h22_auto
// -------------------------------------------------------------
Param par_h22 ;
par_h22.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_h22.add_double(relax_poisson, 0) ; // relaxation parameter
par_h22.add_double(precis_poisson, 1) ; // required precision
par_h22.add_int_mod(niter, 0) ; // number of iterations actually used
par_h22.add_cmp_mod( ssjm1h22 ) ;
// Parameters for the function Scalar::poisson for h32_auto
// -------------------------------------------------------------
Param par_h32 ;
par_h32.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_h32.add_double(relax_poisson, 0) ; // relaxation parameter
par_h32.add_double(precis_poisson, 1) ; // required precision
par_h32.add_int_mod(niter, 0) ; // number of iterations actually used
par_h32.add_cmp_mod( ssjm1h32 ) ;
// Parameters for the function Scalar::poisson for h33_auto
// -------------------------------------------------------------
Param par_h33 ;
par_h33.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_h33.add_double(relax_poisson, 0) ; // relaxation parameter
par_h33.add_double(precis_poisson, 1) ; // required precision
par_h33.add_int_mod(niter, 0) ; // number of iterations actually used
par_h33.add_cmp_mod( ssjm1h33 ) ;
// External potential
// See Eq (99) from Gourgoulhon et al. (2001)
// ------------------
cout << "logn_comp" << norme(logn_comp) << endl ;
cout << "pot_centri" << norme(pot_centri) << endl ;
cout << "loggam" << norme(loggam) << endl ;
Scalar pot_ext = logn_comp + pot_centri + loggam ;
Scalar ent_jm1 = ent ; // Enthalpy at previous step
Scalar source_tot(mp) ; // source term in the equation for hij_auto,
// logn_auto and beta_auto
Vector source_beta(mp, CON, mp.get_bvect_cart()) ; // source term
// in the equation for beta_auto
//=========================================================================
// Start of iteration
//=========================================================================
for(int mer=0 ; mer<mermax ; mer++ ) {
cout << "-----------------------------------------------" << endl ;
cout << "step: " << mer << endl ;
cout << "diff_ent = " << diff_ent << endl ;
//-----------------------------------------------------
// Resolution of the elliptic equation for the velocity
// scalar potential
//-----------------------------------------------------
if (irrotational) {
diff_vel_pot = velocity_potential(mermax_potvit, precis_poisson,
relax_potvit) ;
}
diff_vel_pot = 0. ; // to avoid the warning
//-----------------------------------------------------
// Computation of the new radial scale
//--------------------------------------------------
// alpha_r (r = alpha_r r') is determined so that the enthalpy
// takes the requested value ent_b at the stellar surface
// Values at the center of the star:
double logn_auto_c = logn_auto.val_grid_point(0, 0, 0, 0) ;
double pot_ext_c = pot_ext.val_grid_point(0, 0, 0, 0) ;
// Search for the reference point (theta_*, phi_*) [notation of
// Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
// at the surface of the star
// ------------------------------------------------------------
double alpha_r2 = 0 ;
for (int k=0; k<mg->get_np(l_b); k++) {
for (int j=0; j<mg->get_nt(l_b); j++) {
double pot_ext_b = pot_ext.val_grid_point(l_b, k, j, i_b) ;
double logn_auto_b = logn_auto.val_grid_point(l_b, k, j, i_b) ;
// See Eq (100) from Gourgoulhon et al. (2001)
double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b) /
( logn_auto_b - logn_auto_c ) ;
if (alpha_r2_jk > alpha_r2) {
alpha_r2 = alpha_r2_jk ;
k_b = k ;
j_b = j ;
}
}
}
alpha_r = sqrt(alpha_r2) ;
cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
<< alpha_r << endl ;
// New value of logn_auto
// ----------------------
logn_auto = alpha_r2 * logn_auto ;
logn_auto_c = logn_auto.val_grid_point(0, 0, 0, 0) ;
//------------------------------------------------------------
// Change the values of the inner points of the second domain
// by those of the outer points of the first domain
//------------------------------------------------------------
logn_auto.set_spectral_va().smooth(nzet, logn_auto.set_spectral_va()) ;
//------------------------------------------
// First integral --> enthalpy in all space
// See Eq (98) from Gourgoulhon et al. (2001)
//-------------------------------------------
ent = (ent_c + logn_auto_c + pot_ext_c) - logn_auto - pot_ext ;
cout.precision(8) ;
cout << "pot" << norme(pot_ext) << endl ;
(ent.set_spectral_va()).smooth(nzet, ent.set_spectral_va()) ;
//----------------------------------------------------
// Adaptation of the mapping to the new enthalpy field
//----------------------------------------------------
// Shall the adaptation be performed (cusp) ?
// ------------------------------------------
double dent_eq = ent.dsdr().val_point(ray_eq(),M_PI/2.,0.) ;
double dent_pole = ent.dsdr().val_point(ray_pole(),0.,0.) ;
double rap_dent = fabs( dent_eq / dent_pole ) ;
cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
if ( rap_dent < thres_adapt ) {
adapt_flag = 0 ; // No adaptation of the mapping
cout << "******* FROZEN MAPPING *********" << endl ;
}
else{
adapt_flag = 1 ; // The adaptation of the mapping is to be
// performed
}
ent_limit.set_etat_qcq() ;
for (int l=0; l<nzet; l++) { // loop on domains inside the star
ent_limit.set(l) = ent.val_grid_point(l, k_b, j_b, i_b) ;
}
ent_limit.set(nzet-1) = ent_b ;
Map_et mp_prev = mp_et ;
Cmp ent_cmp(ent) ;
mp.adapt(ent_cmp, par_adapt) ;
ent = ent_cmp ;
// Readjustment of the external boundary of domain l=nzet
// to keep a fixed ratio with respect to star's surface
if (nz>= 5) {
double separation = 2. * fabs(mp.get_ori_x()) ;
double ray_eqq = ray_eq() ;
double ray_eqq_pi = ray_eq_pi() ;
double new_rr_out_2 = (separation - ray_eqq) * 0.95 ;
double new_rr_out_3 = (separation + ray_eqq_pi) * 1.05 ;
double rr_in_1 = mp.val_r(1,-1., M_PI/2, 0.) ;
double rr_out_1 = mp.val_r(1, 1., M_PI/2, 0.) ;
double rr_out_2 = mp.val_r(2, 1., M_PI/2, 0.) ;
double rr_out_3 = mp.val_r(3, 1., M_PI/2, 0.) ;
mp.resize(1, 0.5*(new_rr_out_2 + rr_in_1) / rr_out_1) ;
mp.resize(2, new_rr_out_2 / rr_out_2) ;
mp.resize(3, new_rr_out_3 / rr_out_3) ;
for (int dd=4; dd<=nz-2; dd++) {
mp.resize(dd, new_rr_out_3 * pow(4., dd-3) /
mp.val_r(dd, 1., M_PI/2, 0.)) ;
}
}
else {
cout << "too small number of domains" << endl ;
}
//----------------------------------------------------
// Computation of the enthalpy at the new grid points
//----------------------------------------------------
mp_prev.homothetie(alpha_r) ;
Cmp ent_cmp2 (ent) ;
mp.reevaluate_symy(&mp_prev, nzet+1, ent_cmp2) ;
ent = ent_cmp2 ;
double ent_s_max = -1 ;
int k_s_max = -1 ;
int j_s_max = -1 ;
for (int k=0; k<mg->get_np(l_b); k++) {
for (int j=0; j<mg->get_nt(l_b); j++) {
double xx = fabs( ent.val_grid_point(l_b, k, j, i_b) ) ;
if (xx > ent_s_max) {
ent_s_max = xx ;
k_s_max = k ;
j_s_max = j ;
}
}
}
cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
<< " and nzet : " << endl ;
cout << " " << ent_s_max << " reached for k = " << k_s_max <<
" and j = " << j_s_max << endl ;
//----------------------------------------------------
// Equation of state
//----------------------------------------------------
equation_of_state() ; // computes new values for nbar (n), ener (e)
// and press (p) from the new ent (H)
//---------------------------------------------------------
// Matter source terms in the gravitational field equations
//---------------------------------------------------------
hydro_euler() ; // computes new values for ener_euler (E),
// s_euler (S) and u_euler (U^i)
// -------------------------------
// AUXILIARY QUANTITIES
// -------------------------------
// Derivatives of N and logN
//--------------------------
const Vector dcov_logn_auto = logn_auto.derive_cov(flat) ;
Tensor dcovdcov_logn_auto = (logn_auto.derive_cov(flat))
.derive_cov(flat) ;
dcovdcov_logn_auto.inc_dzpuis() ;
// Derivatives of lnq, phi and Q
//-------------------------------
const Scalar phi (0.5 * (lnq - logn)) ;
const Scalar phi_auto (0.5 * (lnq_auto - logn_auto)) ;
const Vector dcov_phi_auto = phi_auto.derive_cov(flat) ;
const Vector dcov_lnq = 2*dcov_phi + dcov_logn ;
const Vector dcon_lnq = 2*dcon_phi + dcon_logn ;
const Vector dcov_lnq_auto = lnq_auto.derive_cov(flat) ;
Tensor dcovdcov_lnq_auto = dcov_lnq_auto.derive_cov(flat) ;
dcovdcov_lnq_auto.inc_dzpuis() ;
Scalar qq = exp(lnq) ;
qq.std_spectral_base() ;
const Vector& dcov_qq = qq.derive_cov(flat) ;
const Scalar& divbeta_auto = beta_auto.divergence(flat) ;
const Tensor& dcov_beta_auto = beta_auto.derive_cov(flat) ;
Tensor dcovdcov_beta_auto = beta_auto.derive_cov(flat)
.derive_cov(flat) ;
dcovdcov_beta_auto.inc_dzpuis() ;
// Derivatives of hij, gtilde...
//------------------------------
Scalar psi2 (pow(psi4, 0.5)) ;
psi2.std_spectral_base() ;
const Tensor& dcov_hij = hij.derive_cov(flat) ;
const Tensor& dcon_hij = hij.derive_con(flat) ;
const Tensor& dcov_hij_auto = hij_auto.derive_cov(flat) ;
const Sym_tensor& gtilde_cov = gtilde.cov() ;
const Sym_tensor& gtilde_con = gtilde.con() ;
const Tensor& dcov_gtilde = gtilde_cov.derive_cov(flat) ;
Connection gamijk (gtilde, flat) ;
const Tensor& deltaijk = gamijk.get_delta() ;
// H^i and its derivatives ( = O in Dirac gauge)
// ---------------------------------------------
double lambda_dirac = 0. ;
const Vector hdirac = lambda_dirac * hij.divergence(flat) ;
const Vector hdirac_auto = lambda_dirac * hij_auto.divergence(flat) ;
Tensor dcov_hdirac = lambda_dirac * hdirac.derive_cov(flat) ;
dcov_hdirac.inc_dzpuis() ;
Tensor dcov_hdirac_auto = lambda_dirac * hdirac_auto.derive_cov(flat) ;
dcov_hdirac_auto.inc_dzpuis() ;
Tensor dcon_hdirac_auto = lambda_dirac * hdirac_auto.derive_con(flat) ;
dcon_hdirac_auto.inc_dzpuis() ;
//--------------------------------------------------------
// Poisson equation for logn_auto (nu_auto)
//--------------------------------------------------------
// Source
//--------
int nr = mp.get_mg()->get_nr(0) ;
int nt = mp.get_mg()->get_nt(0) ;
int np = mp.get_mg()->get_np(0) ;
Scalar source1(mp) ;
Scalar source2(mp) ;
Scalar source3(mp) ;
Scalar source4(mp) ;
Scalar source5(mp) ;
Scalar source6(mp) ;
Scalar source7(mp) ;
Scalar source8(mp) ;
source1 = qpig * psi4 % (ener_euler + s_euler) ;
source2 = psi4 % (kcar_auto + kcar_comp) ;
source3 = - contract(dcov_logn_auto, 0, dcon_logn, 0, true)
- 2. * contract(contract(gtilde_con, 0, dcov_phi, 0),
0, dcov_logn_auto, 0, true) ;
source4 = - contract(hij, 0, 1, dcovdcov_logn_auto +
dcov_logn_auto*dcov_logn, 0, 1) ;
source5 = - contract(hdirac, 0, dcov_logn_auto, 0) ;
source_tot = source1 + source2 + source3 + source4 + source5 ;
cout << "moyenne de la source 1 pour logn_auto" << endl ;
cout << norme(source1/(nr*nt*np)) << endl ;
cout << "moyenne de la source 2 pour logn_auto" << endl ;
cout << norme(source2/(nr*nt*np)) << endl ;
cout << "moyenne de la source 3 pour logn_auto" << endl ;
cout << norme(source3/(nr*nt*np)) << endl ;
cout << "moyenne de la source 4 pour logn_auto" << endl ;
cout << norme(source4/(nr*nt*np)) << endl ;
cout << "moyenne de la source 5 pour logn_auto" << endl ;
cout << norme(source5/(nr*nt*np)) << endl ;
cout << "moyenne de la source pour logn_auto" << endl ;
cout << norme(source_tot/(nr*nt*np)) << endl ;
// Resolution of the Poisson equation
// ----------------------------------
source_tot.poisson(par_logn, logn_auto) ;
ssjm1_logn = ssjm1logn ;
cout << "logn_auto" << endl << norme(logn_auto/(nr*nt*np)) << endl ;
// Check: has the Poisson equation been correctly solved ?
// -----------------------------------------------------
Tbl tdiff_logn = diffrel(logn_auto.laplacian(), source_tot) ;
cout <<
"Relative error in the resolution of the equation for logn_auto : "
<< endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_logn(l) << " " ;
}
cout << endl ;
diff_logn = max(abs(tdiff_logn)) ;
//--------------------------------------------------------
// Poisson equation for lnq_auto
//--------------------------------------------------------
// Source
//--------
source1 = qpig * psi4 % s_euler ;
source2 = 0.75 * psi4 % (kcar_auto + kcar_comp) ;
source3 = - contract(dcon_lnq, 0, dcov_lnq_auto, 0, true) ;
source4 = 2. * contract(contract(gtilde_con, 0, dcov_phi, 0), 0,
dcov_phi_auto + dcov_logn_auto, 0, true) ;
source5 = 0.0625 * contract(gtilde_con, 0, 1, contract(
dcov_hij_auto, 0, 1, dcov_gtilde, 0, 1), 0, 1) ;
source6 = - 0.125 * contract(gtilde_con, 0, 1, contract(
dcov_hij_auto, 0, 1, dcov_gtilde, 0, 2), 0, 1) ;
source7 = - contract(hij, 0, 1, dcovdcov_lnq_auto + dcov_lnq_auto *
dcov_lnq, 0, 1) ;
source8 = - 0.25 * contract(dcov_hdirac_auto, 0, 1)
- contract(hdirac, 0, dcov_lnq_auto, 0) ;
source_tot = source1 + source2 + source3 + source4 + source5 +
source6 + source7 + source8 ;
cout << "moyenne de la source 1 pour lnq_auto" << endl ;
cout << norme(source1/(nr*nt*np)) << endl ;
cout << "moyenne de la source 2 pour lnq_auto" << endl ;
cout << norme(source2/(nr*nt*np)) << endl ;
cout << "moyenne de la source 3 pour lnq_auto" << endl ;
cout << norme(source3/(nr*nt*np)) << endl ;
cout << "moyenne de la source 4 pour lnq_auto" << endl ;
cout << norme(source4/(nr*nt*np)) << endl ;
cout << "moyenne de la source 5 pour lnq_auto" << endl ;
cout << norme(source5/(nr*nt*np)) << endl ;
cout << "moyenne de la source 6 pour lnq_auto" << endl ;
cout << norme(source6/(nr*nt*np)) << endl ;
cout << "moyenne de la source 7 pour lnq_auto" << endl ;
cout << norme(source7/(nr*nt*np)) << endl ;
cout << "moyenne de la source 8 pour lnq_auto" << endl ;
cout << norme(source8/(nr*nt*np)) << endl ;
cout << "moyenne de la source pour lnq_auto" << endl ;
cout << norme(source_tot/(nr*nt*np)) << endl ;
// Resolution of the Poisson equation
// ----------------------------------
source_tot.poisson(par_lnq, lnq_auto) ;
ssjm1_lnq = ssjm1lnq ;
cout << "lnq_auto" << endl << norme(lnq_auto/(nr*nt*np)) << endl ;
// Check: has the Poisson equation been correctly solved
// -----------------------------------------------------
Tbl tdiff_lnq = diffrel(lnq_auto.laplacian(), source_tot) ;
cout <<
"Relative error in the resolution of the equation for lnq : "
<< endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_lnq(l) << " " ;
}
cout << endl ;
diff_lnq = max(abs(tdiff_lnq)) ;
//--------------------------------------------------------
// Vector Poisson equation for beta_auto
//--------------------------------------------------------
// Source
//--------
Vector source1_beta(mp, CON, mp.get_bvect_cart()) ;
Vector source2_beta(mp, CON, mp.get_bvect_cart()) ;
Vector source3_beta(mp, CON, mp.get_bvect_cart()) ;
Vector source4_beta(mp, CON, mp.get_bvect_cart()) ;
Vector source5_beta(mp, CON, mp.get_bvect_cart()) ;
Vector source6_beta(mp, CON, mp.get_bvect_cart()) ;
Vector source7_beta(mp, CON, mp.get_bvect_cart()) ;
source1_beta = (4.*qpig) * nn % psi4
%(ener_euler + press) * u_euler ;
source2_beta = 2. * nn * contract(tkij_auto, 1,
dcov_logn - 6 * dcov_phi, 0) ;
source3_beta = - 2. * nn * contract(tkij_auto, 0, 1, deltaijk,
1, 2) ;
source4_beta = - contract(hij, 0, 1, dcovdcov_beta_auto, 1, 2) ;
source5_beta = - 0.3333333333333333 * contract(hij, 1, contract(
dcovdcov_beta_auto, 0, 1), 0) ;
source6_beta.set_etat_zero() ; //hdirac_auto.derive_lie(omdsdp) ;
source6_beta.inc_dzpuis() ;
source7_beta = contract(beta, 0, dcov_hdirac_auto, 1) ;
+ 2./3. * hdirac * divbeta_auto
- contract(hdirac, 0, dcov_beta_auto, 1) ;
source_beta = source1_beta + source2_beta + source3_beta
+ source4_beta + source5_beta + source6_beta + source7_beta ;
cout << "moyenne de la source 1 pour beta_auto" << endl ;
cout << norme(source1_beta(1)/(nr*nt*np)) << endl ;
cout << norme(source1_beta(2)/(nr*nt*np)) << endl ;
cout << norme(source1_beta(3)/(nr*nt*np)) << endl ;
cout << "moyenne de la source 2 pour beta_auto" << endl ;
cout << norme(source2_beta(1)/(nr*nt*np)) << endl ;
cout << norme(source2_beta(2)/(nr*nt*np)) << endl ;
cout << norme(source2_beta(3)/(nr*nt*np)) << endl ;
cout << "moyenne de la source 3 pour beta_auto" << endl ;
cout << norme(source3_beta(1)/(nr*nt*np)) << endl ;
cout << norme(source3_beta(2)/(nr*nt*np)) << endl ;
cout << norme(source3_beta(3)/(nr*nt*np)) << endl ;
cout << "moyenne de la source 4 pour beta_auto" << endl ;
cout << norme(source4_beta(1)/(nr*nt*np)) << endl ;
cout << norme(source4_beta(2)/(nr*nt*np)) << endl ;
cout << norme(source4_beta(3)/(nr*nt*np)) << endl ;
cout << "moyenne de la source 5 pour beta_auto" << endl ;
cout << norme(source5_beta(1)/(nr*nt*np)) << endl ;
cout << norme(source5_beta(2)/(nr*nt*np)) << endl ;
cout << norme(source5_beta(3)/(nr*nt*np)) << endl ;
cout << "moyenne de la source 6 pour beta_auto" << endl ;
cout << norme(source6_beta(1)/(nr*nt*np)) << endl ;
cout << norme(source6_beta(2)/(nr*nt*np)) << endl ;
cout << norme(source6_beta(3)/(nr*nt*np)) << endl ;
cout << "moyenne de la source 7 pour beta_auto" << endl ;
cout << norme(source7_beta(1)/(nr*nt*np)) << endl ;
cout << norme(source7_beta(2)/(nr*nt*np)) << endl ;
cout << norme(source7_beta(3)/(nr*nt*np)) << endl ;
cout << "moyenne de la source pour beta_auto" << endl ;
cout << norme(source_beta(1)/(nr*nt*np)) << endl ;
cout << norme(source_beta(2)/(nr*nt*np)) << endl ;
cout << norme(source_beta(3)/(nr*nt*np)) << endl ;
// Resolution of the Poisson equation
// ----------------------------------
// Filter for the source of beta vector
for (int i=1; i<=3; i++) {
if (source_beta(i).get_etat() != ETATZERO)
source_beta.set(i).filtre(4) ;
}
for (int i=1; i<=3; i++) {
if(source_beta(i).dz_nonzero()) {
assert( source_beta(i).get_dzpuis() == 4 ) ;
}
else{
(source_beta.set(i)).set_dzpuis(4) ;
}
}
double lambda = double(1) / double(3) ;
Tenseur source_p(mp, 1, CON, mp.get_bvect_cart() ) ;
source_p.set_etat_qcq() ;
for (int i=0; i<3; i++) {
source_p.set(i) = Cmp(source_beta(i+1)) ;
}
Tenseur vect_auxi (mp, 1, CON, mp.get_bvect_cart()) ;
vect_auxi.set_etat_qcq() ;
for (int i=0; i<3 ;i++){
vect_auxi.set(i) = 0. ;
}
Tenseur scal_auxi (mp) ;
scal_auxi.set_etat_qcq() ;
scal_auxi.set().annule_hard() ;
scal_auxi.set_std_base() ;
Tenseur resu_p(mp, 1, CON, mp.get_bvect_cart() ) ;
resu_p.set_etat_qcq() ;
for (int i=0; i<3 ;i++)
resu_p.set(i).annule_hard() ;
resu_p.set_std_base() ;
//source_p.poisson_vect(lambda, par_beta2, resu_p, vect_auxi,
// scal_auxi) ;
source_p.poisson_vect_oohara(lambda, par_beta2, resu_p, scal_auxi) ;
for (int i=1; i<=3; i++)
beta_auto.set(i) = resu_p(i-1) ;
ssjm1_khi = ssjm1khi ;
for (int i=0; i<3; i++){
ssjm1_wbeta.set(i+1) = ssjm1wbeta(i) ;
}
cout << "beta_auto_x" << endl << norme(beta_auto(1)/(nr*nt*np))
<< endl ;
cout << "beta_auto_y" << endl << norme(beta_auto(2)/(nr*nt*np))
<< endl ;
cout << "beta_auto_z" << endl << norme(beta_auto(3)/(nr*nt*np))
<< endl ;
// Check: has the equation for beta_auto been correctly solved ?
// --------------------------------------------------------------
Vector lap_beta = (beta_auto.derive_con(flat)).divergence(flat)
+ lambda* beta_auto.divergence(flat).derive_con(flat) ;
source_beta.dec_dzpuis() ;
Tbl tdiff_beta_x = diffrel(lap_beta(1), source_beta(1)) ;
Tbl tdiff_beta_y = diffrel(lap_beta(2), source_beta(2)) ;
Tbl tdiff_beta_z = diffrel(lap_beta(3), source_beta(3)) ;
cout <<
"Relative error in the resolution of the equation for beta_auto : "
<< endl ;
cout << "x component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_beta_x(l) << " " ;
}
cout << endl ;
cout << "y component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_beta_y(l) << " " ;
}
cout << endl ;
cout << "z component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_beta_z(l) << " " ;
}
cout << endl ;
diff_beta_x = max(abs(tdiff_beta_x)) ;
diff_beta_y = max(abs(tdiff_beta_y)) ;
diff_beta_z = max(abs(tdiff_beta_z)) ;
if (!conf_flat) {
//--------------------------------------------------------
// Poisson equation for hij
//--------------------------------------------------------
// Declaration of all sources
//---------------------------
Scalar source_tot_hij(mp) ;
Tensor source_Sij(mp, 2, CON, mp.get_bvect_cart()) ;
Tensor source_Rij(mp, 2, CON, mp.get_bvect_cart()) ;
Tensor tens_temp(mp, 2, CON, mp.get_bvect_cart()) ;
Tensor source_1 (mp, 2, CON, mp.get_bvect_cart()) ;
Tensor source_2 (mp, 2, CON, mp.get_bvect_cart()) ;
Tensor source_3a (mp, 2, CON, mp.get_bvect_cart()) ;
Tensor source_3b (mp, 2, CON, mp.get_bvect_cart()) ;
Tensor source_4 (mp, 2, CON, mp.get_bvect_cart()) ;
Tensor source_5 (mp, 2, CON, mp.get_bvect_cart()) ;
Tensor source_6 (mp, 2, CON, mp.get_bvect_cart()) ;
// Sources
//-----------
source_1 = contract(dcon_hij, 1, dcov_lnq_auto, 0) ;
source_2 = - contract(dcon_hij, 2, dcov_lnq_auto, 0)
- 2./3. * contract(hdirac, 0, dcov_lnq_auto, 0) * flat.con() ;
// Lie derivative of A^{ij}
// --------------------------
Scalar decouple_logn = (logn_auto - 1.e-8)/(logn - 2.e-8) ;
// Function exp(-(r-r_0)^2/sigma^2)
// --------------------------------
double r0 = mp.val_r(nz-2, 1, 0, 0) ;
double sigma = 1.*r0 ;
Scalar rr (mp) ;
rr = mp.r ;
Scalar ff (mp) ;
ff = exp( -(rr - r0)*(rr - r0)/sigma/sigma ) ;
for (int ii=0; ii<nz-1; ii++)
ff.set_domain(ii) = 1. ;
ff.set_outer_boundary(nz-1, 0) ;
ff.std_spectral_base() ;
// ff.annule_domain(nz-1) ;
//des_profile(ff, 0, 20, 0, 0) ;
// Construction of Omega d/dphi
// ----------------------------
// Construction of D_k \Phi^i
Itbl type (2) ;
type.set(0) = CON ;
type.set(1) = COV ;
Tensor dcov_omdsdphi (mp, 2, type, mp.get_bvect_cart()) ;
dcov_omdsdphi.set(1,1) = 0. ;
dcov_omdsdphi.set(2,1) = om * ff ;
dcov_omdsdphi.set(3,1) = 0. ;
dcov_omdsdphi.set(2,2) = 0. ;
dcov_omdsdphi.set(3,2) = 0. ;
dcov_omdsdphi.set(3,3) = 0. ;
dcov_omdsdphi.set(1,2) = -om * ff ;
dcov_omdsdphi.set(1,3) = 0. ;
dcov_omdsdphi.set(2,3) = 0. ;
dcov_omdsdphi.std_spectral_base() ;
source_3a = contract(tkij_auto, 0, dcov_omdsdphi, 1) ;
source_3a.inc_dzpuis() ;
// Source 3b
// ------------
Vector omdsdp (mp, CON, mp.get_bvect_cart()) ;
Scalar yya (mp) ;
yya = mp.ya ;
Scalar xxa (mp) ;
xxa = mp.xa ;
Scalar zza (mp) ;
zza = mp.za ;
if (fabs(mp.get_rot_phi()) < 1e-10){
omdsdp.set(1) = - om * yya * ff ;
omdsdp.set(2) = om * xxa * ff ;
omdsdp.set(3).annule_hard() ;
}
else{
omdsdp.set(1) = om * yya * ff ;
omdsdp.set(2) = - om * xxa * ff ;
omdsdp.set(3).annule_hard() ;
}
omdsdp.set(1).set_outer_boundary(nz-1, 0) ;
omdsdp.set(2).set_outer_boundary(nz-1, 0) ;
omdsdp.std_spectral_base() ;
source_3b = - contract(omdsdp, 0, tkij_auto.derive_cov(flat), 2) ;
// Source 4
// ---------
source_4 = - tkij_auto.derive_lie(beta) ;
source_4.inc_dzpuis() ;
source_4 += - 2./3. * beta.divergence(flat) * tkij_auto ;
source_5 = dcon_hdirac_auto ;
source_6 = - 2./3. * hdirac_auto.divergence(flat) * flat.con() ;
source_6.inc_dzpuis() ;
// Source terms for Sij
//---------------------
source_Sij = 8. * nn / psi4 * phi_auto.derive_con(gtilde) *
contract(gtilde_con, 0, dcov_phi, 0) ;
source_Sij += 4. / psi4 * phi_auto.derive_con(gtilde) *
nn * contract(gtilde_con, 0, dcov_logn, 0) +
4. / psi4 * nn * contract(gtilde_con, 0, dcov_logn, 0) *
phi_auto.derive_con(gtilde) ;
source_Sij += - nn / (3.*psi4) * gtilde_con *
( 0.25 * contract(gtilde_con, 0, 1, contract(dcov_hij_auto, 0, 1,
dcov_gtilde, 0, 1), 0, 1)
- 0.5 * contract(gtilde_con, 0, 1, contract(dcov_hij_auto, 0, 1,
dcov_gtilde, 0, 2), 0, 1)) ;
source_Sij += - 8.*nn / (3.*psi4) * gtilde_con *
contract(dcov_phi_auto, 0, contract(gtilde_con, 0, dcov_phi, 0), 0) ;
tens_temp = nn / (3.*psi4) * hdirac.divergence(flat)*hij_auto ;
tens_temp.inc_dzpuis() ;
source_Sij += tens_temp ;
source_Sij += - 8./(3.*psi4) * contract(dcov_phi_auto, 0,
nn*contract(gtilde_con, 0, dcov_logn, 0), 0) * gtilde_con ;
source_Sij += 2.*nn* contract(gtilde_cov, 0, 1, tkij_auto *
(tkij_auto+tkij_comp), 1, 3) ;
source_Sij += - 2. * qpig * nn * ( psi4 * stress_euler
- 0.33333333333333333 * s_euler * gtilde_con ) ;
source_Sij += - 1./(psi4*psi2) * contract(gtilde_con, 1,
contract(gtilde_con, 1, qq*dcovdcov_lnq_auto +
qq*dcov_lnq_auto*dcov_lnq, 0), 1) ;
source_Sij += - 0.5/(psi4*psi2) * contract(contract(hij, 1,
dcov_hij_auto, 2), 1, dcov_qq, 0) -
0.5/(psi4*psi2) * contract(contract(dcov_hij_auto, 2,
hij, 1), 1, dcov_qq, 0) ;
source_Sij += 0.5/(psi4*psi2) * contract(contract(hij, 0,
dcov_hij_auto, 2), 0, dcov_qq, 0) ;
source_Sij += 1./(3.*psi4*psi2)*contract(gtilde_con, 0, 1,
qq*dcovdcov_lnq_auto + qq*dcov_lnq_auto*dcov_lnq, 0, 1)
*gtilde_con ;
source_Sij += 1./(3.*psi4*psi2) * contract(hdirac, 0,
dcov_qq, 0)*hij_auto ;
// Source terms for Rij
//---------------------
source_Rij = contract(hij, 0, 1, dcov_hij_auto.derive_cov(flat),
2, 3) ;
source_Rij.inc_dzpuis() ;
source_Rij += - contract(hij_auto, 1, dcov_hdirac, 1) -
contract(dcov_hdirac, 1, hij_auto, 1) ;
source_Rij += contract(hdirac, 0, dcov_hij_auto, 2) ;
source_Rij += - contract(contract(dcov_hij_auto, 1, dcov_hij, 2),
1, 3) ;
source_Rij += - contract(gtilde_cov, 0, 1, contract(contract(
gtilde_con, 0, dcov_hij_auto, 2), 0, dcov_hij, 2), 1, 3) ;
source_Rij += contract(contract(contract(contract(gtilde_cov, 0,
dcov_hij_auto, 1), 2, gtilde_con, 1), 0, dcov_hij, 1), 0, 3) +
contract(contract(contract(contract(gtilde_cov, 0,
dcov_hij_auto, 1), 0, dcov_hij, 1), 0, 3), 0, gtilde_con, 1) ;
source_Rij += 0.5 * contract(gtilde_con*gtilde_con, 1, 3,
contract(dcov_hij_auto, 0, 1, dcov_gtilde, 0, 1), 0, 1) ;
source_Rij = source_Rij * 0.5 ;
for(int i=1; i<=3; i++)
for(int j=1; j<=i; j++) {
source_tot_hij = source_1(i,j) + source_1(j,i)
+ source_2(i,j) + 2.*psi4/nn * (
source_4(i,j) - source_Sij(i,j))
- 2.* source_Rij(i,j) +
source_5(i,j) + source_5(j,i) + source_6(i,j) ;
source_tot_hij.dec_dzpuis() ;
source3 = 2.*psi4/nn * (source_3a(i,j) + source_3a(j,i) +
source_3b(i,j)) ;
source_tot_hij = source_tot_hij + source3 ;
//source_tot_hij.inc_dzpuis() ;
cout << "source_mat" << endl
<< norme((- 2. * qpig * nn * ( psi4 * stress_euler
- 0.33333333333333333 * s_euler * gtilde_con ))
(i,j))/(nr*nt*np) << endl ;
cout << "max source_mat" << endl
<< max((- 2. * qpig * nn * ( psi4 * stress_euler
- 0.33333333333333333 * s_euler * gtilde_con ))
(i,j)) << endl ;
cout << "source1" << endl
<< norme(source_1(i,j)/(nr*nt*np)) << endl ;
cout << "max source1" << endl
<< max(source_1(i,j)) << endl ;
cout << "source2" << endl
<< norme(source_2(i,j)/(nr*nt*np)) << endl ;
cout << "max source2" << endl
<< max(source_2(i,j)) << endl ;
cout << "source3a" << endl
<< norme(source_3a(i,j)/(nr*nt*np)) << endl ;
cout << "max source3a" << endl
<< max(source_3a(i,j)) << endl ;
cout << "source3b" << endl
<< norme(source_3b(i,j)/(nr*nt*np)) << endl ;
cout << "max source3b" << endl
<< max(source_3b(i,j)) << endl ;
cout << "source4" << endl
<< norme(source_4(i,j)/(nr*nt*np)) << endl ;
cout << "max source4" << endl
<< max(source_4(i,j)) << endl ;
cout << "source5" << endl
<< norme(source_5(i,j)/(nr*nt*np)) << endl ;
cout << "max source5" << endl
<< max(source_5(i,j)) << endl ;
cout << "source6" << endl
<< norme(source_6(i,j)/(nr*nt*np)) << endl ;
cout << "max source6" << endl
<< max(source_6(i,j)) << endl ;
cout << "source_Rij" << endl
<< norme(source_Rij(i,j)/(nr*nt*np)) << endl ;
cout << "max source_Rij" << endl
<< max(source_Rij(i,j)) << endl ;
cout << "source_Sij" << endl
<< norme(source_Sij(i,j)/(nr*nt*np)) << endl ;
cout << "max source_Sij" << endl
<< max(source_Sij(i,j)) << endl ;
cout << "source_tot" << endl
<< norme(source_tot_hij/(nr*nt*np)) << endl ;
cout << "max source_tot" << endl
<< max(source_tot_hij) << endl ;
// Resolution of the Poisson equations and
// Check: has the Poisson equation been correctly solved ?
// -----------------------------------------------------
if(i==1 && j==1) {
source_tot_hij.poisson(par_h11, hij_auto.set(1,1)) ;
Scalar laplac (hij_auto(1,1).laplacian()) ;
laplac.dec_dzpuis() ;
Tbl tdiff_h11 = diffrel(laplac, source_tot_hij) ;
cout << "Relative error in the resolution of the equation for "
<< "h11_auto : " << endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_h11(l) << " " ;
}
cout << endl ;
diff_h11 = max(abs(tdiff_h11)) ;
}
if(i==2 && j==1) {
source_tot_hij.poisson(par_h21, hij_auto.set(2,1)) ;
Scalar laplac (hij_auto(2,1).laplacian()) ;
laplac.dec_dzpuis() ;
Tbl tdiff_h21 = diffrel(laplac, source_tot_hij) ;
cout <<
"Relative error in the resolution of the equation for "
<< "h21_auto : " << endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_h21(l) << " " ;
}
cout << endl ;
diff_h21 = max(abs(tdiff_h21)) ;
}
if(i==3 && j==1) {
source_tot_hij.poisson(par_h31, hij_auto.set(3,1)) ;
Scalar laplac (hij_auto(3,1).laplacian()) ;
laplac.dec_dzpuis() ;
Tbl tdiff_h31 = diffrel(laplac, source_tot_hij) ;
cout <<
"Relative error in the resolution of the equation for "
<< "h31_auto : " << endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_h31(l) << " " ;
}
cout << endl ;
diff_h31 = max(abs(tdiff_h31)) ;
}
if(i==2 && j==2) {
source_tot_hij.poisson(par_h22, hij_auto.set(2,2)) ;
Scalar laplac (hij_auto(2,2).laplacian()) ;
laplac.dec_dzpuis() ;
Tbl tdiff_h22 = diffrel(laplac, source_tot_hij) ;
cout <<
"Relative error in the resolution of the equation for "
<< "h22_auto : " << endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_h22(l) << " " ;
}
cout << endl ;
diff_h22 = max(abs(tdiff_h22)) ;
}
if(i==3 && j==2) {
source_tot_hij.poisson(par_h32, hij_auto.set(3,2)) ;
Scalar laplac (hij_auto(3,2).laplacian()) ;
laplac.dec_dzpuis() ;
Tbl tdiff_h32 = diffrel(laplac, source_tot_hij) ;
cout <<
"Relative error in the resolution of the equation for "
<< "h32_auto : " << endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_h32(l) << " " ;
}
cout << endl ;
diff_h32 = max(abs(tdiff_h32)) ;
}
if(i==3 && j==3) {
source_tot_hij.poisson(par_h33, hij_auto.set(3,3)) ;
Scalar laplac (hij_auto(3,3).laplacian()) ;
laplac.dec_dzpuis() ;
Tbl tdiff_h33 = diffrel(laplac, source_tot_hij) ;
cout <<
"Relative error in the resolution of the equation for "
<< "h33_auto : " << endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_h33(l) << " " ;
}
cout << endl ;
diff_h33 = max(abs(tdiff_h33)) ;
}
}
cout << "Tenseur hij auto in cartesian coordinates" << endl ;
for (int i=1; i<=3; i++)
for (int j=1; j<=i; j++) {
cout << " Comp. " << i << " " << j << " : " ;
for (int l=0; l<nz; l++){
cout << norme(hij_auto(i,j)/(nr*nt*np))(l) << " " ;
}
cout << endl ;
}
cout << endl ;
ssjm1_h11 = ssjm1h11 ;
ssjm1_h21 = ssjm1h21 ;
ssjm1_h31 = ssjm1h31 ;
ssjm1_h22 = ssjm1h22 ;
ssjm1_h32 = ssjm1h32 ;
ssjm1_h33 = ssjm1h33 ;
}
// End of relativistic equations
//-------------------------------------------------
// Relative change in enthalpy
//-------------------------------------------------
Tbl diff_ent_tbl = diffrel( ent, ent_jm1 ) ;
diff_ent = diff_ent_tbl(0) ;
for (int l=1; l<nzet; l++) {
diff_ent += diff_ent_tbl(l) ;
}
diff_ent /= nzet ;
ent_jm1 = ent ;
} // End of main loop
//=========================================================================
// End of iteration
//=========================================================================
}
}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition cmp.C:304
Lorene::Cmp::annule_hard
void annule_hard()
Sets the Cmp to zero in a hard way.
Definition cmp.C:338
Lorene::Connection
Class Connection.
Definition connection.h:113
Lorene::Connection::get_delta
const Tensor_sym & get_delta() const
Returns the tensor which defines the connection with respect to the flat one: is the difference bet...
Definition connection.h:271
Lorene::Itbl
Basic integer array class.
Definition itbl.h:122
Lorene::Itbl::set
int & set(int i)
Read/write of a particular element (index i ) (1D case).
Definition itbl.h:247
Lorene::Map_et
Radial mapping of rather general form.
Definition map.h:2752
Lorene::Map_et::homothetie
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition map_et.C:905
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::add_double
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition param.C:315
Lorene::Param::add_cmp_mod
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Definition param.C:1004
Lorene::Param::add_int_mod
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition param.C:385
Lorene::Param::add_tenseur_mod
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.
Definition param.C:1142
Lorene::Param::add_int
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
Lorene::Param::add_tbl
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
Definition param.C:522
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::filtre
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition scalar_manip.C:151
Lorene::Scalar::poisson
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition scalar_pde.C:136
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::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
Lorene::Scalar::set_outer_boundary
void set_outer_boundary(int l, double x)
Sets the value of the Scalar at the outer boundary of a given domain.
Definition scalar_manip.C:315
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::Star_bin::loggam
Scalar loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition star.h:512
Lorene::Star_bin::lnq_auto
Scalar lnq_auto
Scalar field generated principally by the star.
Definition star.h:543
Lorene::Star_bin::ssjm1_khi
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for khi.
Definition star.h:634
Lorene::Star_bin::dcon_logn
Vector dcon_logn
Contravariant derivative of the total logarithm of the lapse.
Definition star.h:538
Lorene::Star_bin::dcon_phi
Vector dcon_phi
Contravariant derivative of the logarithm of the conformal factor.
Definition star.h:557
Lorene::Star_bin::velocity_potential
double velocity_potential(int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Definition star_bin_vel_pot.C:76
Lorene::Star_bin::ssjm1_h11
Scalar ssjm1_h11
Effective source at the previous step for the resolution of the Poisson equation for h00_auto.
Definition star.h:641
Lorene::Star_bin::ssjm1_logn
Scalar ssjm1_logn
Effective source at the previous step for the resolution of the Poisson equation for logn_auto.
Definition star.h:623
Lorene::Star_bin::dcov_logn
Vector dcov_logn
Covariant derivative of the total logarithm of the lapse.
Definition star.h:535
Lorene::Star_bin::logn_auto
Scalar logn_auto
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principally by...
Definition star.h:527
Lorene::Star_bin::irrotational
bool irrotational
true for an irrotational star, false for a corotating one
Definition star.h:491
Lorene::Star_bin::kcar_comp
Scalar kcar_comp
Part of the scalar generated by beta_auto and beta_comp, i.e.
Definition star.h:618
Lorene::Star_bin::logn_comp
Scalar logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principally by...
Definition star.h:532
Lorene::Star_bin::ssjm1_h21
Scalar ssjm1_h21
Effective source at the previous step for the resolution of the Poisson equation for h10_auto.
Definition star.h:646
Lorene::Star_bin::hydro_euler
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Definition star_bin_hydro.C:66
Lorene::Star_bin::tkij_comp
Sym_tensor tkij_comp
Part of the extrinsic curvature tensor generated by beta_comp.
Definition star.h:606
Lorene::Star_bin::ssjm1_h33
Scalar ssjm1_h33
Effective source at the previous step for the resolution of the Poisson equation for h22_auto.
Definition star.h:666
Lorene::Star_bin::kcar_auto
Scalar kcar_auto
Part of the scalar generated by beta_auto, i.e.
Definition star.h:612
Lorene::Star_bin::beta_auto
Vector beta_auto
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Definition star.h:570
Lorene::Star_bin::hij
Sym_tensor hij
Total deviation of the inverse conformal metric from the inverse flat metric.
Definition star.h:581
Lorene::Star_bin::hij_auto
Sym_tensor hij_auto
Deviation of the inverse conformal metric from the inverse flat metric generated principally by the ...
Definition star.h:588
Lorene::Star_bin::psi4
Scalar psi4
Conformal factor .
Definition star.h:552
Lorene::Star_bin::ssjm1_h31
Scalar ssjm1_h31
Effective source at the previous step for the resolution of the Poisson equation for h20_auto.
Definition star.h:651
Lorene::Star_bin::tkij_auto
Sym_tensor tkij_auto
Part of the extrinsic curvature tensor generated by beta_auto.
Definition star.h:600
Lorene::Star_bin::dcov_phi
Vector dcov_phi
Covariant derivative of the logarithm of the conformal factor.
Definition star.h:555
Lorene::Star_bin::flat
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) .
Definition star.h:562
Lorene::Star_bin::ssjm1_h22
Scalar ssjm1_h22
Effective source at the previous step for the resolution of the Poisson equation for h11_auto.
Definition star.h:656
Lorene::Star_bin::ssjm1_h32
Scalar ssjm1_h32
Effective source at the previous step for the resolution of the Poisson equation for h21_auto.
Definition star.h:661
Lorene::Star_bin::equilibrium
void equilibrium(double ent_c, int mermax, int mermax_potvit, int mermax_poisson, double relax_poisson, double relax_potvit, double thres_adapt, Tbl &diff, double om)
Computes an equilibrium configuration.
Definition star_bin_equilibrium.C:143
Lorene::Star_bin::gtilde
Metric gtilde
Conformal metric .
Definition star.h:565
Lorene::Star_bin::pot_centri
Scalar pot_centri
Centrifugal potential.
Definition star.h:521
Lorene::Star_bin::ssjm1_lnq
Scalar ssjm1_lnq
Effective source at the previous step for the resolution of the Poisson equation for lnq_auto.
Definition star.h:628
Lorene::Star_bin::conf_flat
bool conf_flat
true if the 3-metric is conformally flat, false for a more general metric.
Definition star.h:681
Lorene::Star::logn
Scalar logn
Logarithm of the lapse N .
Definition star.h:222
Lorene::Star::nn
Scalar nn
Lapse function N .
Definition star.h:225
Lorene::Star::ener_euler
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition star.h:198
Lorene::Star::equation_of_state
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition star.C:462
Lorene::Star::ray_eq
double ray_eq() const
Coordinate radius at , [r_unit].
Definition star_global.C:108
Lorene::Star::s_euler
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition star.h:201
Lorene::Star::stress_euler
Sym_tensor stress_euler
Spatial part of the stress-energy tensor with respect to the Eulerian observer.
Definition star.h:212
Lorene::Star::press
Scalar press
Fluid pressure.
Definition star.h:194
Lorene::Star::ray_eq_pi
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Definition star_global.C:186
Lorene::Star::ent
Scalar ent
Log-enthalpy.
Definition star.h:190
Lorene::Star::u_euler
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition star.h:207
Lorene::Star::mp
Map & mp
Mapping associated with the star.
Definition star.h:180
Lorene::Star::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition star.h:183
Lorene::Star::beta
Vector beta
Shift vector.
Definition star.h:228
Lorene::Star::ray_pole
double ray_pole() const
Coordinate radius at [r_unit].
Definition star_global.C:278
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Lorene::Tenseur::set
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
Lorene::Tenseur::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
Lorene::Tenseur::poisson_vect_oohara
void poisson_vect_oohara(double lambda, Param &par, Tenseur &shift, Tenseur &scal) const
Solves the vectorial Poisson equation .
Definition tenseur_pde.C:218
Lorene::Tenseur::set_std_base
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1170
Lorene::Tensor
Tensor handling.
Definition tensor.h:288
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::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::exp
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
Lorene::norme
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:481
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition cmp_math.C:410
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::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::Tensor::set_etat_zero
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition tensor.C:497
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::phi
Coord phi
coordinate centered on the grid
Definition map.h:720
Unites
Standard units of space, time and mass.