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
30char 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
59namespace 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 {
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 ;
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 ;
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
294double 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
330double 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
369double 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
Equation of state for the CompOSE database.
Definition eos_compose.h:77
Equation of state base class.
Definition eos.h:190
Base class for temperature-dependent equations of state (abstract class).
Definition hoteos.h:67
virtual int identify() const =0
Returns a number to identify the sub-classe of Hot_eos the object belongs to.
Eos * p_cold_eos
Corresponding cold Eos.
Definition hoteos.h:108
virtual void sauve(FILE *) const
Save in a file.
Definition hoteos.C:129
double hmin
Lower boundary of the enthalpy interval.
Definition hoteos.h:561
virtual void sauve(FILE *) const
Save in a file.
virtual ~Hoteos_tabul()
Destructor.
double sbmax
Upper boundary of the entropy interval.
Definition hoteos.h:570
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...
Hoteos_tabul(const string &filename)
Standard constructor from a filename.
virtual bool operator!=(const Hot_eos &) const
Comparison operator (difference)
Tbl * hhh
Table of .
Definition hoteos.h:573
Tbl * d2p
Table of .
Definition hoteos.h:588
void set_arrays_0x0()
Sets all the arrays to the null pointer.
virtual ostream & operator>>(ostream &) const
Operator >>
virtual bool operator==(const Hot_eos &) const
Comparison operator (egality)
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...
double hmax
Upper boundary of the enthalpy interval.
Definition hoteos.h:564
virtual const Eos & new_cold_Eos() const
Returns the corresponding cold Eos.
string tablename
Name of the file containing the tabulated data.
Definition hoteos.h:556
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...
Tbl * dpdh
Table of .
Definition hoteos.h:582
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...
virtual int identify() const
Returns a number to identify the sub-classe of Hot_eos the object belongs to.
Tbl * dpds
Table of .
Definition hoteos.h:585
string authors
Authors - reference for the table.
Definition hoteos.h:558
Tbl * s_B
Table of , entropy per baryon (in units of Boltzmann constant).
Definition hoteos.h:576
void read_table()
Reads the file containing the table and initializes in the arrays hhh , s_B, ppp, ....
Tbl * ppp
Table of pressure $P$.
Definition hoteos.h:579
double sbmin
Lower boundary of the entropy interval.
Definition hoteos.h:567
Basic array class.
Definition tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition tbl.h:403
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.