star_bin_equilibrium_xcts.C
/*
* Method of class Star_bin to compute an equilibrium configuration
* (see file star.h for documentation).
*/
/*
* Copyright (c) 2010 Michal Bejger
*
* 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_xcts_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_bin_equilibrium_xcts.C,v 1.13 2014/10/13 08:53:38 j_novak Exp $" ;
/*
* $Id: star_bin_equilibrium_xcts.C,v 1.13 2014/10/13 08:53:38 j_novak Exp $
* $Log: star_bin_equilibrium_xcts.C,v $
* Revision 1.13 2014/10/13 08:53:38 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.12 2014/10/06 15:13:16 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.11 2011/03/25 16:28:12 e_gourgoulhon
* Still in progress
*
* Revision 1.10 2010/12/20 15:42:10 m_bejger
* Various rearrangements of fields in Poissson equations
*
* Revision 1.9 2010/12/14 17:34:42 m_bejger
* Improved iteration for beta_auto poisson_vect()
*
* Revision 1.8 2010/12/09 10:48:06 m_bejger
* Testing the main equations
*
* Revision 1.7 2010/10/26 18:46:28 m_bejger
* Added table fact_resize for domain resizing
*
* Revision 1.6 2010/10/18 19:08:14 m_bejger
* Changed to allow for calculations with more than one domain in the star
*
* Revision 1.5 2010/06/23 20:40:56 m_bejger
* Corrections in equations for Psi_auto, chi_auto and beta_auto
*
* Revision 1.4 2010/06/18 13:28:59 m_bejger
* Adjusted the computation of the first integral, radial scale
*
* Revision 1.3 2010/06/17 17:05:06 m_bejger
* Testing version
*
* Revision 1.2 2010/06/15 08:21:21 m_bejger
* Minor changes; still not working properly
*
* Revision 1.1 2010/05/04 07:51:05 m_bejger
* Initial version
*
* $Header: /cvsroot/Lorene/C++/Source/Star/star_bin_equilibrium_xcts.C,v 1.13 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_xcts::equilibrium(double ent_c,
int mermax,
int mermax_potvit,
int mermax_poisson,
double relax_poisson,
double relax_potvit,
double thres_adapt,
const Tbl& fact_resize,
const Tbl* pent_limit,
Tbl& diff) {
// 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_psi = diff.set(2) ;
double& diff_chi = diff.set(3) ;
double& diff_beta_x = diff.set(4) ;
double& diff_beta_y = diff.set(5) ;
double& diff_beta_z = diff.set(6) ;
// Parameters for the 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
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
// Enthalpy values for the adaptation
Tbl ent_limit(nzet) ;
if (pent_limit != 0x0) ent_limit = *pent_limit ;
par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
// to define the isosurfaces.
double precis_poisson = 1.e-16 ;
Cmp ssjm1psi (ssjm1_psi) ;
Cmp ssjm1chi (ssjm1_chi) ;
// Parameters for the function Scalar::poisson for Psi_auto
// ---------------------------------------------------------------
Param par_psi ;
par_psi.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_psi.add_double(relax_poisson, 0) ; // relaxation parameter
par_psi.add_double(precis_poisson, 1) ; // required precision
par_psi.add_int_mod(niter, 0) ; // number of iterations actually used
par_psi.add_cmp_mod( ssjm1psi ) ;
// Parameters for the function Scalar::poisson for chi_auto
// ---------------------------------------------------------------
Param par_chi ;
par_chi.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_chi.add_double(relax_poisson, 0) ; // relaxation parameter
par_chi.add_double(precis_poisson, 1) ; // required precision
par_chi.add_int_mod(niter, 0) ; // number of iterations actually used
par_chi.add_cmp_mod( ssjm1chi ) ;
// Parameters for the function Vector::poisson for beta
// ----------------------------------------------------
Param par_beta ;
par_beta.add_int(mermax_poisson, 0) ; // maximum number of
// iterations
par_beta.add_double(relax_poisson, 0) ; // relaxation parameter
par_beta.add_double(precis_poisson, 1) ; // required precision
par_beta.add_int_mod(niter, 0) ; // number of iterations
// actually used
// Sources at the previous step, for a poisson_vect() solver
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_beta.add_cmp_mod(ssjm1khi) ;
par_beta.add_tenseur_mod(ssjm1wbeta) ;
// Redefinition of external potential
// See Eq (99) from Gourgoulhon et al. (2001)
// logN = logN_auto + logn_ac_rest = log(chi_auto + 1.)
// - log(Psi_auto + 1.) + logn_ac_rest
//------------------------------------
Scalar Psi_auto_p1 = Psi_auto + 1. ;
Scalar chi_auto_p1 = chi_auto + 1. ;
Scalar logn_auto = log(chi_auto_p1) - log(Psi_auto_p1) ;
logn_auto.std_spectral_base() ;
Scalar logn_ac_rest = log(1. + chi_comp/chi_auto_p1)
- log(1. + Psi_comp/Psi_auto_p1) ;
logn_ac_rest.std_spectral_base() ;
//cout << "logn_auto" << norme(logn_auto) << endl ;
//cout << "logn_ac_rest" << norme(logn_ac_rest) << endl ;
//cout << "pot_centri" << norme(pot_centri) << endl ;
//cout << "loggam" << norme(loggam) << endl ;
Scalar pot_ext = logn_ac_rest + pot_centri + loggam ;
//cout << "pot_ext" << norme(pot_ext) << endl ;
Scalar ent_jm1 = ent ; // Enthalpy at previous step
Scalar source_tot(mp) ; // source term in equations Psi_auto
// and chi_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 ) ;
// cout << "k, j, alpha_r2_jk : " << k << " " << j << " " << alpha_r2_jk << endl ;
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
// ----------------------
Psi_auto = pow(Psi_auto +1.,alpha_r2) - 1. ;
chi_auto = pow(chi_auto +1.,alpha_r2) - 1. ;
Psi_auto.std_spectral_base() ;
chi_auto.std_spectral_base() ;
//------------------------------------------------------------
// Change the values of the inner points of the domain adjascent
// to the surface of the star by those of the outer points of
// the domain under the surface
//------------------------------------------------------------
Psi_auto.set_spectral_va().smooth(nzet, Psi_auto.set_spectral_va()) ;
chi_auto.set_spectral_va().smooth(nzet, chi_auto.set_spectral_va()) ;
logn_auto = log(chi_auto + 1.) - log(Psi_auto + 1.) ;
logn_auto.std_spectral_base() ;
logn_auto_c = logn_auto.val_grid_point(0, 0, 0, 0) ;
//------------------------------------------
// 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
double rr_in_1 = mp.val_r(nzet, -1., M_PI/2., 0.) ;
// Resizes the outer boundary of the shell including the comp. NS
double rr_out_nm2 = mp.val_r(nz-2, 1., M_PI/2., 0.) ;
mp.resize(nz-2, rr_in_1/rr_out_nm2 * fact_resize(1)) ;
// Resizes the inner boundary of the shell including the comp. NS
double rr_out_nm3 = mp.val_r(nz-3, 1., M_PI/2., 0.) ;
mp.resize(nz-3, rr_in_1/rr_out_nm3 * fact_resize(0)) ;
if (nz > nzet+3) {
// Resize of the domain from 1(nzet) to N-4
double rr_out_nm3_new = mp.val_r(nz-3, 1., M_PI/2., 0.) ;
for (int i=nzet-1; i<nz-4; i++) {
double rr_out_i = mp.val_r(i, 1., M_PI/2., 0.) ;
double rr_mid = rr_out_i
+ (rr_out_nm3_new - rr_out_i) / double(nz - 3 - i) ;
double rr_2timesi = 2. * rr_out_i ;
if (rr_2timesi < rr_mid) {
double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
mp.resize(i+1, rr_2timesi / rr_out_ip1) ;
} else {
double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
mp.resize(i+1, rr_mid / rr_out_ip1) ;
} // End of else
} // End of i loop
} // End of (nz > nzet+3) loop
// 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 ;
mp.reevaluate_symy(&mp_prev, nzet+1, ent) ;
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
// -------------------------------
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 Psi3 = psi4 / Psi ;
Psi3.std_spectral_base() ;
Scalar sPsi7 = Psi * pow(psi4, -2.) ;
sPsi7.std_spectral_base() ;
//------------------------------------------------------------------
// Poisson equation for Psi_auto (Eq. 8.127 of arXiv:gr-qc/0703035)
//------------------------------------------------------------------
// Source
//--------
source_tot = - 0.5 * qpig * psi4 % Psi % ener_euler
- 0.125 * sPsi7 * (hacar_auto + hacar_comp) ;
source_tot.std_spectral_base() ;
// Resolution of the Poisson equation
// ----------------------------------
cout << "Resolution of the Poisson equation for Psi_auto:" << endl ;
source_tot.poisson(par_psi, Psi_auto) ;
ssjm1_psi = ssjm1psi ;
cout << "Psi_auto: " << endl << norme(Psi_auto/(nr*nt*np)) << endl ;
// Check: has the Poisson equation been correctly solved ?
// -----------------------------------------------------
Tbl tdiff_psi = diffrel(Psi_auto.laplacian(), source_tot) ;
cout <<
"Relative error in the resolution of the equation for Psi_auto : "
<< endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_psi(l) << " " ;
}
cout << endl
<< "==========================================================="
<< endl ;
diff_psi = max(abs(tdiff_psi)) ;
//------------------------------------------------------------------
// Poisson equation for chi_auto (Eq. 8.129 of arXiv:gr-qc/0703035)
//------------------------------------------------------------------
// Source
//--------
source_tot = chi * 0.5 * qpig * psi4 % (ener_euler + 2.*s_euler)
+ 0.875 * nn % sPsi7 * (hacar_auto + hacar_comp) ;
source_tot.std_spectral_base() ;
// Resolution of the Poisson equation
// ----------------------------------
cout << "Resolution of the Poisson equation for chi_auto:" << endl ;
source_tot.poisson(par_chi, chi_auto) ;
ssjm1_chi = ssjm1chi ;
cout << "chi_auto: " << endl << norme(chi_auto/(nr*nt*np)) << endl ;
// Check: has the Poisson equation been correctly solved ?
// -----------------------------------------------------
Tbl tdiff_chi = diffrel(chi_auto.laplacian(), source_tot) ;
cout <<
"Relative error in the resolution of the equation for chi_auto : "
<< endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_chi(l) << " " ;
}
cout << endl
<< "==========================================================="
<< endl ;
diff_chi = max(abs(tdiff_chi)) ;
//------------------------------------------------------------------
// Vector Poisson equation for beta_auto
// (Eq. 8.128 of arXiv:gr-qc/0703035)
//------------------------------------------------------------------
// Source
//--------
source_beta = 4.* qpig * chi % Psi3 * (ener_euler + press) * u_euler
+ 2. * sPsi7 * contract(haij_auto, 1, dcov_chi, 0)
-14. * nn % sPsi7 * contract(haij_auto, 1, dcov_Psi, 0) ;
source_beta.std_spectral_base() ;
// Resolution of the Poisson equation
// ----------------------------------
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)) ;
//source_p.set(i).filtre(4) ;
}
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) = Cmp(w_beta(i+1)) ;
vect_auxi.set_std_base() ;
Tenseur scal_auxi (mp) ;
scal_auxi.set_etat_qcq() ;
scal_auxi.set() = Cmp(khi) ;
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) = Cmp(beta_auto.set(i+1));
resu_p.set_std_base() ;
// Resolution of the vector Poisson equation
// -----------------------------------------
double lambda = double(1) / double(3) ;
source_p.poisson_vect(lambda, par_beta, resu_p, vect_auxi, scal_auxi) ;
//source_p.poisson_vect_oohara(lambda, par_beta, resu_p, scal_auxi) ;
for (int i=1; i<=3; i++) {
beta_auto.set(i) = resu_p(i-1) ;
w_beta.set(i) = vect_auxi(i-1) ;
}
khi = scal_auxi() ;
// Values of sources from the previous step
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)) ;
// 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::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::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::Star_bin_xcts::Psi
Scalar Psi
Total conformal factor .
Definition star.h:1152
Lorene::Star_bin_xcts::pot_centri
Scalar pot_centri
Centrifugal potential.
Definition star.h:1129
Lorene::Star_bin_xcts::haij_auto
Sym_tensor haij_auto
Part of the extrinsic curvature tensor generated by beta_auto.
Definition star.h:1193
Lorene::Star_bin_xcts::psi4
Scalar psi4
Conformal factor .
Definition star.h:1158
Lorene::Star_bin_xcts::flat
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) .
Definition star.h:1177
Lorene::Star_bin_xcts::equilibrium
void equilibrium(double ent_c, int mermax, int mermax_potvit, int mermax_poisson, double relax_poisson, double relax_potvit, double thres_adapt, const Tbl &fact, const Tbl *pent_limit, Tbl &diff)
Computes an equilibrium configuration.
Definition star_bin_equilibrium_xcts.C:91
Lorene::Star_bin_xcts::ssjm1_khi
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for khi.
Definition star.h:1227
Lorene::Star_bin_xcts::khi
Scalar khi
Solution for the scalar part of the vector Poisson equation for .
Definition star.h:1168
Lorene::Star_bin_xcts::ssjm1_wbeta
Vector ssjm1_wbeta
Effective source at the previous step for wbeta of the vector Poisson equation for wbeta (needed for ...
Definition star.h:1234
Lorene::Star_bin_xcts::Psi_comp
Scalar Psi_comp
Scalar field generated principally by the companion star.
Definition star.h:1149
Lorene::Star_bin_xcts::ssjm1_chi
Scalar ssjm1_chi
Effective source at the previous step for the resolution of the Poisson equation for \chi_auto .
Definition star.h:1216
Lorene::Star_bin_xcts::dcov_Psi
Vector dcov_Psi
Covariant derivative of the conformal factor .
Definition star.h:1171
Lorene::Star_bin_xcts::beta_auto
Vector beta_auto
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Definition star.h:1182
Lorene::Star_bin_xcts::Psi_auto
Scalar Psi_auto
Scalar field generated principally by the star.
Definition star.h:1144
Lorene::Star_bin_xcts::hacar_comp
Scalar hacar_comp
Part of the scalar generated by beta_auto and beta_comp, i.e.
Definition star.h:1211
Lorene::Star_bin_xcts::chi
Scalar chi
Total function .
Definition star.h:1155
Lorene::Star_bin_xcts::loggam
Scalar loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition star.h:1120
Lorene::Star_bin_xcts::dcov_chi
Vector dcov_chi
Covariant derivative of the function .
Definition star.h:1174
Lorene::Star_bin_xcts::chi_auto
Scalar chi_auto
Scalar field generated principally by the star.
Definition star.h:1134
Lorene::Star_bin_xcts::w_beta
Vector w_beta
Solution for the vector part of the vector Poisson equation for .
Definition star.h:1163
Lorene::Star_bin_xcts::chi_comp
Scalar chi_comp
Scalar field generated principally by the companion star.
Definition star.h:1139
Lorene::Star_bin_xcts::ssjm1_psi
Scalar ssjm1_psi
Effective source at the previous step for the resolution of the Poisson equation for \Psi_auto .
Definition star.h:1221
Lorene::Star_bin_xcts::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_xcts.C:49
Lorene::Star_bin_xcts::irrotational
bool irrotational
true for an irrotational star, false for a corotating one
Definition star.h:1099
Lorene::Star_bin_xcts::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_xcts.C:70
Lorene::Star_bin_xcts::hacar_auto
Scalar hacar_auto
Part of the scalar generated by beta_auto, i.e.
Definition star.h:1205
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::press
Scalar press
Fluid pressure.
Definition star.h:194
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::ray_pole
double ray_pole() const
Coordinate radius at [r_unit].
Definition star_global.C:278
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::set_std_base
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1170
Lorene::Tenseur::poisson_vect
void poisson_vect(double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
Definition tenseur_pde.C:118
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::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
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::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
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::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
Unites
Standard units of space, time and mass.