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 
30 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 $ " ;
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 
52 namespace 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 // ----------------------------
71 {}
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 
182 bool 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 
195 bool 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 
208 double 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 
242 double 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 
276 double 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 
311 ostream& 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
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
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.
Definition: integrate_1D.C:47
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.