et_bin_nsbh_equilibrium.C
/*
* Method of class Etoile to compute a static spherical configuration
* of a neutron star in a NS-BH binary system.
*
* (see file etoile.h for documentation).
*
*/
/*
* Copyright (c) 2003 Keisuke Taniguchi
*
* 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 et_bin_nsbh_equilibrium_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_nsbh_equilibrium.C,v 1.13 2014/10/13 08:52:56 j_novak Exp $" ;
/*
* $Id: et_bin_nsbh_equilibrium.C,v 1.13 2014/10/13 08:52:56 j_novak Exp $
* $Log: et_bin_nsbh_equilibrium.C,v $
* Revision 1.13 2014/10/13 08:52:56 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.12 2014/10/06 15:13:08 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.11 2008/09/26 08:38:45 p_grandclement
* get rid of desaliasing
*
* Revision 1.10 2006/09/05 13:39:45 p_grandclement
* update of the bin_ns_bh project
*
* Revision 1.9 2006/06/01 12:47:53 p_grandclement
* update of the Bin_ns_bh project
*
* Revision 1.8 2006/04/25 07:21:58 p_grandclement
* Various changes for the NS_BH project
*
* Revision 1.7 2006/03/30 07:33:47 p_grandclement
* *** empty log message ***
*
* Revision 1.6 2005/10/18 13:12:33 p_grandclement
* update of the mixted binary codes
*
* Revision 1.5 2005/08/29 15:10:17 p_grandclement
* Addition of things needed :
* 1) For BBH with different masses
* 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
* WORKING YET !!!
*
* Revision 1.4 2004/03/25 10:29:04 j_novak
* All LORENE's units are now defined in the namespace Unites (in file unites.h).
*
* Revision 1.3 2003/10/24 12:34:06 k_taniguchi
* Change the notation as it should be
*
* Revision 1.2 2003/10/21 11:49:33 k_taniguchi
* Change the class from Etoile_bin to sub-class Et_bin_nsbh.
*
* Revision 1.1 2003/10/20 15:01:55 k_taniguchi
* Computation of an equilibrium configuration of a neutron star
* in a NS-BH binary system.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_nsbh_equilibrium.C,v 1.13 2014/10/13 08:52:56 j_novak Exp $
*
*/
// Headers C
#include <cmath>
// Headers Lorene
#include "etoile.h"
#include "map.h"
#include "nbr_spx.h"
#include "et_bin_nsbh.h"
#include "param.h"
#include "graphique.h"
#include "utilitaires.h"
#include "unites.h"
namespace Lorene {
void Et_bin_nsbh::equilibrium_nsbh(bool adapt, double ent_c, int& niter, int mermax,
int mermax_poisson, double relax_poisson,
int mermax_potvit, double relax_potvit,
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) ;
// Error indicators
// ----------------
double& diff_ent = diff.set(0) ;
double& diff_vel_pot = diff.set(1) ;
double& diff_lapse = diff.set(2) ;
double& diff_confpsi = diff.set(3) ;
double& diff_shift_x = diff.set(4) ;
double& diff_shift_y = diff.set(5) ;
double& diff_shift_z = diff.set(6) ;
// Parameters for the function Map_et::poisson for n_auto
// ------------------------------------------------------
double precis_poisson = 1.e-16 ;
Param par_poisson1 ;
par_poisson1.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_poisson1.add_double(relax_poisson, 0) ; // relaxation parameter
par_poisson1.add_double(precis_poisson, 1) ; // required precision
par_poisson1.add_int_mod(niter, 0) ; // number of iterations actually used
par_poisson1.add_cmp_mod( ssjm1_lapse ) ;
// Parameters for the function Map_et::poisson for confpsi_auto
// ------------------------------------------------------------
Param par_poisson2 ;
par_poisson2.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_poisson2.add_double(relax_poisson, 0) ; // relaxation parameter
par_poisson2.add_double(precis_poisson, 1) ; // required precision
par_poisson2.add_int_mod(niter, 0) ; // number of iterations actually used
par_poisson2.add_cmp_mod( ssjm1_confpsi ) ;
// Parameters for the function Tenseur::poisson_vect
// -------------------------------------------------
Param par_poisson_vect ;
par_poisson_vect.add_int(mermax_poisson, 0) ;
// maximum number of iterations
par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
par_poisson_vect.add_cmp_mod( ssjm1_khi ) ;
par_poisson_vect.add_tenseur_mod( ssjm1_wshift ) ;
par_poisson_vect.add_int_mod(niter, 0) ;
// Parameters for the adaptation
Param par_adapt ;
int nitermax = 100 ;
int niter_adapt ;
int adapt_flag = (adapt) ? 1 : 0 ;
int nz_search = nzet + 1 ;
double precis_secant = 1.e-14 ;
double alpha_r ;
double reg_map = 1. ;
int k_b ;
int j_b ;
Tbl ent_limit(nzet) ;
par_adapt.add_int(nitermax, 0) ;
par_adapt.add_int(nzet, 1) ;
par_adapt.add_int(nz_search, 2) ;
par_adapt.add_int(adapt_flag, 3) ;
par_adapt.add_int(j_b, 4) ;
par_adapt.add_int(k_b, 5) ;
par_adapt.add_int_mod(niter_adapt, 0) ;
par_adapt.add_double(precis_secant, 0) ;
par_adapt.add_double(reg_map, 1) ;
par_adapt.add_double(alpha_r, 2) ;
par_adapt.add_tbl(ent_limit, 0) ;
// External potential
// See Eq (99) from Gourgoulhon et al. (2001)
// -----------------------------------------
Tenseur ent_jm1 = ent ; // Enthalpy at previous step
Tenseur source(mp) ; // source term in the equation for logn_auto
// and beta_auto
Tenseur source_shift(mp, 1, CON, ref_triad) ; // source term in the
// equation for shift_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) ;
}
// Equation de la surface
//if (adapt) {
// Rescaling of the radius : (Be carefull !)
int nt = mg->get_nt(nzet-1) ;
int np = mg->get_np(nzet-1) ;
int nr = mg->get_nr(nzet-1) ;
// valeurs au centre
double hc = exp(ent_c) ;
double gamma_c = exp(loggam())(0,0,0,0) ;
double gamma_0_c = exp(-pot_centri())(0,0,0,0) ;
double n_auto_c = n_auto()(0,0,0,0) ;
double n_comp_c = n_comp()(0,0,0,0) ;
double alpha_square = 0 ;
double constante = 0;
for (int k=0; k<np; k++) {
for (int j=0; j<nt; j++) {
// valeurs au bord
double gamma_b = exp(loggam())(nzet-1,k,j,nr-1) ;
double gamma_0_b = exp(-pot_centri())(nzet-1,k,j,nr-1) ;
double n_auto_b = n_auto()(nzet-1,k,j,nr-1) ;
double n_comp_b = n_comp()(nzet-1,k,j,nr-1) ;
// Les solutions :
double alpha_square_courant = (gamma_0_c*gamma_b*n_comp_b - hc*gamma_c*gamma_0_b*n_comp_c) /
(hc*gamma_c*gamma_0_b*n_auto_c-gamma_0_c*gamma_b*n_auto_b) ;
double constante_courant = gamma_b*(n_comp_b+alpha_square_courant*n_auto_b)/gamma_0_b ;
if (alpha_square_courant > alpha_square) {
alpha_square = alpha_square_courant ;
k_b = k ;
j_b = j ;
constante = constante_courant ;
}
}
}
alpha_r = sqrt(alpha_square) ;
cout << "Adaptation : " << k_b << " " << j_b << " " << alpha_r << endl ;
// Le potentiel :
Tenseur potentiel (constante*exp(-loggam-pot_centri)/(n_auto*alpha_square+n_comp)) ;
potentiel.set_std_base() ;
for (int l=nzet+1 ; l<nz ; l++)
potentiel.set().va.set(l) = 1 ;
Map_et mp_prev = mp_et ;
ent = log(potentiel) ;
ent.set_std_base() ;
ent().va.smooth(nzet, (ent.set().va)) ;
ent_limit.set_etat_qcq() ;
for (int l=0; l<nzet; l++) { // loop on domains inside the star
ent_limit.set(l) = ent()(l, k_b, j_b, nr-1) ;
}
// On adapte :
mp.adapt(ent(), par_adapt, 4) ;
mp_prev.homothetie(alpha_r) ;
for (int l=nzet ; l<nz-1 ; l++)
mp.resize(l, 1./alpha_r) ;
mp.reevaluate_symy (&mp_prev, nzet, ent.set()) ;
// 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)
//-------------------------------------------------
// 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 ;
//--------------------------------------------------------
// Poisson equation for n_auto
//--------------------------------------------------------
// Source
// See Eq (50) from Gourgoulhon et al. (2001)
// ------------------------------------------
Tenseur confpsi_q = pow(confpsi, 4.) ;
Tenseur confpsi_c = pow(confpsi, 5.) ;
if (relativistic) {
Tenseur tmp = flat_scalar_prod(tkij_tot, tkij_auto) ;
Tenseur kk (mp) ;
kk = 0 ;
Tenseur tmp2(mp) ;
tmp2.set_etat_qcq() ;
for (int i=0 ; i<3 ; i++) {
tmp2.set() = tmp(i, i) ;
kk = kk + tmp2 ;
}
source = qpig * nnn * confpsi_q * (ener_euler + s_euler)
+ nnn * confpsi_q * kk
- 2.*flat_scalar_prod(d_confpsi_auto+d_confpsi_comp, d_n_auto) /
confpsi ;
}
else {
cout <<
"WARNING : Et_bin_nsbh is for the relativistic calculation"
<< endl ;
abort() ;
}
source.set_std_base() ;
// Resolution of the Poisson equation
// ----------------------------------
Cmp n_auto_old (n_auto()) ;
source().poisson(par_poisson1, n_auto.set()) ;
n_auto.set() = n_auto() + 0.5 ;
// Difference pas précédent
// -----------------------------------------------------
Tbl tdiff_lapse = diffrel(n_auto(), n_auto_old) ;
cout <<
"Relative difference on n_auto : "
<< endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_lapse(l) << " " ;
}
cout << endl ;
diff_lapse = max(abs(tdiff_lapse)) ;
if (relativistic) {
//--------------------------------------------------------
// Poisson equation for confpsi_auto
//--------------------------------------------------------
// Source
// See Eq (51) from Gourgoulhon et al. (2001)
// ------------------------------------------
Tenseur tmp = flat_scalar_prod(tkij_tot, tkij_auto) ;
Tenseur kk (mp) ;
kk = 0 ;
Tenseur tmp2(mp) ;
tmp2.set_etat_qcq() ;
for (int i=0 ; i<3 ; i++) {
tmp2.set() = tmp(i, i) ;
kk = kk + tmp2 ;
}
source = -0.5 * qpig * confpsi_c * ener_euler
- 0.125 * confpsi_c * kk ;
source.set_std_base() ;
// Resolution of the Poisson equation
// ----------------------------------
Cmp psi_old (confpsi_auto()) ;
source().poisson(par_poisson2, confpsi_auto.set()) ;
confpsi_auto.set() = confpsi_auto() + 0.5 ;
// Check: has the Poisson equation been correctly solved ?
// -----------------------------------------------------
Tbl tdiff_confpsi = diffrel(confpsi_auto(), psi_old) ;
cout <<
"Relative difference on confpsi_auto : "
<< endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_confpsi(l) << " " ;
}
cout << endl ;
diff_confpsi = max(abs(tdiff_confpsi)) ;
//--------------------------------------------------------
// Vector Poisson equation for shift_auto
//--------------------------------------------------------
// Source
// See Eq (52) from Gourgoulhon et al. (2001)
// ------
Tenseur vtmp = d_n_auto -6. * nnn * d_confpsi_auto / confpsi ;
source_shift = 4.*qpig * nnn *confpsi_q *(ener_euler + press)
* u_euler ;
if (tkij_tot.get_etat() != ETATZERO)
source_shift = source_shift + 2.* flat_scalar_prod(tkij_tot, vtmp) ;
source_shift.set_std_base() ;
// Resolution of the Poisson equation
// ----------------------------------
// Filter for the source of shift vector
for (int i=0 ; i<3 ; i++)
if (source_shift(i).get_etat() != ETATZERO)
source_shift.set(i).va.coef_i() ;
for (int i=0; i<3; i++)
if ((source_shift(i).get_etat() != ETATZERO) && (source_shift(i).va.c->t[nz-1]->get_etat() != ETATZERO))
source_shift.set(i).filtre(4) ;
for (int i=0; i<3; i++) {
if(source_shift(i).dz_nonzero()) {
assert( source_shift(i).get_dzpuis() == 4 ) ;
}
else{
(source_shift.set(i)).set_dzpuis(4) ;
}
}
//##
// source_shift.dec2_dzpuis() ; // dzpuis 4 -> 2
double lambda_shift = double(1) / double(3) ;
// ON DOIT CHANGER DE TRIADE
source_shift.change_triad(mp.get_bvect_cart()) ;
Tenseur shift_old (shift_auto) ;
source_shift.poisson_vect_oohara(lambda_shift, par_poisson_vect,
shift_auto, khi_shift) ;
shift_auto.change_triad(ref_triad) ;
// Check: has the equation for shift_auto been correctly solved ?
// --------------------------------------------------------------
Tbl tdiff_shift_x = diffrel(shift_auto(0), shift_old(0)) ;
Tbl tdiff_shift_y = diffrel(shift_auto(1), shift_old(1)) ;
Tbl tdiff_shift_z = diffrel(shift_auto(2), shift_old(2)) ;
cout <<
"Relative difference on shift_auto : "
<< endl ;
cout << "x component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_shift_x(l) << " " ;
}
cout << endl ;
cout << "y component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_shift_y(l) << " " ;
}
cout << endl ;
cout << "z component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_shift_z(l) << " " ;
}
cout << endl ;
diff_shift_x = max(abs(tdiff_shift_x)) ;
diff_shift_y = max(abs(tdiff_shift_y)) ;
diff_shift_z = max(abs(tdiff_shift_z)) ;
} // End of relativistic equations
} // End of main loop
//=========================================================================
// End of iteration
//=========================================================================
}
// Truc pourri
void Et_bin_nsbh::equilibrium_nsbh (double, int, int, double,
int, double, double, const Tbl&, Tbl&) {
cout << "Not implemented !" << endl ;
abort() ;
}
}
Lorene::Et_bin_nsbh::confpsi_auto
Tenseur confpsi_auto
Part of the conformal factor $\Psi$ generated principaly by the star.
Definition et_bin_nsbh.h:104
Lorene::Et_bin_nsbh::d_confpsi_comp
Tenseur d_confpsi_comp
Gradient of {\tt confpsi_comp} (Cartesian components with respect to {\tt ref_triad}...
Definition et_bin_nsbh.h:119
Lorene::Et_bin_nsbh::ssjm1_confpsi
Cmp ssjm1_confpsi
Effective source at the previous step for the resolution of the Poisson equation for {\tt confpsi_aut...
Definition et_bin_nsbh.h:164
Lorene::Et_bin_nsbh::tkij_tot
Tenseur_sym tkij_tot
Total extrinsic curvature tensor $K^{ij}$ generated by {\tt shift_auto} and {\tt shift_comp}...
Definition et_bin_nsbh.h:152
Lorene::Et_bin_nsbh::equilibrium_nsbh
virtual void equilibrium_nsbh(double ent_c, int mermax, int mermax_poisson, double relax_poisson, int mermax_potvit, double relax_potvit, double thres_adapt, const Tbl &fact, Tbl &diff)
Computes an equilibrium configuration in a NS-BH binary system.
Definition et_bin_nsbh_equilibrium.C:494
Lorene::Et_bin_nsbh::d_n_auto
Tenseur d_n_auto
Gradient of {\tt n_auto} (Cartesian components with respect to {\tt ref_triad}).
Definition et_bin_nsbh.h:93
Lorene::Et_bin_nsbh::ssjm1_lapse
Cmp ssjm1_lapse
Effective source at the previous step for the resolution of the Poisson equation for {\tt n_auto}...
Definition et_bin_nsbh.h:158
Lorene::Et_bin_nsbh::tkij_auto
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor $K^{ij}$ generated by {\tt shift_auto}.
Definition et_bin_nsbh.h:146
Lorene::Et_bin_nsbh::n_auto
Tenseur n_auto
Part of the lapse {\it N} generated principaly by the star.
Definition et_bin_nsbh.h:85
Lorene::Et_bin_nsbh::n_comp
Tenseur n_comp
Part of the lapse {\it N} generated principaly by the companion star.
Definition et_bin_nsbh.h:88
Lorene::Et_bin_nsbh::confpsi
Tenseur confpsi
Total conformal factor $\Psi$.
Definition et_bin_nsbh.h:101
Lorene::Et_bin_nsbh::d_confpsi_auto
Tenseur d_confpsi_auto
Gradient of {\tt confpsi_auto} (Cartesian components with respect to {\tt ref_triad}...
Definition et_bin_nsbh.h:114
Lorene::Etoile_bin::ssjm1_khi
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition etoile.h:973
Lorene::Etoile_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 et_bin_vel_pot.C:142
Lorene::Etoile_bin::ref_triad
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur 's...
Definition etoile.h:828
Lorene::Etoile_bin::irrotational
bool irrotational
true for an irrotational star, false for a corotating one
Definition etoile.h:822
Lorene::Etoile_bin::hydro_euler
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Definition et_bin_hydro.C:106
Lorene::Etoile_bin::ssjm1_wshift
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for by means...
Definition etoile.h:983
Lorene::Etoile_bin::pot_centri
Tenseur pot_centri
Centrifugal potential.
Definition etoile.h:953
Lorene::Etoile_bin::loggam
Tenseur loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition etoile.h:849
Lorene::Etoile_bin::shift_auto
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
Definition etoile.h:889
Lorene::Etoile_bin::khi_shift
Tenseur khi_shift
Scalar used in the decomposition of shift_auto , following Shibata's prescription [Prog.
Definition etoile.h:918
Lorene::Etoile::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
Lorene::Etoile::nnn
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Lorene::Etoile::equation_of_state
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition etoile.C:566
Lorene::Etoile::u_euler
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition etoile.h:474
Lorene::Etoile::mp
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Lorene::Etoile::press
Tenseur press
Fluid pressure.
Definition etoile.h:461
Lorene::Etoile::relativistic
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition etoile.h:437
Lorene::Etoile::ener_euler
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition etoile.h:465
Lorene::Etoile::s_euler
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition etoile.h:468
Lorene::Etoile::ent
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case).
Definition etoile.h:457
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tenseur::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
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::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::flat_scalar_prod
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
Definition tenseur_operateur.C:653
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::adapt
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.