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 
31 
32 /*
33  * $Id: eos_consistent.C,v 1.7 2022/03/10 16:38:39 j_novak Exp $
34  * $Log: eos_consistent.C,v $
35  * Revision 1.7 2022/03/10 16:38:39 j_novak
36  * log(cs^2) is tabulated instead of cs^2.
37  *
38  * Revision 1.6 2021/05/31 11:31:23 g_servignat
39  * Added csound_square_ent routine to calculate the sound speed from enthalpy to Eos_consistent and corrected error outputs
40  *
41  * Revision 1.5 2019/03/28 13:41:02 j_novak
42  * Improved managed of saved EoS (functions sauve and constructor form FILE*)
43  *
44  * Revision 1.4 2017/08/11 13:42:04 j_novak
45  * Suppression of spurious output
46  *
47  * Revision 1.3 2017/08/10 15:14:27 j_novak
48  * Now Eos_consistent is also reading Lorene format tables.
49  *
50  * Revision 1.2 2016/12/05 16:17:51 j_novak
51  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
52  *
53  * Revision 1.1 2015/08/04 14:41:29 j_novak
54  * Back to previous version for Eos_CompOSE. Enthalpy-consistent EoS can be accessed using Eos_consistent class (derived from Eos_CompOSE).
55  *
56  *
57  * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_consistent.C,v 1.7 2022/03/10 16:38:39 j_novak Exp $
58  *
59  */
60 
61 #include <string>
62 
63 // Headers Lorene
64 #include "headcpp.h"
65 #include "eos.h"
66 #include "tbl.h"
67 #include "utilitaires.h"
68 #include "unites.h"
69 #include<string>
70 
71 namespace Lorene {
72  void interpol_herm(const Tbl& , const Tbl&, const Tbl&, double, int&,
73  double&, double& ) ;
74  void interpol_linear(const Tbl&, const Tbl&, double, int&, double&) ;
75 
76  //----------------------------//
77  // Constructors //
78  //----------------------------//
79 
80 // Standard constructor
81 // --------------------
82  Eos_consistent::Eos_consistent(const char* file_name)
83  : Eos_CompOSE(file_name) { }
84 
85 
86 // Constructor from binary file
87 // ----------------------------
89  if (format == 0) read_table() ;
90  else read_compose_data() ;
91  }
92 
93 
94 // Constructor from a formatted file
95 // ---------------------------------
96  Eos_consistent::Eos_consistent(ifstream& fich) : Eos_CompOSE(fich) {
97  read_table() ;
98  }
99 
100 
101 // Constructor from CompOSE data files
102 // ------------------------------------
103  Eos_consistent::Eos_consistent(const string& files) : Eos_CompOSE(files)
104  {
106  }
107 
109 
110  using namespace Unites ;
111 
112  cout << "Computing a thermodynamically - consistent table." << endl ;
113 
114  ifstream tfich(tablename.data()) ;
115 
116  for (int i=0; i<5; i++) { // jump over the file
117  tfich.ignore(1000, '\n') ; // header
118  } //
119 
120  int nbp ;
121  tfich >> nbp ; tfich.ignore(1000, '\n') ; // number of data
122  if (nbp<=0) {
123  cerr << "Eos_consistent::read_table(): " << endl ;
124  cerr << "Wrong value for the number of lines!" << endl ;
125  cerr << "nbp = " << nbp << endl ;
126  cerr << "Aborting..." << endl ;
127  abort() ;
128  }
129 
130  for (int i=0; i<3; i++) { // jump over the table
131  tfich.ignore(1000, '\n') ; // description
132  }
133 
134  press = new double[nbp] ;
135  nb = new double[nbp] ;
136  ro = new double[nbp] ;
137 
138  double rhonuc_cgs = rhonuc_si * 1e-3 ;
139  double c2_cgs = c_si * c_si * 1e4 ;
140 
141  int no ;
142  double nb_fm3, rho_cgs, p_cgs ;
143 
144  for (int i=0; i<nbp; i++) {
145 
146  tfich >> no ;
147  tfich >> nb_fm3 ;
148  tfich >> rho_cgs ;
149  tfich >> p_cgs ; tfich.ignore(1000,'\n') ;
150  if ( (nb_fm3<0) || (rho_cgs<0) || (p_cgs < 0) ){
151  cout << "Eos_consistent(ifstream&): " << endl ;
152  cout << "Negative value in table!" << endl ;
153  cout << "nb = " << nb_fm3 << ", rho = " << rho_cgs <<
154  ", p = " << p_cgs << endl ;
155  cout << "Aborting..." << endl ;
156  abort() ;
157  }
158 
159  press[i] = p_cgs / c2_cgs ;
160  nb[i] = nb_fm3 ;
161  ro[i] = rho_cgs ;
162  }
163 
164  Tbl pp(nbp) ; pp.set_etat_qcq() ;
165  Tbl dh(nbp) ; dh.set_etat_qcq() ;
166  for (int i=0; i<nbp; i++) {
167  pp.set(i) = log(press[i] / rhonuc_cgs) ;
168  dh.set(i) = press[i] / (ro[i] + press[i]) ;
169  }
170 
171  Tbl hh = integ1D(pp, dh) + 1.e-14 ;
172 
173  for (int i=0; i<nbp; i++) {
174  logh->set(i) = log10( hh(i) ) ;
175  logp->set(i) = log10( press[i] / rhonuc_cgs ) ;
176  dlpsdlh->set(i) = hh(i) / dh(i) ;
177  lognb->set(i) = log10(nb[i]) ;
178  }
179 
180  hmin = pow( double(10), (*logh)(0) ) ;
181  hmax = pow( double(10), (*logh)(nbp-1) ) ;
182 
183  // Cleaning
184  //---------
185  delete [] press ;
186  delete [] nb ;
187  delete [] ro ;
188 
189  }
190 
191 
193 
194  using namespace Unites ;
195 
196  cout << "Computing a thermodynamically - consistent table." << endl ;
197 
198  // Files containing data and a test
199  //---------------------------------
200  string file_nb = tablename + ".nb" ;
201  string file_thermo = tablename + ".thermo" ;
202 
203  ifstream in_nb(file_nb.data()) ;
204 
205  // obtaining the size of the tables for memory allocation
206  //-------------------------------------------------------
207  int index1, index2 ;
208  in_nb >> index1 >> index2 ;
209  int nbp = index2 - index1 + 1 ;
210  assert( nbp == logh->get_taille() ) ;
211 
212  press = new double[nbp] ;
213  nb = new double[nbp] ;
214  ro = new double[nbp] ;
215 
216  // Variables and conversion
217  //-------------------------
218  double nb_fm3, rho_cgs, p_cgs, p_over_nb_comp, eps_comp ;
219  double dummy_x ;
220  int dummy_n ;
221 
222  double rhonuc_cgs = rhonuc_si * 1e-3 ;
223  double c2_cgs = c_si * c_si * 1e4 ;
224  double m_neutron_MeV, m_proton_MeV ;
225 
226  ifstream in_p_rho (file_thermo.data()) ;
227  in_p_rho >> m_neutron_MeV >> m_proton_MeV ; //Neutron and proton masses
228  in_p_rho.ignore(1000, '\n') ;
229 
230  double p_convert = mev_si * 1.e45 * 10. ; // Conversion from MeV/fm^3 to cgs
231  double eps_convert = mev_si * 1.e42 / (c_si*c_si) ; //From meV/fm^3 to g/cm^3
232 
233  // Main loop reading the table
234  //----------------------------
235  for (int i=0; i<nbp; i++) {
236  in_nb >> nb_fm3 ;
237  in_p_rho >> dummy_n >> dummy_n >> dummy_n >> p_over_nb_comp ;
238  in_p_rho >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> eps_comp ;
239  in_p_rho.ignore(1000, '\n') ;
240  p_cgs = p_over_nb_comp * nb_fm3 * p_convert ;
241  rho_cgs = ( eps_comp + 1. ) * m_neutron_MeV * nb_fm3 * eps_convert ;
242 
243  press[i] = p_cgs / c2_cgs ;
244  nb[i] = nb_fm3 ;
245  ro[i] = rho_cgs ;
246  }
247 
248  Tbl pp(nbp) ; pp.set_etat_qcq() ;
249  Tbl dh(nbp) ; dh.set_etat_qcq() ;
250  for (int i=0; i<nbp; i++) {
251  pp.set(i) = log(press[i] / rhonuc_cgs) ;
252  dh.set(i) = press[i] / (ro[i] + press[i]) ;
253  }
254 
255  Tbl hh = integ1D(pp, dh) + 1.e-14 ;
256 
257  for (int i=0; i<nbp; i++) {
258  logh->set(i) = log10( hh(i) ) ;
259  logp->set(i) = log10( press[i] / rhonuc_cgs ) ;
260  dlpsdlh->set(i) = hh(i) / dh(i) ;
261  lognb->set(i) = log10(nb[i]) ;
262  }
263 
264  hmin = pow( double(10), (*logh)(0) ) ;
265  hmax = pow( double(10), (*logh)(nbp-1) ) ;
266 
267  // Cleaning
268  //---------
269  delete [] press ;
270  delete [] nb ;
271  delete [] ro ;
272 
273  }
274 
275 
276  //--------------//
277  // Destructor //
278  //--------------//
279 
281 
282  // does nothing
283 
284 }
285 
286 
287  //------------------------//
288  // Comparison operators //
289  //------------------------//
290 
291 
292 bool Eos_consistent::operator==(const Eos& eos_i) const {
293 
294  bool resu = true ;
295 
296  if ( eos_i.identify() != identify() ) {
297  cout << "The second EOS is not of type Eos_consistent !" << endl ;
298  resu = false ;
299  }
300 
301  return resu ;
302 
303 }
304 
305 bool Eos_consistent::operator!=(const Eos& eos_i) const {
306 
307  return !(operator==(eos_i)) ;
308 
309 }
310 
311  //------------------------------//
312  // Computational routines //
313  //------------------------------//
314 
315 // Baryon density from enthalpy
316 //------------------------------
317 
318 double Eos_consistent::nbar_ent_p(double ent, const Param* ) const {
319 
320  static int i_near = logh->get_taille() / 2 ;
321 
322  if ( ent > hmin ) {
323  if (ent > hmax) {
324  cout << "Eos_consistent::nbar_ent_p : ent > hmax !" << endl ;
325  abort() ;
326  }
327  double logh0 = log10( ent ) ;
328  double logp0 ;
329  double dlpsdlh0 ;
330  interpol_herm(*logh, *logp, *dlpsdlh, logh0, i_near, logp0,
331  dlpsdlh0) ;
332 
333  double pp = pow(double(10), logp0) ;
334 
335  double resu = pp / ent * dlpsdlh0 * exp(-ent) ;
336  if (i_near == 0)
337  { // Use of linear interpolation for the first interval
338  double pp_near = pow(double(10), (*logp)(i_near)) ;
339  double ent_near = pow(double(10), (*logh)(i_near)) ;
340  resu = pp_near / ent_near * (*dlpsdlh)(i_near) * exp(-ent_near) ;
341  }
342  return resu ;
343  }
344  else{
345  return 0 ;
346  }
347 }
348 
349 // Energy density from enthalpy
350 //------------------------------
351 
352 double Eos_consistent::ener_ent_p(double ent, const Param* ) const {
353 
354  static int i_near = logh->get_taille() / 2 ;
355 
356  if ( ent > hmin ) {
357  if (ent > hmax) {
358  cout << "Eos_consistent::ener_ent_p : ent > hmax !" << endl ;
359  abort() ;
360  }
361  double logh0 = log10( ent ) ;
362  double logp0 ;
363  double dlpsdlh0 ;
364  interpol_herm(*logh, *logp, *dlpsdlh, logh0, i_near, logp0,
365  dlpsdlh0) ;
366 
367  double pp = pow(double(10), logp0) ;
368 
369  double resu = pp / ent * dlpsdlh0 - pp ;
370  if (i_near == 0)
371  {
372  double p_near = pow(double(10), (*logp)(i_near)) ;
373  double p_nearp1 = pow(double(10), (*logp)(i_near+1)) ;
374  double ener_near = p_near/ pow(double(10), (*logh)(i_near))
375  * (*dlpsdlh)(i_near) - p_near ;
376  double ener_nearp1 = p_nearp1/ pow(double(10), (*logh)(i_near+1))
377  * (*dlpsdlh)(i_near+1) - p_nearp1 ;
378  double delta = (*logh)(i_near+1) - (*logh)(i_near) ;
379  resu = (ener_nearp1*(logh0 - (*logh)(i_near))
380  - ener_near*(logh0 - (*logh)(i_near+1))) / delta ;
381  }
382  return resu ;
383  }
384  else{
385  return 0 ;
386  }
387 }
388 
389 // Pressure from enthalpy
390 //------------------------
391 
392 double Eos_consistent::press_ent_p(double ent, const Param* ) const {
393 
394  static int i_near = logh->get_taille() / 2 ;
395 
396  if ( ent > hmin ) {
397  if (ent > hmax) {
398  cout << "Eos_consistent::press_ent_p : ent > hmax !" << endl ;
399  abort() ;
400  }
401 
402  double logh0 = log10( ent ) ;
403  double logp0 ;
404  double dlpsdlh0 ;
405  interpol_herm(*logh, *logp, *dlpsdlh, logh0, i_near, logp0,
406  dlpsdlh0) ;
407  if (i_near == 0)
408  {
409  double logp_near = (*logp)(i_near) ;
410  double logp_nearp1 = (*logp)(i_near+1) ;
411  double delta = (*logh)(i_near+1) - (*logh)(i_near) ;
412  logp0 = (logp_nearp1*(logh0 - (*logh)(i_near))
413  - logp_near*(logh0 - (*logh)(i_near+1))) / delta ;
414  }
415  return pow(double(10), logp0) ;
416  }
417  else{
418  return 0 ;
419  }
420 }
421 
422 // Square of sound speed from enthalpy
423 double Eos_consistent::csound_square_ent_p(double ent, const Param*) const {
424  static int i_near = lognb->get_taille() / 2 ;
425 
426  if ( ent > hmin ) {
427  if (ent > hmax) {
428  cout << "Eos_consistent::csound_square_ent_p : ent>hmax !" << endl ;
429  abort() ;
430  }
431  double log_ent0 = log10( ent ) ;
432  double log_csound0 ;
433 
434  interpol_linear(*logh, *log_cs2, log_ent0, i_near, log_csound0) ;
435 
436  return pow(10., log_csound0) ;
437  }
438  else
439  {
440  return pow(10., (*log_cs2)(0)) ;
441  }
442 }
443 
444  //------------//
445  // Outputs //
446  //------------//
447 
448 
449 ostream& Eos_consistent::operator>>(ostream & ost) const {
450 
451  ost << "EOS of class Eos_consistent." << endl ;
452  ost << "Built from file " << tablename << endl ;
453  ost << "Authors : " << authors << endl ;
454  ost << "Number of points in file : " << logh->get_taille() << endl ;
455  ost << "Table eventually slightly modified to ensure the relation" << endl ;
456  ost << "dp = (e+p) dh" << endl ;
457  return ost ;
458 }
459 
460 
461 }
virtual ostream & operator>>(ostream &) const
Operator >>
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
virtual ~Eos_consistent()
Destructor.
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
Equation of state base class.
Definition: eos.h:209
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual void read_compose_data()
Reads the files containing the table and initializes in the arrays logh , logp and dlpsdlh (CompOSE f...
Tbl * logp
Table of .
Definition: eos_tabul.h:206
Tbl * log_cs2
Table of .
Definition: eos_tabul.h:218
double hmin
Lower boundary of the enthalpy interval.
Definition: eos_tabul.h:197
virtual void read_table()
Reads the file containing the table and initializes in the arrays logh , logp and dlpsdlh ...
Tbl * dlpsdlh
Table of .
Definition: eos_tabul.h:209
int format
0 for standard (old) LORENE format, 1 for CompOSE format
Definition: eos_compose.h:99
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
virtual double csound_square_ent_p(double, const Param *) const
Computes the sound speed squared from the enthapy with extra parameters (virtual function implemente...
string tablename
Name of the file containing the tabulated data.
Definition: eos_tabul.h:192
string authors
Authors - reference for the table.
Definition: eos_tabul.h:194
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
Tbl * lognb
Table of .
Definition: eos_tabul.h:212
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
double hmax
Upper boundary of the enthalpy interval.
Definition: eos_tabul.h:200
Parameter storage.
Definition: param.h:125
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
Tbl integ1D(const Tbl &xx, const Tbl &ff)
Integrates a function defined on an unequally-spaced grid, approximating it by piecewise parabolae...
Definition: integrate_1D.C:50
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:325
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:417
Basic array class.
Definition: tbl.h:164
virtual bool operator==(const Eos &) const
Comparison operator (egality)
Tbl * logh
Table of .
Definition: eos_tabul.h:203
Eos_consistent(const string &files_path)
Constructor from CompOSE data.
Equation of state for the CompOSE database.
Definition: eos_compose.h:93