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