eos_consistent.C
/*
* Methods of class Eos_consistent
*
* (see file eos_compose.h for documentation).
*
*/
/*
* Copyright (c) 2015 Jerome Novak
*
* 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 eos_consistent_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/eos_consistent.C,v 1.1 2015/08/04 14:41:29 j_novak Exp $ " ;
/*
* $Id: eos_consistent.C,v 1.1 2015/08/04 14:41:29 j_novak Exp $
* $Log: eos_consistent.C,v $
* Revision 1.1 2015/08/04 14:41:29 j_novak
* Back to previous version for Eos_CompOSE. Enthalpy-consistent EoS can be accessed using Eos_consistent class (derived from Eos_CompOSE).
*
*
* $Header: /cvsroot/Lorene/C++/Source/Eos/eos_consistent.C,v 1.1 2015/08/04 14:41:29 j_novak Exp $
*
*/
#include <string>
// Headers Lorene
#include "headcpp.h"
#include "eos.h"
#include "tbl.h"
#include "utilitaires.h"
#include "unites.h"
namespace Lorene {
void interpol_herm(const Tbl& , const Tbl&, const Tbl&, double, int&,
double&, double& ) ;
//----------------------------//
// Constructors //
//----------------------------//
// Standard constructor
// --------------------
Eos_consistent::Eos_consistent(const char* file_name)
: Eos_CompOSE(file_name)
{}
// Constructor from binary file
// ----------------------------
Eos_consistent::Eos_consistent(FILE* fich) : Eos_CompOSE(fich)
{}
// Constructor from a formatted file
// ---------------------------------
Eos_consistent::Eos_consistent(ifstream& fich) : Eos_CompOSE(fich)
{}
// Constructor from CompOSE data files
// ------------------------------------
Eos_consistent::Eos_consistent(const string& files) : Eos_CompOSE(files)
{
using namespace Unites ;
// Files containing data and a test
//---------------------------------
string file_nb = files + ".nb" ;
string file_thermo = files + ".thermo" ;
ifstream in_nb(file_nb.data()) ;
// obtaining the size of the tables for memory allocation
//-------------------------------------------------------
int index1, index2 ;
in_nb >> index1 >> index2 ;
int nbp = index2 - index1 + 1 ;
assert( nbp == logh->get_taille() ) ;
press = new double[nbp] ;
nb = new double[nbp] ;
ro = new double[nbp] ;
// Variables and conversion
//-------------------------
double nb_fm3, rho_cgs, p_cgs, p_over_nb_comp, eps_comp ;
double dummy_x ;
int dummy_n ;
double rhonuc_cgs = rhonuc_si * 1e-3 ;
double c2_cgs = c_si * c_si * 1e4 ;
double m_neutron_MeV, m_proton_MeV ;
ifstream in_p_rho (file_thermo.data()) ;
in_p_rho >> m_neutron_MeV >> m_proton_MeV ; //Neutron and proton masses
in_p_rho.ignore(1000, '\n') ;
double p_convert = mev_si * 1.e45 * 10. ; // Conversion from MeV/fm^3 to cgs
double eps_convert = mev_si * 1.e42 / (c_si*c_si) ; //From meV/fm^3 to g/cm^3
// Main loop reading the table
//----------------------------
for (int i=0; i<nbp; i++) {
in_nb >> nb_fm3 ;
in_p_rho >> dummy_n >> dummy_n >> dummy_n >> p_over_nb_comp ;
in_p_rho >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> eps_comp ;
in_p_rho.ignore(1000, '\n') ;
p_cgs = p_over_nb_comp * nb_fm3 * p_convert ;
rho_cgs = ( eps_comp + 1. ) * m_neutron_MeV * nb_fm3 * eps_convert ;
press[i] = p_cgs / c2_cgs ;
nb[i] = nb_fm3 ;
ro[i] = rho_cgs ;
}
Tbl pp(nbp) ; pp.set_etat_qcq() ;
Tbl dh(nbp) ; dh.set_etat_qcq() ;
for (int i=0; i<nbp; i++) {
pp.set(i) = log(press[i] / rhonuc_cgs) ;
dh.set(i) = press[i] / (ro[i] + press[i]) ;
}
Tbl hh = integ1D(pp, dh) + 1.e-14 ;
for (int i=0; i<nbp; i++) {
logh->set(i) = log10( hh(i) ) ;
logp->set(i) = log10( press[i] / rhonuc_cgs ) ;
dlpsdlh->set(i) = hh(i) / dh(i) ;
lognb->set(i) = log10(nb[i]) ;
}
hmin = pow( double(10), (*logh)(0) ) ;
hmax = pow( double(10), (*logh)(nbp-1) ) ;
// Cleaning
//---------
delete [] press ;
delete [] nb ;
delete [] ro ;
}
//--------------//
// Destructor //
//--------------//
Eos_consistent::~Eos_consistent(){
// does nothing
}
//------------------------//
// Comparison operators //
//------------------------//
bool Eos_consistent::operator==(const Eos& eos_i) const {
bool resu = true ;
if ( eos_i.identify() != identify() ) {
cout << "The second EOS is not of type Eos_consistent !" << endl ;
resu = false ;
}
return resu ;
}
bool Eos_consistent::operator!=(const Eos& eos_i) const {
return !(operator==(eos_i)) ;
}
//------------------------------//
// Computational routines //
//------------------------------//
// Baryon density from enthalpy
//------------------------------
double Eos_consistent::nbar_ent_p(double ent, const Param* ) const {
static int i_near = logh->get_taille() / 2 ;
if ( ent > hmin ) {
if (ent > hmax) {
cout << "Eos_consistent::nbar_ent_p : ent > hmax !" << endl ;
abort() ;
}
double logh0 = log10( ent ) ;
double logp0 ;
double dlpsdlh0 ;
interpol_herm(*logh, *logp, *dlpsdlh, logh0, i_near, logp0,
dlpsdlh0) ;
double pp = pow(double(10), logp0) ;
double resu = pp / ent * dlpsdlh0 * exp(-ent) ;
if (i_near == 0)
{ // Use of linear interpolation for the first interval
double pp_near = pow(double(10), (*logp)(i_near)) ;
double ent_near = pow(double(10), (*logh)(i_near)) ;
resu = pp_near / ent_near * (*dlpsdlh)(i_near) * exp(-ent_near) ;
}
return resu ;
}
else{
return 0 ;
}
}
// Energy density from enthalpy
//------------------------------
double Eos_consistent::ener_ent_p(double ent, const Param* ) const {
static int i_near = logh->get_taille() / 2 ;
if ( ent > hmin ) {
if (ent > hmax) {
cout << "Eos_consistent::ener_ent_p : ent > hmax !" << endl ;
abort() ;
}
double logh0 = log10( ent ) ;
double logp0 ;
double dlpsdlh0 ;
interpol_herm(*logh, *logp, *dlpsdlh, logh0, i_near, logp0,
dlpsdlh0) ;
double pp = pow(double(10), logp0) ;
double resu = pp / ent * dlpsdlh0 - pp ;
if (i_near == 0)
{
double pp_near = pow(double(10), (*logp)(i_near)) ;
double ent_near = pow(double(10), (*logh)(i_near)) ;
resu = pp_near / ent_near * (*dlpsdlh)(i_near) - pp_near ;
}
return resu ;
}
else{
return 0 ;
}
}
// Pressure from enthalpy
//------------------------
double Eos_consistent::press_ent_p(double ent, const Param* ) const {
static int i_near = logh->get_taille() / 2 ;
if ( ent > hmin ) {
if (ent > hmax) {
cout << "Eos_consistent::press_ent_p : ent > hmax !" << endl ;
abort() ;
}
double logh0 = log10( ent ) ;
double logp0 ;
double dlpsdlh0 ;
interpol_herm(*logh, *logp, *dlpsdlh, logh0, i_near, logp0,
dlpsdlh0) ;
if (i_near == 0)
{
double logp_near = (*logp)(i_near) ;
double logp_nearp1 = (*logp)(i_near+1) ;
double delta = (*logh)(i_near+1) - (*logh)(i_near) ;
logp0 = (logp_nearp1*(logh0 - (*logh)(i_near))
- logp_near*(logh0 - (*logh)(i_near+1))) / delta ;
}
return pow(double(10), logp0) ;
}
else{
return 0 ;
}
}
//------------//
// Outputs //
//------------//
ostream& Eos_consistent::operator>>(ostream & ost) const {
ost << "EOS of class Eos_consistent." << endl ;
ost << "Built from file " << tablename << endl ;
ost << "Authors : " << authors << endl ;
ost << "Number of points in file : " << logh->get_taille() << endl ;
ost << "Table eventually slightly modified to ensure the relation" << endl ;
ost << "dp = (e+p) dh" << endl ;
return ost ;
}
}
Lorene::Eos_CompOSE::Eos_CompOSE
Eos_CompOSE(const string &files_path)
Constructor from CompOSE data.
Definition eos_compose.C:99
Lorene::Eos_consistent::operator==
virtual bool operator==(const Eos &) const
Comparison operator (egality).
Definition eos_consistent.C:182
Lorene::Eos_consistent::press_ent_p
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
Definition eos_consistent.C:276
Lorene::Eos_consistent::operator!=
virtual bool operator!=(const Eos &) const
Comparison operator (difference).
Definition eos_consistent.C:195
Lorene::Eos_consistent::operator>>
virtual ostream & operator>>(ostream &) const
Operator >>.
Definition eos_consistent.C:311
Lorene::Eos_consistent::Eos_consistent
Eos_consistent(const string &files_path)
Constructor from CompOSE data.
Definition eos_consistent.C:82
Lorene::Eos_consistent::nbar_ent_p
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
Definition eos_consistent.C:208
Lorene::Eos_consistent::identify
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
Definition eos_from_file.C:161
Lorene::Eos_consistent::~Eos_consistent
virtual ~Eos_consistent()
Destructor.
Definition eos_consistent.C:170
Lorene::Eos_consistent::ener_ent_p
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
Definition eos_consistent.C:242
Lorene::Eos_tabul::lognb
Tbl * lognb
Table of .
Definition eos_tabul.h:176
Lorene::Eos_tabul::dlpsdlh
Tbl * dlpsdlh
Table of .
Definition eos_tabul.h:173
Lorene::Eos_tabul::logh
Tbl * logh
Table of .
Definition eos_tabul.h:167
Lorene::Eos_tabul::tablename
string tablename
Name of the file containing the tabulated data.
Definition eos_tabul.h:156
Lorene::Eos_tabul::hmin
double hmin
Lower boundary of the enthalpy interval.
Definition eos_tabul.h:161
Lorene::Eos_tabul::logp
Tbl * logp
Table of .
Definition eos_tabul.h:170
Lorene::Eos_tabul::authors
string authors
Authors - reference for the table.
Definition eos_tabul.h:158
Lorene::Eos_tabul::hmax
double hmax
Upper boundary of the enthalpy interval.
Definition eos_tabul.h:164
Lorene::Eos
Equation of state base class.
Definition eos.h:190
Lorene::Eos::identify
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
Lorene::Param
Parameter storage.
Definition param.h:125
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::log10
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition cmp_math.C:322
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
Lorene::integ1D
Tbl integ1D(const Tbl &xx, const Tbl &ff)
Integrates a function defined on an unequally-spaced grid, approximating it by piece parabolae.
Definition integrate_1D.C:47
Lorene
Lorene prototypes.
Definition app_hor.h:64
Unites
Standard units of space, time and mass.