LORENE
eos_consistent.C
1/*
2 * Methods of class Eos_consistent
3 *
4 * (see file eos_compose.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2015 Jerome Novak
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License as published by
15 * the Free Software Foundation; either version 2 of the License, or
16 * (at your option) any later version.
17 *
18 * LORENE is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with LORENE; if not, write to the Free Software
25 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 *
27 */
28
29
30char eos_consistent_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/eos_consistent.C,v 1.1 2015/08/04 14:41:29 j_novak Exp $ " ;
31
32/*
33 * $Id: eos_consistent.C,v 1.1 2015/08/04 14:41:29 j_novak Exp $
34 * $Log: eos_consistent.C,v $
35 * Revision 1.1 2015/08/04 14:41:29 j_novak
36 * Back to previous version for Eos_CompOSE. Enthalpy-consistent EoS can be accessed using Eos_consistent class (derived from Eos_CompOSE).
37 *
38 *
39 * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_consistent.C,v 1.1 2015/08/04 14:41:29 j_novak Exp $
40 *
41 */
42
43#include <string>
44
45// Headers Lorene
46#include "headcpp.h"
47#include "eos.h"
48#include "tbl.h"
49#include "utilitaires.h"
50#include "unites.h"
51
52namespace Lorene {
53 void interpol_herm(const Tbl& , const Tbl&, const Tbl&, double, int&,
54 double&, double& ) ;
55
56
57 //----------------------------//
58 // Constructors //
59 //----------------------------//
60
61// Standard constructor
62// --------------------
63 Eos_consistent::Eos_consistent(const char* file_name)
64 : Eos_CompOSE(file_name)
65{}
66
67
68// Constructor from binary file
69// ----------------------------
72
73
74// Constructor from a formatted file
75// ---------------------------------
77{}
78
79
80// Constructor from CompOSE data files
81// ------------------------------------
82 Eos_consistent::Eos_consistent(const string& files) : Eos_CompOSE(files)
83{
84 using namespace Unites ;
85
86 // Files containing data and a test
87 //---------------------------------
88 string file_nb = files + ".nb" ;
89 string file_thermo = files + ".thermo" ;
90
91 ifstream in_nb(file_nb.data()) ;
92
93 // obtaining the size of the tables for memory allocation
94 //-------------------------------------------------------
95 int index1, index2 ;
96 in_nb >> index1 >> index2 ;
97 int nbp = index2 - index1 + 1 ;
98 assert( nbp == logh->get_taille() ) ;
99
100 press = new double[nbp] ;
101 nb = new double[nbp] ;
102 ro = new double[nbp] ;
103
104 // Variables and conversion
105 //-------------------------
106 double nb_fm3, rho_cgs, p_cgs, p_over_nb_comp, eps_comp ;
107 double dummy_x ;
108 int dummy_n ;
109
110 double rhonuc_cgs = rhonuc_si * 1e-3 ;
111 double c2_cgs = c_si * c_si * 1e4 ;
112 double m_neutron_MeV, m_proton_MeV ;
113
114 ifstream in_p_rho (file_thermo.data()) ;
115 in_p_rho >> m_neutron_MeV >> m_proton_MeV ; //Neutron and proton masses
116 in_p_rho.ignore(1000, '\n') ;
117
118 double p_convert = mev_si * 1.e45 * 10. ; // Conversion from MeV/fm^3 to cgs
119 double eps_convert = mev_si * 1.e42 / (c_si*c_si) ; //From meV/fm^3 to g/cm^3
120
121 // Main loop reading the table
122 //----------------------------
123 for (int i=0; i<nbp; i++) {
124 in_nb >> nb_fm3 ;
125 in_p_rho >> dummy_n >> dummy_n >> dummy_n >> p_over_nb_comp ;
126 in_p_rho >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> eps_comp ;
127 in_p_rho.ignore(1000, '\n') ;
128 p_cgs = p_over_nb_comp * nb_fm3 * p_convert ;
129 rho_cgs = ( eps_comp + 1. ) * m_neutron_MeV * nb_fm3 * eps_convert ;
130
131 press[i] = p_cgs / c2_cgs ;
132 nb[i] = nb_fm3 ;
133 ro[i] = rho_cgs ;
134 }
135
136 Tbl pp(nbp) ; pp.set_etat_qcq() ;
137 Tbl dh(nbp) ; dh.set_etat_qcq() ;
138 for (int i=0; i<nbp; i++) {
139 pp.set(i) = log(press[i] / rhonuc_cgs) ;
140 dh.set(i) = press[i] / (ro[i] + press[i]) ;
141 }
142
143 Tbl hh = integ1D(pp, dh) + 1.e-14 ;
144
145 for (int i=0; i<nbp; i++) {
146 logh->set(i) = log10( hh(i) ) ;
147 logp->set(i) = log10( press[i] / rhonuc_cgs ) ;
148 dlpsdlh->set(i) = hh(i) / dh(i) ;
149 lognb->set(i) = log10(nb[i]) ;
150 }
151
152 hmin = pow( double(10), (*logh)(0) ) ;
153 hmax = pow( double(10), (*logh)(nbp-1) ) ;
154
155 // Cleaning
156 //---------
157 delete [] press ;
158 delete [] nb ;
159 delete [] ro ;
160
161
162}
163
164
165
166 //--------------//
167 // Destructor //
168 //--------------//
169
171
172 // does nothing
173
174}
175
176
177 //------------------------//
178 // Comparison operators //
179 //------------------------//
180
181
182bool Eos_consistent::operator==(const Eos& eos_i) const {
183
184 bool resu = true ;
185
186 if ( eos_i.identify() != identify() ) {
187 cout << "The second EOS is not of type Eos_consistent !" << endl ;
188 resu = false ;
189 }
190
191 return resu ;
192
193}
194
195bool Eos_consistent::operator!=(const Eos& eos_i) const {
196
197 return !(operator==(eos_i)) ;
198
199}
200
201 //------------------------------//
202 // Computational routines //
203 //------------------------------//
204
205// Baryon density from enthalpy
206//------------------------------
207
208double Eos_consistent::nbar_ent_p(double ent, const Param* ) const {
209
210 static int i_near = logh->get_taille() / 2 ;
211
212 if ( ent > hmin ) {
213 if (ent > hmax) {
214 cout << "Eos_consistent::nbar_ent_p : ent > hmax !" << endl ;
215 abort() ;
216 }
217 double logh0 = log10( ent ) ;
218 double logp0 ;
219 double dlpsdlh0 ;
220 interpol_herm(*logh, *logp, *dlpsdlh, logh0, i_near, logp0,
221 dlpsdlh0) ;
222
223 double pp = pow(double(10), logp0) ;
224
225 double resu = pp / ent * dlpsdlh0 * exp(-ent) ;
226 if (i_near == 0)
227 { // Use of linear interpolation for the first interval
228 double pp_near = pow(double(10), (*logp)(i_near)) ;
229 double ent_near = pow(double(10), (*logh)(i_near)) ;
230 resu = pp_near / ent_near * (*dlpsdlh)(i_near) * exp(-ent_near) ;
231 }
232 return resu ;
233 }
234 else{
235 return 0 ;
236 }
237}
238
239// Energy density from enthalpy
240//------------------------------
241
242double Eos_consistent::ener_ent_p(double ent, const Param* ) const {
243
244 static int i_near = logh->get_taille() / 2 ;
245
246 if ( ent > hmin ) {
247 if (ent > hmax) {
248 cout << "Eos_consistent::ener_ent_p : ent > hmax !" << endl ;
249 abort() ;
250 }
251 double logh0 = log10( ent ) ;
252 double logp0 ;
253 double dlpsdlh0 ;
254 interpol_herm(*logh, *logp, *dlpsdlh, logh0, i_near, logp0,
255 dlpsdlh0) ;
256
257 double pp = pow(double(10), logp0) ;
258
259 double resu = pp / ent * dlpsdlh0 - pp ;
260 if (i_near == 0)
261 {
262 double pp_near = pow(double(10), (*logp)(i_near)) ;
263 double ent_near = pow(double(10), (*logh)(i_near)) ;
264 resu = pp_near / ent_near * (*dlpsdlh)(i_near) - pp_near ;
265 }
266 return resu ;
267 }
268 else{
269 return 0 ;
270 }
271}
272
273// Pressure from enthalpy
274//------------------------
275
276double Eos_consistent::press_ent_p(double ent, const Param* ) const {
277
278 static int i_near = logh->get_taille() / 2 ;
279
280 if ( ent > hmin ) {
281 if (ent > hmax) {
282 cout << "Eos_consistent::press_ent_p : ent > hmax !" << endl ;
283 abort() ;
284 }
285
286 double logh0 = log10( ent ) ;
287 double logp0 ;
288 double dlpsdlh0 ;
289 interpol_herm(*logh, *logp, *dlpsdlh, logh0, i_near, logp0,
290 dlpsdlh0) ;
291 if (i_near == 0)
292 {
293 double logp_near = (*logp)(i_near) ;
294 double logp_nearp1 = (*logp)(i_near+1) ;
295 double delta = (*logh)(i_near+1) - (*logh)(i_near) ;
296 logp0 = (logp_nearp1*(logh0 - (*logh)(i_near))
297 - logp_near*(logh0 - (*logh)(i_near+1))) / delta ;
298 }
299 return pow(double(10), logp0) ;
300 }
301 else{
302 return 0 ;
303 }
304}
305
306 //------------//
307 // Outputs //
308 //------------//
309
310
311ostream& Eos_consistent::operator>>(ostream & ost) const {
312
313 ost << "EOS of class Eos_consistent." << endl ;
314 ost << "Built from file " << tablename << endl ;
315 ost << "Authors : " << authors << endl ;
316 ost << "Number of points in file : " << logh->get_taille() << endl ;
317 ost << "Table eventually slightly modified to ensure the relation" << endl ;
318 ost << "dp = (e+p) dh" << endl ;
319 return ost ;
320}
321
322
323}
Equation of state for the CompOSE database.
Definition eos_compose.h:77
virtual bool operator==(const Eos &) const
Comparison operator (egality)
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
virtual ostream & operator>>(ostream &) const
Operator >>
Eos_consistent(const string &files_path)
Constructor from CompOSE data.
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual ~Eos_consistent()
Destructor.
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
Tbl * lognb
Table of .
Definition eos_tabul.h:176
Tbl * dlpsdlh
Table of .
Definition eos_tabul.h:173
Tbl * logh
Table of .
Definition eos_tabul.h:167
string tablename
Name of the file containing the tabulated data.
Definition eos_tabul.h:156
double hmin
Lower boundary of the enthalpy interval.
Definition eos_tabul.h:161
Tbl * logp
Table of .
Definition eos_tabul.h:170
string authors
Authors - reference for the table.
Definition eos_tabul.h:158
double hmax
Upper boundary of the enthalpy interval.
Definition eos_tabul.h:164
Equation of state base class.
Definition eos.h:190
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
Parameter storage.
Definition param.h:125
Basic array class.
Definition tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
int get_taille() const
Gives the total size (ie dim.taille)
Definition tbl.h:397
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition cmp_math.C:322
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
Tbl integ1D(const Tbl &xx, const Tbl &ff)
Integrates a function defined on an unequally-spaced grid, approximating it by piece parabolae.
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.