star_bin_vel_pot.C
/*
* Method of class Star_bin to compute the velocity scalar potential $\psi$
* by solving the continuity equation.
*
* (see file star.h for documentation).
*
*/
/*
* Copyright (c) 2004 Francois Limousin
* 2000-2001 Eric Gourgoulhon (for precedent version)
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
char star_bin_vel_pot_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_bin_vel_pot.C,v 1.7 2014/10/13 08:53:39 j_novak Exp $" ;
/*
* $Id: star_bin_vel_pot.C,v 1.7 2014/10/13 08:53:39 j_novak Exp $
* $Log: star_bin_vel_pot.C,v $
* Revision 1.7 2014/10/13 08:53:39 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.6 2005/09/13 19:38:31 f_limousin
* Reintroduction of the resolution of the equations in cartesian coordinates.
*
* Revision 1.5 2005/02/24 16:07:57 f_limousin
* Change the name of some variables (for instance dcov_logn --> dlogn).
*
* Revision 1.4 2005/02/17 17:33:11 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.3 2004/06/22 12:53:09 f_limousin
* Change qq, qq_auto and qq_comp to beta, beta_auto and beta_comp.
*
* Revision 1.2 2004/06/07 16:26:01 f_limousin
* Many modif...
*
* Revision 1.1 2004/04/08 16:34:08 f_limousin
* First version
*
*
* $Header: /cvsroot/Lorene/C++/Source/Star/star_bin_vel_pot.C,v 1.7 2014/10/13 08:53:39 j_novak Exp $
*
*/
// Headers Lorene
#include "star.h"
#include "eos.h"
#include "param.h"
#include "cmp.h"
#include "tenseur.h"
#include "graphique.h"
#include "utilitaires.h"
// Local prototype
namespace Lorene {
Cmp raccord_c1(const Cmp& uu, int l1) ;
double Star_bin::velocity_potential(int mermax, double precis, double relax) {
int nzm1 = mp.get_mg()->get_nzone() - 1 ;
//----------------------------------
// Specific relativistic enthalpy ---> hhh
//----------------------------------
Scalar hhh = exp(ent) ; // = 1 at the Newtonian limit
hhh.std_spectral_base() ;
//----------------------------------------------
// Computation of W^i = - A^2 h Gamma_n B^i/N
// See Eq (62) from Gourgoulhon et al. (2001)
//----------------------------------------------
Vector www = hhh * gam_euler * bsn * psi4 ;
www.change_triad( mp.get_bvect_cart() ) ; // components on the mapping
// Cartesian basis
//-------------------------------------------------
// Constant value of W^i at the center of the star
//-------------------------------------------------
Vector v_orb(mp, CON, mp.get_bvect_cart()) ;
v_orb.set_etat_qcq() ;
for (int i=1; i<=3; i++) {
v_orb.set(i) = www(i).val_grid_point(0, 0, 0, 0) ;
}
v_orb.set_triad( *(www.get_triad()) ) ;
v_orb.std_spectral_base() ;
//-------------------------------------------------
// Source and coefficients a,b for poisson_compact (independent from psi0)
//-------------------------------------------------
Scalar dndh_log = eos.der_nbar_ent(ent, nzet) ;
// In order to avoid any division by zero in the computation of zeta_h
// the value of dndh_log is set to 1 in the external domains:
for (int l=nzet; l <= nzm1; l++) {
dndh_log.set_domain(l) = 1 ;
}
Scalar zeta_h( ent / dndh_log ) ;
zeta_h.std_spectral_base() ;
Metric_flat flat_spher (mp.flat_met_spher()) ;
Vector bb = (1 - zeta_h) * ent.derive_con(flat_spher) +
zeta_h * lnq.derive_con(flat_spher) ;
Scalar entmb = ent - lnq ;
www.change_triad(mp.get_bvect_cart()) ;
v_orb.change_triad(mp.get_bvect_cart()) ;
Tensor dcovdcov_psi0 = psi0.derive_cov(flat).derive_cov(flat) ;
dcovdcov_psi0.inc_dzpuis() ;
// See Eq (63) from Gourgoulhon et al. (2001)
Scalar source = contract(www - v_orb, 0, ent.derive_cov(flat), 0)
+ zeta_h * ( contract(v_orb, 0, entmb.derive_cov(flat), 0) +
contract(www/gam_euler, 0, gam_euler.derive_cov(flat), 0)
+ contract(hij, 0, 1, (psi0.derive_cov(flat)
+ v_orb.down(0, flat))*entmb.derive_cov(flat),
0, 1)
- contract(hij, 0, 1, dcovdcov_psi0, 0, 1))
- contract(hij, 0, 1, ent.derive_cov(flat) * (psi0.derive_cov(flat)
+ v_orb.down(0, flat)),
0, 1) ;
/*
des_meridian(zeta_h,0., 4., "zeta_h", 10) ;
arrete() ;
des_meridian(bb(1),0., 4., "bb(1)", 10) ;
arrete() ;
des_meridian(bb(2),0., 4., "bb(2)", 10) ;
arrete() ;
des_meridian(bb(3),0., 4., "bb(3)", 10) ;
arrete() ;
des_meridian(psi0,0., 4., "psi0", 10) ;
arrete() ;
des_meridian(source,0., 4., "source", 10) ;
arrete() ;
*/
www.change_triad(mp.get_bvect_cart()) ;
v_orb.change_triad(mp.get_bvect_cart()) ;
source.annule(nzet, nzm1) ;
//---------------------------------------------------
// Resolution by means of Map_radial::poisson_compact
//---------------------------------------------------
Param par ;
int niter ;
par.add_int(mermax) ;
par.add_double(precis, 0) ;
par.add_double(relax, 1) ;
par.add_int_mod(niter) ;
if (psi0.get_etat() == ETATZERO) {
psi0 = 0 ;
}
Cmp source_cmp (source) ;
Cmp zeta_h_cmp (zeta_h) ;
Cmp psi0_cmp (psi0) ;
Tenseur bb_cmp(mp, 1, CON, mp.get_bvect_spher()) ;
bb_cmp.set_etat_qcq() ;
Cmp bb_cmp1 (bb(1)) ;
Cmp bb_cmp2 (bb(2)) ;
Cmp bb_cmp3 (bb(3)) ;
bb_cmp.set(0) = bb_cmp1 ;
bb_cmp.set(1) = bb_cmp2 ;
bb_cmp.set(2) = bb_cmp3 ;
source_cmp.va.ylm() ;
cout << "source" << endl << norme(source_cmp)<< endl ;
cout << "zeta_h " << endl << norme(zeta_h_cmp) << endl ;
cout << "bb(1)" << endl << norme(bb_cmp(0)) << endl ;
cout << "bb(2)" << endl << norme(bb_cmp(1)) << endl ;
cout << "bb(3)" << endl << norme(bb_cmp(2)) << endl ;
cout << "psiO" << endl << norme(psi0_cmp) << endl ;
mp.poisson_compact(source_cmp, zeta_h_cmp, bb_cmp, par, psi0_cmp ) ;
psi0 = psi0_cmp ;
cout << "psiO apres" << endl << norme(psi0) << endl ;
//---------------------------------------------------
// Check of the solution
//---------------------------------------------------
Scalar bb_dpsi0 = contract(bb, 0, psi0.derive_cov(flat_spher), 0) ;
Scalar oper = zeta_h * psi0.laplacian() + bb_dpsi0 ;
source.set_spectral_va().ylm_i() ;
double erreur = diffrel(oper, source)(0) ;
cout << "Check of the resolution of the continuity equation : "
<< endl ;
cout << " norme(source) : " << norme(source)(0)
<< " diff oper/source : " << erreur << endl ;
//--------------------------------
// Computation of grad(psi)
//--------------------------------
v_orb.change_triad(mp.get_bvect_cart()) ;
d_psi.change_triad(mp.get_bvect_cart()) ;
for (int i=1; i<=3; i++)
d_psi.set(i) = (psi0.derive_cov(flat))(i) + v_orb(i) ;
v_orb.change_triad(mp.get_bvect_cart()) ;
d_psi.change_triad(mp.get_bvect_cart()) ;
// C^1 continuation of d_psi outside the star
// (to ensure a smooth enthalpy field accross the stellar surface)
// ----------------------------------------------------------------
d_psi.annule(nzet, nzm1) ;
for (int i=1; i<=3; i++) {
Cmp d_psi_i (d_psi(i)) ;
d_psi_i = raccord_c1(d_psi_i, nzet) ;
d_psi.set(i) = d_psi_i ;
}
return erreur ;
}
}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp.
Definition cmp.h:464
Lorene::Metric_flat
Flat metric for tensor calculation.
Definition metric.h:261
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_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_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::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::derive_cov
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this.
Definition scalar_deriv.C:390
Lorene::Scalar::std_spectral_base
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
Lorene::Scalar::set_domain
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition scalar.h:615
Lorene::Scalar::annule
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:604
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::Star_bin::bsn
Vector bsn
3-vector shift, divided by N, of the rotating coordinates, .
Definition star.h:518
Lorene::Star_bin::d_psi
Vector d_psi
Gradient of (in the irrotational case) (Spherical components with respect to the mapping of the star...
Definition star.h:501
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::psi0
Scalar psi0
Scalar potential of the non-translational part of fluid 4-velocity (in the irrotational case).
Definition star.h:496
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::psi4
Scalar psi4
Conformal factor .
Definition star.h:552
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::eos
const Eos & eos
Equation of state of the stellar matter.
Definition star.h:185
Lorene::Star::gam_euler
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition star.h:204
Lorene::Star::ent
Scalar ent
Log-enthalpy.
Definition star.h:190
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::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::Tensor
Tensor handling.
Definition tensor.h:288
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm().
Definition valeur_ylm_i.C:131
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::change_triad
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition vector_change_triad.C:75
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
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::Tensor::set_triad
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition tensor.C:519
Lorene::Tensor::down
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.
Definition tensor_calculus.C:265
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::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:866
Lorene::Tensor::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:481
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