LORENE
dyneos_cons.C
1 /*
2  * Methods of class Dyn_eos_cons
3  *
4  * (see file dyneos.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2019 Jerome Novak
10  * (c) 2000 Eric Gourgoulhon for Eos classes
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 
32 /*
33  * $Id: dyneos_cons.C,v 1.3 2022/03/22 13:36:00 j_novak Exp $
34  * $Log: dyneos_cons.C,v $
35  * Revision 1.3 2022/03/22 13:36:00 j_novak
36  * Added declaration of compute_derivative to utilitaires.h
37  *
38  * Revision 1.2 2020/12/17 17:00:27 j_novak
39  * Output of sound speed squared, instead of sound speed.
40  *
41  * Revision 1.1 2019/12/06 14:30:50 j_novak
42  * New classes Dyn_eos... for cold Eos's with baryon density as input.
43  *
44  *
45  *
46  * $Header: /cvsroot/Lorene/C++/Source/Eos/dyneos_cons.C,v 1.3 2022/03/22 13:36:00 j_novak Exp $
47  *
48  */
49 
50 // Headers Lorene
51 #include "dyneos.h"
52 #include "tbl.h"
53 #include "utilitaires.h"
54 #include "unites.h"
55 
56 
57 namespace Lorene {
58 
59  //----------------------------//
60  // Constructors //
61  //----------------------------//
62 
63 // Standard constructor
64 // --------------------
65  Dyn_eos_cons::Dyn_eos_cons(const string& name_i, const string& tablename_i,
66  bool compose) : Dyn_eos_tab()
67  {
68  name = name_i ;
69  tablename = tablename_i ;
70  compose_format = compose ;
71  if (compose_format)
73  else
75  }
76 
77 // Constructor from binary file
78 // ----------------------------
80  {
81  const int nc1 = 100 ;
82  const int nc2 = 160 ;
83  char tmp_s1[nc1] ;
84  size_t ret = fread(tmp_s1, sizeof(char), nc1, fich) ;
85  if (int(ret) == nc1)
86  name = tmp_s1 ;
87  char tmp_s2[nc2] ;
88  ret = fread(tmp_s2, sizeof(char), nc2, fich) ;
89  if (ret == nc2)
90  tablename = tmp_s2 ;
91  int comp ;
92  fread_be(&comp, sizeof(int), 1, fich) ;
93  compose_format = comp ;
94 
95  if (compose_format)
97  else
99  }
100 
101 
102 // Constructor from a formatted file
103 // ---------------------------------
105  {
106  fich.seekg(0, fich.beg) ;
107  fich.ignore(1000, '\n') ;
108  fich >> compose_format ;
109  fich.ignore(1000, '\n') ;
110  getline(fich, name, '\n') ;
111  fich >> tablename ;
112 
113  if (compose_format)
115  else
117  }
118 
119  //--------------//
120  // Destructor //
121  //--------------//
122 
124  {
125  }
126 
127  //------------------------//
128  // Comparison operators //
129  //------------------------//
130 
131 
132  bool Dyn_eos_cons::operator==(const Dyn_eos& eos_i) const {
133 
134  bool resu = true ;
135 
136  if ( eos_i.identify() != identify() ) {
137  cout << "The second EOS is not of type Dyn_eos_cons !" << endl ;
138  resu = false ;
139  }
140 
141  return resu ;
142 
143  }
144 
145  bool Dyn_eos_cons::operator!=(const Dyn_eos& eos_i) const {
146 
147  return !(operator==(eos_i)) ;
148 
149  }
150 
151  //------------//
152  // Outputs //
153  //------------//
154 
155  ostream& Dyn_eos_cons::operator>>(ostream & ost) const
156  {
157  ost << "EOS of class Dyn_eos_cons." << endl ;
158  ost << "Built from file " << tablename << endl ;
159  ost << ((compose_format == 0) ? "LORENE format" : "CompOSE format") << endl ;
160  ost << "Authors : " << authors << endl ;
161  ost << "Number of points in file : " << lognb->get_taille() << endl ;
162  ost << "Table eventually slightly modified to ensure the relation" << endl ;
163  ost << "de/dn = (e+p)/n" << endl ;
164  return ost ;
165  }
166  //------------------------//
167  // Reading of the table //
168  //------------------------//
169 
170  // LORENE format
171  //---------------
173 
174  using namespace Unites ;
175 
176  ifstream fich(tablename.data()) ;
177 
178  if (!fich) {
179  cerr << "Dyn_eos_cons::read_table_lorene(): " << endl ;
180  cerr << "Problem in opening the EOS file!" << endl ;
181  cerr << "While trying to open " << tablename << endl ;
182  cerr << "Aborting..." << endl ;
183  abort() ;
184  }
185 
186  fich.ignore(1000, '\n') ; // Jump over the first header
187  fich.ignore(1) ;
188  getline(fich, authors, '\n') ;
189  for (int i=0; i<3; i++) { // jump over the file
190  fich.ignore(1000, '\n') ; // header
191  } //
192 
193  int nbp ;
194  fich >> nbp ; fich.ignore(1000, '\n') ; // number of data
195  if (nbp<=0) {
196  cerr << "Dyn_eos_cons::read_table_lorene(): " << endl ;
197  cerr << "Wrong value for the number of lines!" << endl ;
198  cerr << "nbp = " << nbp << endl ;
199  cerr << "Aborting..." << endl ;
200  abort() ;
201  }
202 
203  for (int i=0; i<3; i++) { // jump over the table
204  fich.ignore(1000, '\n') ; // description
205  }
206 
207  Tbl press(nbp) ; press.set_etat_qcq() ;
208  Tbl nb(nbp) ; nb.set_etat_qcq() ;
209  Tbl ro(nbp) ; ro.set_etat_qcq() ;
210 
211  lognb = new Tbl(nbp) ;
212  loge = new Tbl(nbp) ;
213  dlesdlnb = new Tbl(nbp) ;
214 
215  lognb->set_etat_qcq() ;
216  loge->set_etat_qcq() ;
218 
219  double rhonuc_cgs = rhonuc_si * 1e-3 ;
220  double c2_cgs = c_si * c_si * 1e4 ;
221 
222  int no ;
223  double nb_fm3, rho_cgs, p_cgs ;
224 
225  cout << "Dyn_eos_cons: reading Lorene format table from the file "
226  << tablename << endl ;
227  for (int i=0; i<nbp; i++) {
228 
229  fich >> no ;
230  fich >> nb_fm3 ;
231  fich >> rho_cgs ;
232  fich >> p_cgs ; fich.ignore(1000,'\n') ;
233  if ( (nb_fm3<0) || (rho_cgs<0) || (p_cgs < 0) ){
234  cout << "Dyn_eos_cons::read_table_lorene(): " << endl ;
235  cout << "Negative value in table!" << endl ;
236  cout << "nb = " << nb_fm3 << ", rho = " << rho_cgs <<
237  ", p = " << p_cgs << endl ;
238  cout << "Aborting..." << endl ;
239  abort() ;
240  }
241 
242  press.set(i) = p_cgs / (c2_cgs*rhonuc_cgs) ;
243  nb.set(i) = 10.*nb_fm3 ; // Units 0.1 fm^-3
244  ro.set(i) = rho_cgs / rhonuc_cgs ;
245  }
246 
247  Tbl ee(nbp) ; ee.set_etat_qcq() ;
248  Tbl dn(nbp) ; dn.set_etat_qcq() ;
249  ee = log(ro) ;
250  dn = ro / ( ro + press ) ;
251 
252  Tbl nb2 = exp( integ1D(ee, dn) + log(nb(0)) ) ;
253  double maxdif = diffrelmax(nb2, nb) ;
254  cout << "Dyn_eos_cons: max. relative difference between new (consistent)" << endl ;
255  cout << "and old density data: " << maxdif << endl ;
256 
257  *lognb = log10(nb2) ;
258  *loge = log10(ro) ;
259  *dlesdlnb = (ro + press) / ro ;
260  Tbl tmp(nbp) ; tmp.set_etat_qcq() ;
261  compute_derivative(ro, press, tmp) ;
262  c_sound2 = new Tbl(tmp) ; // c_s^2 = dp/de
263 
264  nbmin = pow( double(10), (*lognb)(0) ) ;
265  nbmax = pow( double(10), (*lognb)(nbp-1) ) ;
266 
267  cout << "nbmin, nbmax : " << 0.1*nbmin << " " << 0.1*nbmax << " fm^-3" << endl ;
268 
269  fich.close();
270  }
271 
272  // CompOSE format
273  //----------------
275  {
276  using namespace Unites ;
277 
278  // Files containing data and a test
279  //---------------------------------
280  string file_nb = tablename + ".nb" ;
281  string file_thermo = tablename + ".thermo" ;
282 
283  ifstream in_nb(file_nb.data()) ;
284  if (!in_nb) {
285  cerr << "Dyn_eos_cons::read_table_compose(): " << endl ;
286  cerr << "Problem in opening the EOS file!" << endl ;
287  cerr << "While trying to open " << file_nb << endl ;
288  cerr << "Aborting..." << endl ;
289  abort() ;
290  }
291 
292  // obtaining the size of the tables for memory allocation
293  //-------------------------------------------------------
294  int index1, index2 ;
295  in_nb >> index1 >> index2 ;
296  int nbp = index2 - index1 + 1 ;
297  assert(nbp > 0) ;
298 
299  Tbl press(nbp) ; press.set_etat_qcq() ;
300  Tbl nb(nbp) ; nb.set_etat_qcq() ;
301  Tbl ro(nbp) ; ro.set_etat_qcq() ;
302 
303  lognb = new Tbl(nbp) ;
304  loge = new Tbl(nbp) ;
305  dlesdlnb = new Tbl(nbp) ;
306 
307  lognb->set_etat_qcq() ;
308  loge->set_etat_qcq() ;
310 
311  // Variables and conversion
312  //-------------------------
313  double nb_fm3, rho_cgs, p_cgs, p_over_nb_comp, eps_comp ;
314  double dummy_x ;
315  int dummy_n ;
316 
317  double rhonuc_cgs = rhonuc_si * 1e-3 ;
318  double c2_cgs = c_si * c_si * 1e4 ;
319  double m_neutron_MeV, m_proton_MeV ;
320 
321  ifstream in_p_rho (file_thermo.data()) ;
322  if (!in_p_rho) {
323  cerr << "Dyn_eos_cons::read_table_compose(): " << endl ;
324  cerr << "Problem in opening the EOS file!" << endl ;
325  cerr << "While trying to open " << file_thermo << endl ;
326  cerr << "Aborting..." << endl ;
327  abort() ;
328  }
329  in_p_rho >> m_neutron_MeV >> m_proton_MeV ; //Neutron and proton masses
330  in_p_rho.ignore(1000, '\n') ;
331 
332  double p_convert = mev_si * 1.e45 * 10. ; // Conversion from MeV/fm^3 to cgs
333  double eps_convert = mev_si * 1.e42 / (c_si*c_si) ; //From meV/fm^3 to g/cm^3
334 
335  // Main loop reading the table
336  //----------------------------
337  cout << "Dyn_eos_cons: reading CompOSE format table from the file "
338  << tablename << ".thermo" << endl ;
339  for (int i=0; i<nbp; i++) {
340  in_nb >> nb_fm3 ;
341  in_p_rho >> dummy_n >> dummy_n >> dummy_n >> p_over_nb_comp ;
342  in_p_rho >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> eps_comp ;
343  in_p_rho.ignore(1000, '\n') ;
344  p_cgs = p_over_nb_comp * nb_fm3 * p_convert ;
345  rho_cgs = ( eps_comp + 1. ) * m_neutron_MeV * nb_fm3 * eps_convert ;
346 
347  if ( (nb_fm3<0) || (rho_cgs<0) || (p_cgs < 0) ){
348  cout << "Dyn_eos_cons::read_table_compose(): " << endl ;
349  cout << "Negative value in table!" << endl ;
350  cout << "nb = " << nb_fm3 << ", rho = " << rho_cgs <<
351  ", p = " << p_cgs << endl ;
352  cout << "Aborting..." << endl ;
353  abort() ;
354  }
355 
356  press.set(i) = (p_cgs / c2_cgs) / rhonuc_cgs ;
357  nb.set(i) = 10. * nb_fm3 ; // Units : 0.1 fm^-3
358  ro.set(i) = rho_cgs / rhonuc_cgs ;
359  }
360 
361  Tbl ee(nbp) ; ee.set_etat_qcq() ;
362  Tbl dn(nbp) ; dn.set_etat_qcq() ;
363  ee = log(ro) ;
364  dn = ro / ( ro + press ) ;
365 
366  Tbl nb2 = exp( integ1D(ee, dn) + log(nb(0)) ) ;
367  double maxdif = diffrelmax(nb2, nb) ;
368  cout << "Dyn_eos_cons: max. relative difference between new (consistent)" << endl ;
369  cout << "and old density data: " << maxdif << endl ;
370 
371  *lognb = log10(nb2) ;
372  *loge = log10(ro) ;
373  *dlesdlnb = (ro + press) / ro ;
374  Tbl tmp(nbp) ; tmp.set_etat_qcq() ;
375  compute_derivative(ro, press, tmp) ;
376  c_sound2 = new Tbl(tmp) ; // c_s^2 = dp/de
377 
378  nbmin = pow( double(10), (*lognb)(0) ) ;
379  nbmax = pow( double(10), (*lognb)(nbp-1) ) ;
380 
381  cout << "nbmin, nbmax : " << nbmin << " " << nbmax << endl ;
382  }
383 
384 }
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
virtual int identify() const =0
Returns a number to identify the sub-classe of Dyn_eos the object belongs to.
virtual ostream & operator>>(ostream &) const
Operator >>
Definition: dyneos_cons.C:155
virtual bool operator==(const Dyn_eos &) const
Comparison operator (egality)
Definition: dyneos_cons.C:132
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
string tablename
Name of the file containing the tabulated data.
Definition: dyneos.h:655
Tbl * loge
Table of .
Definition: dyneos.h:671
Tbl * dlesdlnb
Table of .
Definition: dyneos.h:674
Tbl * c_sound2
Table of .
Definition: dyneos.h:677
double nbmin
Lower boundary of the baryon density interval.
Definition: dyneos.h:662
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Equation of state for use in dynamical code base class.
Definition: dyneos.h:75
Class for tabulated equations of state for use in dynamical code.
Definition: dyneos.h:648
virtual ~Dyn_eos_cons()
Destructor.
Definition: dyneos_cons.C:123
Tbl * lognb
Table of .
Definition: dyneos.h:668
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
void compute_derivative(const Tbl &xx, const Tbl &ff, Tbl &dfdx)
Derives a function defined on an unequally-spaced grid, approximating it by piecewise parabolae...
Definition: misc.C:64
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
string name
EOS name.
Definition: dyneos.h:81
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition: fread_be.C:72
string authors
Authors - reference for the table.
Definition: dyneos.h:657
virtual void read_table_lorene()
Reads the file containing the table in LORENE format and initializes the arrays lognb ...
Definition: dyneos_cons.C:172
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:325
virtual int identify() const
Returns a number to identify the sub-classe of Dyn_eos the object belongs to.
Definition: dyneos.C:317
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:417
virtual bool operator!=(const Dyn_eos &) const
Comparison operator (difference)
Definition: dyneos_cons.C:145
Basic array class.
Definition: tbl.h:164
bool compose_format
Are(is) the table(s) in CompOSE format?
Definition: dyneos.h:659
double nbmax
Upper boundary of the baryon density interval.
Definition: dyneos.h:665
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:542
virtual void read_table_compose()
Reads the files .nb and .thermo containing the table in CompOSE format and initializes the arrays log...
Definition: dyneos_cons.C:274
Dyn_eos_cons(const string &name_i, const string &table_name, bool compose=true)
Standard constructor.
Definition: dyneos_cons.C:65