LORENE
eos_compose.C
1 /*
2  * Methods of class Eos_CompOSE
3  *
4  * (see file eos_compose.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2010, 2014 Jerome Novak
10  * (c) 2019 Lorenzo Sala
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 /*
34  * $Id: eos_compose.C,v 1.13 2022/03/22 13:36:00 j_novak Exp $
35  * $Log: eos_compose.C,v $
36  * Revision 1.13 2022/03/22 13:36:00 j_novak
37  * Added declaration of compute_derivative to utilitaires.h
38  *
39  * Revision 1.12 2022/03/10 16:38:39 j_novak
40  * log(cs^2) is tabulated instead of cs^2.
41  *
42  * Revision 1.11 2021/05/14 15:39:22 g_servignat
43  * Added sound speed computation from enthalpy to Eos class and tabulated+polytropic derived classes
44  *
45  * Revision 1.10 2019/03/28 13:41:02 j_novak
46  * Improved managed of saved EoS (functions sauve and constructor form FILE*)
47  *
48  * Revision 1.9 2019/03/25 14:21:24 j_novak
49  * Improved constructor from a FILE*, following patch by Lorenzo Sala.
50  *
51  * Revision 1.8 2016/12/05 16:17:51 j_novak
52  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
53  *
54  * Revision 1.7 2015/12/04 16:27:05 j_novak
55  * Correction of constructor calling.
56  *
57  * Revision 1.6 2015/08/04 14:41:29 j_novak
58  * Back to previous version for Eos_CompOSE. Enthalpy-consistent EoS can be accessed using Eos_consistent class (derived from Eos_CompOSE).
59  *
60  * Revision 1.5 2015/01/27 14:22:38 j_novak
61  * New methods in Eos_tabul to correct for EoS themro consistency (optional).
62  *
63  * Revision 1.4 2014/10/13 08:52:52 j_novak
64  * Lorene classes and functions now belong to the namespace Lorene.
65  *
66  * Revision 1.3 2014/07/01 09:26:21 j_novak
67  * Improvement of comments
68  *
69  * Revision 1.2 2014/06/30 16:13:18 j_novak
70  * New methods for reading directly from CompOSE files.
71  *
72  * Revision 1.1 2014/03/06 15:53:34 j_novak
73  * Eos_compstar is now Eos_compOSE. Eos_tabul uses strings and contains informations about authors.
74  *
75  * Revision 1.1 2010/02/03 14:56:45 j_novak
76  * *** empty log message ***
77  *
78  *
79  * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_compose.C,v 1.13 2022/03/22 13:36:00 j_novak Exp $
80  *
81  */
82 
83 #include <string>
84 
85 // Headers Lorene
86 #include "headcpp.h"
87 #include "eos.h"
88 #include "tbl.h"
89 #include "utilitaires.h"
90 #include "unites.h"
91 
92  //----------------------------//
93  // Constructors //
94  //----------------------------//
95 
96 // Standard constructor
97 // --------------------
98 namespace Lorene {
99 
100 // Standard constructor with only the name of the file
101 //------------------------------------------------------
102  Eos_CompOSE::Eos_CompOSE(const char* file_name)
103  : Eos_tabul("CompOSE EoS"), format(0) {
104  tablename = file_name ;
105  }
106 
107 
108 // Constructor from binary file
109 // ----------------------------
110  Eos_CompOSE::Eos_CompOSE(FILE* fich) : Eos_tabul("CompOSE Eos"), format(0) {
111 
112  fread(name, sizeof(char), 100, fich) ;
113  char tmp_string[160] ;
114  fread(tmp_string, sizeof(char), 160, fich) ;
115  tablename = tmp_string;
116 
117  int format_read = int(fread_be(&format, sizeof(int), 1, fich)) ;
118 
119  if (format_read != 1) format = 0 ;
120 
121  if (format == 0) read_table() ;
122  else read_compose_data() ;
123  }
124 
125 
126 
127 // Constructor from a formatted file
128 // ---------------------------------
129  Eos_CompOSE::Eos_CompOSE(ifstream& fich) : Eos_tabul(fich), format(0)
130  {}
131 
132 
133 // Constructor from CompOSE data files
134 // ------------------------------------
135  Eos_CompOSE::Eos_CompOSE(const string& files) :
136  Eos_tabul("CompOSE Eos"), format(1) {
137 
138  tablename = files ;
140  }
141 
142 
143 // Reading the data files in CompOSE format
144 //-----------------------------------------
146 
147  using namespace Unites ;
148 
149  // Files containing data and a test
150  //---------------------------------
151  string file_nb = tablename + ".nb" ;
152  string file_thermo = tablename + ".thermo" ;
153 
154  ifstream in_nb(file_nb.data()) ;
155  if (!in_nb) {
156  cerr << "Eos_CompOSE::read_compose_data(): " << endl ;
157  cerr << "Problem in opening the EOS file!" << endl ;
158  cerr << "While trying to open " << file_nb << endl ;
159  cerr << "Aborting..." << endl ;
160  abort() ;
161  }
162 
163  // obtaining the size of the tables for memory allocation
164  //-------------------------------------------------------
165  int index1, index2 ;
166  in_nb >> index1 >> index2 ;
167  int nbp = index2 - index1 + 1 ;
168  assert(nbp > 0) ;
169 
170  press = new double[nbp] ;
171  nb = new double[nbp] ;
172  ro = new double[nbp] ;
173 
174  Tbl press_tbl(nbp) ; press_tbl.set_etat_qcq() ;
175  Tbl nb_tbl(nbp) ; nb_tbl.set_etat_qcq() ;
176  Tbl ro_tbl(nbp) ; ro_tbl.set_etat_qcq() ;
177 
178  logh = new Tbl(nbp) ;
179  logp = new Tbl(nbp) ;
180  dlpsdlh = new Tbl(nbp) ;
181  lognb = new Tbl(nbp) ;
182  dlpsdlnb = new Tbl(nbp) ;
183 
184  logh->set_etat_qcq() ;
185  logp->set_etat_qcq() ;
186  dlpsdlh->set_etat_qcq() ;
187  lognb->set_etat_qcq() ;
189 
190  // Variables and conversion
191  //-------------------------
192  double nb_fm3, rho_cgs, p_cgs, p_over_nb_comp, eps_comp ;
193  double dummy_x ;
194  int dummy_n ;
195 
196  double rhonuc_cgs = rhonuc_si * 1e-3 ;
197  double c2_cgs = c_si * c_si * 1e4 ;
198  double m_neutron_MeV, m_proton_MeV ;
199 
200  ifstream in_p_rho (file_thermo.data()) ;
201  if (!in_p_rho) {
202  cerr << "Eos_CompOSE::read_compose_data(): " << endl ;
203  cerr << "Problem in opening the EOS file!" << endl ;
204  cerr << "While trying to open " << file_thermo << endl ;
205  cerr << "Aborting..." << endl ;
206  abort() ;
207  }
208  in_p_rho >> m_neutron_MeV >> m_proton_MeV ; //Neutron and proton masses
209  in_p_rho.ignore(1000, '\n') ;
210 
211  double p_convert = mev_si * 1.e45 * 10. ; // Conversion from MeV/fm^3 to cgs
212  double eps_convert = mev_si * 1.e42 / (c_si*c_si) ; //From meV/fm^3 to g/cm^3
213 
214  // Main loop reading the table
215  //----------------------------
216  for (int i=0; i<nbp; i++) {
217  in_nb >> nb_fm3 ;
218  in_p_rho >> dummy_n >> dummy_n >> dummy_n >> p_over_nb_comp ;
219  in_p_rho >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> eps_comp ;
220  in_p_rho.ignore(1000, '\n') ;
221  p_cgs = p_over_nb_comp * nb_fm3 * p_convert ;
222  rho_cgs = ( eps_comp + 1. ) * m_neutron_MeV * nb_fm3 * eps_convert ;
223 
224  if ( (nb_fm3<0) || (rho_cgs<0) || (p_cgs < 0) ){
225  cout << "Eos_CompOSE::read_compose_data(): " << endl ;
226  cout << "Negative value in table!" << endl ;
227  cout << "nb = " << nb_fm3 << ", rho = " << rho_cgs <<
228  ", p = " << p_cgs << endl ;
229  cout << "Aborting..." << endl ;
230  abort() ;
231  }
232 
233  press[i] = p_cgs / c2_cgs ;
234  nb[i] = nb_fm3 ;
235  ro[i] = rho_cgs ;
236 
237  press_tbl.set(i) = press[i] ;
238  nb_tbl.set(i) = nb[i] ;
239  ro_tbl.set(i) = ro[i] ;
240  }
241 
242  double ww = 0. ;
243  for (int i=0; i<nbp; i++) {
244  double h = log( (ro[i] + press[i]) / (10 * nb[i] * rhonuc_cgs) ) ;
245 
246  if (i==0) { ww = h ; }
247  h = h - ww + 1.e-14 ;
248 
249  logh->set(i) = log10( h ) ;
250  logp->set(i) = log10( press[i] / rhonuc_cgs ) ;
251  dlpsdlh->set(i) = h * (ro[i] + press[i]) / press[i] ;
252  lognb->set(i) = log10(nb[i]) ;
253  }
254 
255  // Computation of dpdnb
256  //---------------------
257  Tbl logn0 = *lognb * log(10.) ;
258  Tbl logp0 = ((*logp) + log10(rhonuc_cgs)) * log(10.) ;
259  compute_derivative(logn0, logp0, *dlpsdlnb) ;
260 
261  // Computation of the sound speed
262  //-------------------------------
263  Tbl tmp(nbp) ; tmp.set_etat_qcq() ;
264  compute_derivative(ro_tbl,press_tbl,tmp) ;
265  for (int i=0; i<nbp; i++) {
266  if (tmp(i) < 0.) {
267  tmp.set(i) = - tmp(i) ;
268  }
269  }
270  log_cs2 = new Tbl(log10(tmp)) ;
271 
272  hmin = pow( double(10), (*logh)(0) ) ;
273  hmax = pow( double(10), (*logh)(nbp-1) ) ;
274  cout << "hmin, hmax : " << hmin << " " << hmax << endl ;
275 
276  // Cleaning
277  //---------
278  delete [] press ;
279  delete [] nb ;
280  delete [] ro ;
281 
282  }
283 
284  //--------------//
285  // Destructor //
286  //--------------//
287 
289 
290  // does nothing
291 
292  }
293 
294 
295  //------------------------//
296  // Comparison operators //
297  //------------------------//
298 
299 
300  bool Eos_CompOSE::operator==(const Eos& eos_i) const {
301 
302  bool resu = true ;
303 
304  if ( eos_i.identify() != identify() ) {
305  cout << "The second EOS is not of type Eos_CompOSE !" << endl ;
306  resu = false ;
307  }
308 
309  return resu ;
310 
311  }
312 
313  bool Eos_CompOSE::operator!=(const Eos& eos_i) const {
314 
315  return !(operator==(eos_i)) ;
316 
317  }
318 
319  //------------//
320  // Outputs //
321  //------------//
322  void Eos_CompOSE::sauve(FILE* fich) const {
323 
324  Eos_tabul::sauve(fich) ;
325 
326  fwrite_be(&format, sizeof(int), 1, fich) ;
327  }
328 
329 
330 ostream& Eos_CompOSE::operator>>(ostream & ost) const {
331 
332  ost << "EOS of class Eos_CompOSE." << endl ;
333  ost << "Built from file " << tablename << endl ;
334  ost << ((format == 0) ? "Old LORENE format" : "CompOSE format") << endl ;
335  ost << "Authors : " << authors << endl ;
336  ost << "Number of points in file : " << logh->get_taille() << endl ;
337  return ost ;
338 
339 }
340 
341 
342 }
virtual bool operator==(const Eos &) const
Comparison operator (egality)
Definition: eos_compose.C:300
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
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
Eos_CompOSE(const string &files_path)
Constructor from CompOSE data.
Definition: eos_compose.C:135
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
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
Tbl * dlpsdlnb
Table of .
Definition: eos_tabul.h:215
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
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 void sauve(FILE *) const
Save in a file.
Definition: eos_tabul.C:248
Tbl * lognb
Table of .
Definition: eos_tabul.h:212
virtual void sauve(FILE *) const
Save in a file.
Definition: eos_compose.C:322
double hmax
Upper boundary of the enthalpy interval.
Definition: eos_tabul.h:200
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
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition: fwrite_be.C:73
Base class for tabulated equations of state.
Definition: eos_tabul.h:185
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
virtual void read_compose_data()
Reads the files containing the table and initializes in the arrays logh , logp and dlpsdlh (CompOSE f...
Definition: eos_compose.C:145
virtual void read_table()
Reads the file containing the table and initializes in the arrays logh , logp and dlpsdlh ...
Definition: eos_tabul.C:261
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
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:325
virtual ostream & operator>>(ostream &) const
Operator >>
Definition: eos_compose.C:330
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:417
char name[100]
EOS name.
Definition: eos.h:215
Basic array class.
Definition: tbl.h:164
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
Definition: eos_compose.C:313
virtual ~Eos_CompOSE()
Destructor.
Definition: eos_compose.C:288
Tbl * logh
Table of .
Definition: eos_tabul.h:203