LORENE
hoteos_tabul.C
1 /*
2  * Methods of class Hoteos_tabul
3  *
4  * (see file hoteos.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 hoteos_tabul_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/hoteos_tabul.C,v 1.2 2015/12/08 15:42:17 j_novak Exp $" ;
31 
32 /*
33  * $Id: hoteos_tabul.C,v 1.2 2015/12/08 15:42:17 j_novak Exp $
34  * $Log: hoteos_tabul.C,v $
35  * Revision 1.2 2015/12/08 15:42:17 j_novak
36  * Low/zero entropy is set to the lowest value in the table in computational functions.
37  *
38  * Revision 1.1 2015/12/08 10:52:18 j_novak
39  * New class Hoteos_tabul for tabulated temperature-dependent EoSs.
40  *
41  *
42  *
43  * $Header: /cvsroot/Lorene/C++/Source/Eos/hoteos_tabul.C,v 1.2 2015/12/08 15:42:17 j_novak Exp $
44  *
45  */
46 
47 // headers C
48 #include <cstdlib>
49 #include <cstring>
50 #include <cmath>
51 
52 // Headers Lorene
53 #include "hoteos.h"
54 #include "eos.h"
55 #include "utilitaires.h"
56 #include "unites.h"
57 
58 
59 namespace Lorene {
60  void interpol_herm_2d(const Tbl&, const Tbl&, const Tbl&, const Tbl&, const Tbl&,
61  const Tbl&, double, double, double&, double&, double&) ;
62 
63 
64  //----------------------------//
65  // Constructors //
66  //----------------------------//
67 
68  // Standard constructor
69  // --------------------
70  Hoteos_tabul::Hoteos_tabul(const string& filename):
71  Hot_eos("Tabulated hot EoS"), tablename(filename), authors("Unknown"),
72  hmin(0.), hmax(1.), sbmin(0.), sbmax(1.)
73  {
74  set_arrays_0x0() ;
75  read_table() ;
76  }
77 
78 
79  // Constructor from binary file
80  // ----------------------------
81  Hoteos_tabul::Hoteos_tabul(FILE* fich) : Hot_eos(fich) {
82 
83  char tmp_string[160] ;
84  fread(tmp_string, sizeof(char), 160, fich) ;
85  tablename = tmp_string ;
86  set_arrays_0x0() ;
87  read_table() ;
88  }
89 
90  // Constructor from a formatted file
91  // ---------------------------------
92  Hoteos_tabul::Hoteos_tabul(ifstream& fich) :
93  Hot_eos(fich) {
94 
95  fich >> tablename ;
96  set_arrays_0x0() ;
97  read_table() ;
98  }
99 
100  //Sets the arrays to the null pointer
102  hhh = 0x0 ;
103  s_B = 0x0 ;
104  ppp = 0x0 ;
105  dpdh = 0x0 ;
106  dpds = 0x0 ;
107  d2p = 0x0 ;
108  }
109 
110  //--------------//
111  // Destructor //
112  //--------------//
113 
115  if (hhh != 0x0) delete hhh ;
116  if (s_B != 0x0) delete s_B ;
117  if (ppp != 0x0) delete ppp ;
118  if (dpdh != 0x0) delete dpdh ;
119  if (dpds != 0x0) delete dpds ;
120  if (d2p != 0x0) delete d2p ;
121  }
122 
123  //------------//
124  // Outputs //
125  //------------//
126 
127  void Hoteos_tabul::sauve(FILE* fich) const {
128 
129  Hot_eos::sauve(fich) ;
130 
131  char tmp_string[160] ;
132  strcpy(tmp_string, tablename.c_str()) ;
133  fwrite(tmp_string, sizeof(char), 160, fich) ;
134  }
135 
136  ostream& Hoteos_tabul::operator>>(ostream & ost) const {
137 
138  ost << "Hot EOS of class Hoteos_tabul (tabulated hot beta-equilibrium EoS) : "
139  << endl ;
140  ost << "Built from file " << tablename << endl ;
141  ost << "Authors : " << authors << endl ;
142  ost << "Number of points in file : " << hhh->get_dim(0)
143  << " in enthalpy, and " << hhh->get_dim(1)
144  << " in entropy." << endl ;
145 
146  return ost ;
147 }
148 
149  //------------------------//
150  // Reading of the table //
151  //------------------------//
152 
154 
155  using namespace Unites ;
156 
157  ifstream fich(tablename.data()) ;
158 
159  if (!fich) {
160  cerr << "Hoteos_tabul::read_table(): " << endl ;
161  cerr << "Problem in opening the EOS file!" << endl ;
162  cerr << "While trying to open " << tablename << endl ;
163  cerr << "Aborting..." << endl ;
164  abort() ;
165  }
166 
167  fich.ignore(1000, '\n') ; // Jump over the first header
168  fich.ignore(1) ;
169  getline(fich, authors, '\n') ;
170  for (int i=0; i<3; i++) { // jump over the file
171  fich.ignore(1000, '\n') ; // header
172  } //
173 
174  int nbp1, nbp2 ;
175  fich >> nbp1 >> nbp2 ; fich.ignore(1000, '\n') ; // number of data
176  if ( (nbp1<=0) || (nbp2<=0) ) {
177  cerr << "Hoteos_tabul::read_table(): " << endl ;
178  cerr << "Wrong value for the number of lines!" << endl ;
179  cerr << "nbp1 = " << nbp1 << ", nbp2 = " << nbp2 << endl ;
180  cerr << "Aborting..." << endl ;
181  abort() ;
182  }
183 
184  for (int i=0; i<3; i++) { // jump over the table
185  fich.ignore(1000, '\n') ; // description
186  }
187 
188  ppp = new Tbl(nbp2, nbp1) ;
189  hhh = new Tbl(nbp2, nbp1) ;
190  s_B = new Tbl(nbp2, nbp1) ;
191  dpdh = new Tbl(nbp2, nbp1) ;
192  dpds = new Tbl(nbp2, nbp1) ;
193  d2p = new Tbl(nbp2, nbp1) ;
194 
195  ppp->set_etat_qcq() ;
196  hhh->set_etat_qcq() ;
197  s_B->set_etat_qcq() ;
198  dpdh->set_etat_qcq() ;
199  dpds->set_etat_qcq() ;
200  d2p->set_etat_qcq() ;
201 
202  double c2 = c_si * c_si ;
203  double dummy, nb_fm3, rho_cgs, p_cgs, mu_MeV, entr, temp, der2 ;
204  double ww = 0. ;
205 
206  for (int j=0; j<nbp2; j++) {
207  for (int i=0; i<nbp1; i++) {
208  fich >> mu_MeV >> entr >> nb_fm3 >> temp >> p_cgs >> der2
209  >> dummy >> rho_cgs ;
210  fich.ignore(1000,'\n') ;
211  if ( (nb_fm3<0) || (rho_cgs<0) || (p_cgs < 0) ){
212  cerr << "Eos_mag::read_table(): " << endl ;
213  cerr << "Negative value in table!" << endl ;
214  cerr << "nb = " << nb_fm3 << ", rho = " << rho_cgs <<
215  ", p = " << p_cgs << endl ;
216  cerr << "Aborting..." << endl ;
217  abort() ;
218  }
219 
220  double psc2 = 0.1*p_cgs/c2 ; // in kg m^-3
221  double rho_si = rho_cgs*1000. ; // in kg m^-3
222 
223  double h_read = log(mu_MeV) ;
224  if ( (i==0) && (j==0) ) ww = h_read ;
225  double h_new = h_read - ww ;
226 
227  ppp->set(j, i) = psc2/rhonuc_si ;
228  hhh->set(j, i) = h_new ;
229  s_B->set(j, i) = entr ; // in Lorene units (k_B)
230  dpdh->set(j, i) = (rho_si + psc2)/rhonuc_si ;
231  dpds->set(j, i) = -temp*nb_fm3*mevpfm3 ;
232  d2p->set(j, i) = 0.1*der2*mu_MeV/(c2*rhonuc_si) ;
233  }
234  }
235 
236  hmin = (*hhh)(0, 0) ;
237  hmax = (*hhh)(0, nbp1-1) ;
238 
239  sbmin = (*s_B)(0, 0) ;
240  sbmax = (*s_B)(nbp2-1, 0) ;
241 
242  cout << "hmin: " << hmin << ", hmax: " << hmax << endl ;
243  cout << "sbmin: " << sbmin << ", sbmax: " << sbmax << endl ;
244 
245  fich.close();
246 }
247 
248  //-------------------------------//
249  // The corresponding cold Eos //
250  //-------------------------------//
251 
253 
254  if (p_cold_eos == 0x0) {
255  cerr << "Warning: Hoteos_tabul::new_cold_Eos " <<
256  "The corresponding cold EoS is likely not to function." << endl ;
257  p_cold_eos = new Eos_CompOSE(tablename.c_str()) ;
258  }
259 
260  return *p_cold_eos ;
261  }
262 
263 
264  //------------------------//
265  // Comparison operators //
266  //------------------------//
267 
268 
269  bool Hoteos_tabul::operator==(const Hot_eos& eos_i) const {
270 
271  bool resu = true ;
272 
273  if ( eos_i.identify() != identify() ) {
274  cout << "The second EOS is not of type Hoteos_tabul !" << endl ;
275  resu = false ;
276  }
277 
278  return resu ;
279  }
280 
281 
282  bool Hoteos_tabul::operator!=(const Hot_eos& eos_i) const {
283  return !(operator==(eos_i)) ;
284  }
285 
286 
287  //------------------------------//
288  // Computational routines //
289  //------------------------------//
290 
291 // Baryon density from enthalpy and entropy
292 //------------------------------------------
293 
294 double Hoteos_tabul::nbar_Hs_p(double ent, double sb) const {
295 
296  if ((ent > hmin - 1.e-12) && (ent < hmin))
297  ent = hmin ;
298 
299  if (sb < sbmin) sb = sbmin ;
300 
301  if ( ent >= hmin ) {
302  if (ent > hmax) {
303  cout << "Hoteos_tabul::nbar_Hs_p : ent > hmax !" << endl ;
304  abort() ;
305  }
306 
307  if (sb > sbmax) {
308  cerr << "Hoteos_tabul::nbar_Hs_p : s_B not in the tabulated interval !"
309  << endl ;
310  cerr << "s_B = " << sb << ", sbmin = " << sbmin << ", sbmax = " << sbmax
311  << endl ;
312  abort() ;
313  }
314 
315  double p_int, dpds_int, dpdh_int ;
316  interpol_herm_2d(*s_B, *hhh, *ppp, *dpds, *dpdh, *d2p, sb, ent, p_int,
317  dpds_int, dpdh_int) ;
318 
319  double nbar_int = dpdh_int * exp(-ent) ;
320  return nbar_int ;
321  }
322  else{
323  return 0 ;
324  }
325 }
326 
327 // Energy density from enthalpy and entropy
328 //-----------------------------------------
329 
330 double Hoteos_tabul::ener_Hs_p(double ent, double sb) const {
331 
332  if ((ent > hmin - 1.e-12) && (ent < hmin))
333  ent = hmin ;
334 
335  if (sb < sbmin) sb = sbmin ;
336 
337  if ( ent >= hmin ) {
338  if (ent > hmax) {
339  cout << "Hoteos_tabul::ener_Hs_p : ent > hmax !" << endl ;
340  abort() ;
341  }
342 
343  if (sb > sbmax) {
344  cerr << "Hoteos_tabul::ener_Hs_p : s_B not in the tabulated interval !"
345  << endl ;
346  cerr << "s_B = " << sb << ", sbmin = " << sbmin << ", sbmax = " << sbmax
347  << endl ;
348  abort() ;
349  }
350 
351  double p_int, dpds_int, dpdh_int ;
352  interpol_herm_2d(*s_B, *hhh, *ppp, *dpds, *dpdh, *d2p, sb, ent, p_int,
353  dpds_int, dpdh_int) ;
354 
355  double nbar_int = dpdh_int * exp(-ent) ;
356 
357  double f_int = - p_int + exp(ent) * nbar_int;
358 
359  return f_int ;
360  }
361  else{
362  return 0 ;
363  }
364 }
365 
366 // Pressure from enthalpy and entropy
367 //-----------------------------------
368 
369 double Hoteos_tabul::press_Hs_p(double ent, double sb) const {
370 
371  if ((ent > hmin - 1.e-12) && (ent < hmin))
372  ent = hmin ;
373 
374  if (sb < sbmin) sb = sbmin ;
375 
376  if ( ent >= hmin ) {
377  if (ent > hmax) {
378  cout << "Hoteos_tabul::press_Hs_p : ent > hmax !" << endl ;
379  abort() ;
380  }
381 
382  if (sb > sbmax) {
383  cerr << "Hoteos_tabul::press_Hs_p : s_B not in the tabulated interval !"
384  << endl ;
385  cerr << "s_B = " << sb << ", sbmin = " << sbmin << ", sbmax = " << sbmax
386  << endl ;
387  abort() ;
388  }
389 
390  double p_int, dpds_int, dpdh_int ;
391  interpol_herm_2d(*s_B, *hhh, *ppp, *dpds, *dpdh, *d2p, sb, ent, p_int,
392  dpds_int, dpdh_int) ;
393 
394  return p_int ;
395  }
396  else{
397  return 0 ;
398  }
399 }
400 
401  // Temperature from enthalpy and entropy
402  //--------------------------------------
403  double Hoteos_tabul::temp_Hs_p(double ent, double sb) const {
404 
405  if ((ent > hmin - 1.e-12) && (ent < hmin))
406  ent = hmin ;
407 
408  if (sb < sbmin) sb = sbmin ;
409 
410  if ( ent >= hmin ) {
411  if (ent > hmax) {
412  cout << "Hoteos_tabul::temp_Hs_p : ent > hmax !" << endl ;
413  abort() ;
414  }
415 
416  if (sb > sbmax) {
417  cerr << "Hoteos_tabul::temp_Hs_p : s_B not in the tabulated interval !"
418  << endl ;
419  cerr << "s_B = " << sb << ", sbmin = " << sbmin << ", sbmax = " << sbmax
420  << endl ;
421  abort() ;
422  }
423 
424  double p_int, dpds_int, dpdh_int ;
425  interpol_herm_2d(*s_B, *hhh, *ppp, *dpds, *dpdh, *d2p, sb, ent, p_int,
426  dpds_int, dpdh_int) ;
427 
428  double nbar_int = dpdh_int * exp(-ent) ;
429 
430  double temp_int = -dpds_int / nbar_int ;
431 
432  return temp_int ;
433  }
434  else {
435  return 0 ;
436  }
437  }
438 
439 } //End of namespace Lorene
Lorene::Hoteos_tabul::operator==
virtual bool operator==(const Hot_eos &) const
Comparison operator (egality)
Definition: hoteos_tabul.C:269
Lorene::Hoteos_tabul::Hoteos_tabul
Hoteos_tabul(const string &filename)
Standard constructor from a filename.
Definition: hoteos_tabul.C:70
Lorene::Hoteos_tabul::ppp
Tbl * ppp
Table of pressure $P$.
Definition: hoteos.h:579
Lorene::Hoteos_tabul::operator>>
virtual ostream & operator>>(ostream &) const
Operator >>
Definition: hoteos_tabul.C:136
Lorene::Hoteos_tabul::sbmax
double sbmax
Upper boundary of the entropy interval.
Definition: hoteos.h:570
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Lorene::Hoteos_tabul::hmin
double hmin
Lower boundary of the enthalpy interval.
Definition: hoteos.h:561
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
Lorene::Hoteos_tabul::d2p
Tbl * d2p
Table of .
Definition: hoteos.h:588
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Hoteos_tabul::hmax
double hmax
Upper boundary of the enthalpy interval.
Definition: hoteos.h:564
Lorene::Hoteos_tabul::dpdh
Tbl * dpdh
Table of .
Definition: hoteos.h:582
Lorene::Eos
Equation of state base class.
Definition: eos.h:190
Lorene::Hoteos_tabul::press_Hs_p
virtual double press_Hs_p(double ent, double sb) const
Computes the pressure from the log-enthalpy and entropy per baryon (virtual function implemented in t...
Definition: hoteos_tabul.C:369
Lorene::Hoteos_tabul::read_table
void read_table()
Reads the file containing the table and initializes in the arrays hhh , s_B, ppp, ....
Definition: hoteos_tabul.C:153
Lorene::Hoteos_tabul::identify
virtual int identify() const
Returns a number to identify the sub-classe of Hot_eos the object belongs to.
Definition: hoteos_from_file.C:60
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::Tbl::get_dim
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:403
Lorene::Hoteos_tabul::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition: hoteos_tabul.C:127
Lorene::Hot_eos::identify
virtual int identify() const =0
Returns a number to identify the sub-classe of Hot_eos the object belongs to.
Lorene::Hoteos_tabul::~Hoteos_tabul
virtual ~Hoteos_tabul()
Destructor.
Definition: hoteos_tabul.C:114
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Hoteos_tabul::s_B
Tbl * s_B
Table of , entropy per baryon (in units of Boltzmann constant).
Definition: hoteos.h:576
Lorene::Hoteos_tabul::hhh
Tbl * hhh
Table of .
Definition: hoteos.h:573
Lorene::Hoteos_tabul::authors
string authors
Authors - reference for the table.
Definition: hoteos.h:558
Lorene::Hoteos_tabul::set_arrays_0x0
void set_arrays_0x0()
Sets all the arrays to the null pointer.
Definition: hoteos_tabul.C:101
Lorene::Hoteos_tabul::temp_Hs_p
virtual double temp_Hs_p(double ent, double sb) const
Computes the temperature from the log-enthalpy and entropy per baryon (virtual function implemented i...
Definition: hoteos_tabul.C:403
Lorene::Hoteos_tabul::nbar_Hs_p
virtual double nbar_Hs_p(double ent, double sb) const
Computes the baryon density from the log-enthalpy and entropy per baryon (virtual function implemente...
Definition: hoteos_tabul.C:294
Lorene::Hoteos_tabul::sbmin
double sbmin
Lower boundary of the entropy interval.
Definition: hoteos.h:567
Lorene::Hoteos_tabul::tablename
string tablename
Name of the file containing the tabulated data.
Definition: hoteos.h:556
Lorene::Hot_eos::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition: hoteos.C:129
Lorene::Hoteos_tabul::operator!=
virtual bool operator!=(const Hot_eos &) const
Comparison operator (difference)
Definition: hoteos_tabul.C:282
Lorene::Hoteos_tabul::dpds
Tbl * dpds
Table of .
Definition: hoteos.h:585
Lorene::Hoteos_tabul::ener_Hs_p
virtual double ener_Hs_p(double ent, double sb) const
Computes the total energy density from the log-enthalpy and entropy per baryon (virtual function impl...
Definition: hoteos_tabul.C:330
Lorene::Hot_eos::p_cold_eos
Eos * p_cold_eos
Corresponding cold Eos.
Definition: hoteos.h:108
Lorene::Hoteos_tabul::new_cold_Eos
virtual const Eos & new_cold_Eos() const
Returns the corresponding cold Eos.
Definition: hoteos_tabul.C:252
Lorene::Hot_eos
Base class for temperature-dependent equations of state (abstract class).
Definition: hoteos.h:67
Lorene::Eos_CompOSE
Equation of state for the CompOSE database.
Definition: eos_compose.h:77