LORENE
eos.C
1 /*
2  * Methods of class Eos.
3  *
4  * (see file eos.h for documentation).
5  */
6 
7 /*
8  * Copyright (c) 2000-2001 Eric Gourgoulhon
9  *
10  * This file is part of LORENE.
11  *
12  * LORENE is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License as published by
14  * the Free Software Foundation; either version 2 of the License, or
15  * (at your option) any later version.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 
31 /*
32  * $Id: eos.C,v 1.10 2021/05/14 15:39:22 g_servignat Exp $
33  * $Log: eos.C,v $
34  * Revision 1.10 2021/05/14 15:39:22 g_servignat
35  * Added sound speed computation from enthalpy to Eos class and tabulated+polytropic derived classes
36  *
37  * Revision 1.9 2017/12/15 15:36:38 j_novak
38  * Improvement of the MEos class. Implementation of automatic offset computation accross different EoSs/domains.
39  *
40  * Revision 1.8 2016/12/05 16:17:50 j_novak
41  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
42  *
43  * Revision 1.7 2014/10/13 08:52:51 j_novak
44  * Lorene classes and functions now belong to the namespace Lorene.
45  *
46  * Revision 1.6 2014/10/06 15:13:06 j_novak
47  * Modified #include directives to use c++ syntax.
48  *
49  * Revision 1.5 2004/01/14 15:59:42 f_limousin
50  * Add methos calcule, nbar_ent, der_bar_ent, der_press_ent and press_ent
51  * for Scalar's.
52  *
53  * Revision 1.4 2002/10/16 14:36:34 j_novak
54  * Reorganization of #include instructions of standard C++, in order to
55  * use experimental version 3 of gcc.
56  *
57  * Revision 1.3 2002/04/09 14:32:15 e_gourgoulhon
58  * 1/ Added extra parameters in EOS computational functions (argument par)
59  * 2/ New class MEos for multi-domain EOS
60  *
61  * Revision 1.2 2001/12/04 21:27:53 e_gourgoulhon
62  *
63  * All writing/reading to a binary file are now performed according to
64  * the big endian convention, whatever the system is big endian or
65  * small endian, thanks to the functions fwrite_be and fread_be
66  *
67  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
68  * LORENE
69  *
70  * Revision 2.5 2001/10/10 13:45:06 eric
71  * Modif Joachim : &(Eos::der_press_ent_p) -> &Eos::der_press_ent_p, etc...
72  * pour conformite au compilateur HP.
73  *
74  * Revision 2.4 2001/02/07 09:50:49 eric
75  * Suppression de la fonction derent_ent_p.
76  * Ajout des fonctions donnant les derivees de l'EOS:
77  * der_nbar_ent_p
78  * der_ener_ent_p
79  * der_press_ent_p
80  * ainsi que des fonctions Cmp associees.
81  *
82  * Revision 2.3 2000/02/14 14:33:05 eric
83  * Ajout des constructeurs par lecture de fichier formate.
84  *
85  * Revision 2.2 2000/01/21 15:17:28 eric
86  * fonction sauve: on ecrit en premier l'identificateur.
87  *
88  * Revision 2.1 2000/01/18 13:47:07 eric
89  * Premiere version operationnelle
90  *
91  * Revision 2.0 2000/01/18 10:46:15 eric
92  * /
93  *
94  *
95  * $Header: /cvsroot/Lorene/C++/Source/Eos/eos.C,v 1.10 2021/05/14 15:39:22 g_servignat Exp $
96  *
97  */
98 
99 // Headers C
100 #include <cstdlib>
101 #include <cstring>
102 
103 // Headers Lorene
104 #include "eos.h"
105 #include "cmp.h"
106 #include "scalar.h"
107 #include "utilitaires.h"
108 #include "param.h"
109 
110  //--------------//
111  // Constructors //
112  //--------------//
113 
114 
115 // Standard constructor without name
116 // ---------------------------------
117 namespace Lorene {
119 
120  set_name("") ;
121 
122 }
123 
124 // Standard constructor with name
125 // ---------------------------------
126 Eos::Eos(const char* name_i){
127 
128  set_name(name_i) ;
129 
130 }
131 
132 // Copy constructor
133 // ----------------
134 Eos::Eos(const Eos& eos_i){
135 
136  set_name(eos_i.name) ;
137 
138 }
139 
140 // Constructor from a binary file
141 // ------------------------------
142 Eos::Eos(FILE* fich){
143 
144  fread(name, sizeof(char), 100, fich) ;
145 
146 }
147 
148 // Constructor from a formatted file
149 // ---------------------------------
150 Eos::Eos(ifstream& fich){
151 
152  fich.getline(name, 100) ;
153 
154 }
155 
156 
157 
158  //--------------//
159  // Destructor //
160  //--------------//
161 
163 
164  // does nothing
165 
166 }
167 
168  //-------------------------//
169  // Manipulation of name //
170  //-------------------------//
171 
172 
173 void Eos::set_name(const char* name_i) {
174 
175  strncpy(name, name_i, 100) ;
176 
177 }
178 
179 const char* Eos::get_name() const {
180 
181  return name ;
182 
183 }
184 
185  //------------//
186  // Outputs //
187  //------------//
188 
189 void Eos::sauve(FILE* fich) const {
190 
191  int ident = identify() ;
192  fwrite_be(&ident, sizeof(int), 1, fich) ;
193 
194  fwrite(name, sizeof(char), 100, fich) ;
195 
196 }
197 
198 
199 
200 
201 ostream& operator<<(ostream& ost, const Eos& eqetat) {
202  ost << eqetat.get_name() << endl ;
203  eqetat >> ost ;
204  return ost ;
205 }
206 
207 
208  //-------------------------------//
209  // Generic computational routine //
210  //-------------------------------//
211 
212 
213 void Eos::calcule(const Cmp& ent, int nzet, int l_min,
214  double (Eos::*fait)(double, const Param*) const, Param* par, Cmp& resu) const {
215 
216  assert(ent.get_etat() != ETATNONDEF) ;
217 
218  const Map* mp = ent.get_mp() ; // Mapping
219 
220 
221  if (ent.get_etat() == ETATZERO) {
222  resu.set_etat_zero() ;
223  return ;
224  }
225 
226  assert(ent.get_etat() == ETATQCQ) ;
227  const MEos* tryMEos = dynamic_cast<const MEos*>(this) ;
228  bool isMEos = (tryMEos != 0x0) ;
229  int l_index = 0 ; // Index of the current domain in case of MEos
230  if (isMEos) par->add_int_mod(l_index) ;
231 
232  const Valeur& vent = ent.va ;
233  vent.coef_i() ; // the values in the configuration space are required
234 
235  const Mg3d* mg = mp->get_mg() ; // Multi-grid
236 
237  int nz = mg->get_nzone() ; // total number of domains
238 
239  // Preparations for a point by point computation:
240  resu.set_etat_qcq() ;
241  Valeur& vresu = resu.va ;
242  vresu.set_etat_c_qcq() ;
243  vresu.c->set_etat_qcq() ;
244 
245  // Loop on domains where the computation has to be done :
246  for (int l = l_min; l< l_min + nzet; l++) {
247 
248  assert(l>=0) ;
249  assert(l<nz) ;
250 
251  l_index = l ; // The domain index is passed to the 'fait' function
252 
253  Tbl* tent = vent.c->t[l] ;
254  Tbl* tresu = vresu.c->t[l] ;
255 
256  if (tent->get_etat() == ETATZERO) {
257  tresu->set_etat_zero() ;
258  }
259  else {
260  assert( tent->get_etat() == ETATQCQ ) ;
261  tresu->set_etat_qcq() ;
262 
263  for (int i=0; i<tent->get_taille(); i++) {
264 
265  tresu->t[i] = (this->*fait)( tent->t[i], par ) ;
266  }
267 
268  } // End of the case where ent != 0 in the considered domain
269 
270  } // End of the loop on domains where the computation had to be done
271 
272  // resu is set to zero in the other domains :
273 
274  if (l_min > 0) {
275  resu.annule(0, l_min-1) ;
276  }
277 
278  if (l_min + nzet < nz) {
279  resu.annule(l_min + nzet, nz - 1) ;
280  }
281 }
282 
283 
284 
285 void Eos::calcule(const Scalar& ent, int nzet, int l_min,
286  double (Eos::*fait)(double, const Param*) const, Param* par, Scalar& resu) const {
287 
288  assert(ent.get_etat() != ETATNONDEF) ;
289 
290  const Map* mp = &(ent.get_mp()) ; // Mapping
291 
292 
293  if (ent.get_etat() == ETATZERO) {
294  resu.set_etat_zero() ;
295  return ;
296  }
297 
298  assert(ent.get_etat() == ETATQCQ) ;
299  const MEos* tryMEos = dynamic_cast<const MEos*>(this) ;
300  bool isMEos = (tryMEos != 0x0) ;
301  int l_index = 0 ; // Index of the current domain in case of MEos
302  if (isMEos) par->add_int_mod(l_index) ;
303 
304  const Valeur& vent = ent.get_spectral_va() ;
305  vent.coef_i() ; // the values in the configuration space are required
306 
307  const Mg3d* mg = mp->get_mg() ; // Multi-grid
308 
309  int nz = mg->get_nzone() ; // total number of domains
310 
311  // Preparations for a point by point computation:
312  resu.set_etat_qcq() ;
313  Valeur& vresu = resu.set_spectral_va() ;
314  vresu.set_etat_c_qcq() ;
315  vresu.c->set_etat_qcq() ;
316 
317  // Loop on domains where the computation has to be done :
318  for (int l = l_min; l< l_min + nzet; l++) {
319 
320  assert(l>=0) ;
321  assert(l<nz) ;
322 
323  l_index = l ; // The domain index is passed to the 'fait' function
324 
325  Tbl* tent = vent.c->t[l] ;
326  Tbl* tresu = vresu.c->t[l] ;
327 
328  if (tent->get_etat() == ETATZERO) {
329  tresu->set_etat_zero() ;
330  }
331  else {
332  assert( tent->get_etat() == ETATQCQ ) ;
333  tresu->set_etat_qcq() ;
334 
335  for (int i=0; i<tent->get_taille(); i++) {
336 
337  tresu->t[i] = (this->*fait)( tent->t[i], par ) ;
338  }
339 
340  } // End of the case where ent != 0 in the considered domain
341 
342  } // End of the loop on domains where the computation had to be done
343 
344  // resu is set to zero in the other domains :
345 
346  if (l_min > 0) {
347  resu.annule(0, l_min-1) ;
348  }
349 
350  if (l_min + nzet < nz) {
351  resu.annule(l_min + nzet, nz - 1) ;
352  }
353 }
354  //---------------------------------//
355  // Public functions //
356  //---------------------------------//
357 
358 
359 // Baryon density from enthalpy
360 //------------------------------
361 
362 Cmp Eos::nbar_ent(const Cmp& ent, int nzet, int l_min, Param* par) const {
363 
364  Cmp resu(ent.get_mp()) ;
365 
366  calcule(ent, nzet, l_min, &Eos::nbar_ent_p, par, resu) ;
367 
368  return resu ;
369 
370 }
371 
372 Scalar Eos::nbar_ent(const Scalar& ent, int nzet, int l_min, Param* par) const {
373 
374  Scalar resu(ent.get_mp()) ;
375 
376  calcule(ent, nzet, l_min, &Eos::nbar_ent_p, par, resu) ;
377 
378  return resu ;
379 
380 }
381 
382 
383 
384 // Energy density from enthalpy
385 //------------------------------
386 
387 Cmp Eos::ener_ent(const Cmp& ent, int nzet, int l_min, Param* par) const {
388 
389  Cmp resu(ent.get_mp()) ;
390 
391  calcule(ent, nzet, l_min, &Eos::ener_ent_p, par, resu) ;
392 
393  return resu ;
394 
395 }
396 
397 Scalar Eos::ener_ent(const Scalar& ent, int nzet, int l_min, Param* par) const {
398 
399  Scalar resu(ent.get_mp()) ;
400 
401  calcule(ent, nzet, l_min, &Eos::ener_ent_p, par, resu) ;
402 
403  return resu ;
404 
405 }
406 // Pressure from enthalpy
407 //-----------------------
408 
409 Cmp Eos::press_ent(const Cmp& ent, int nzet, int l_min, Param* par) const {
410 
411  Cmp resu(ent.get_mp()) ;
412 
413  calcule(ent, nzet, l_min, &Eos::press_ent_p, par, resu) ;
414 
415  return resu ;
416 
417 }
418 
419 Scalar Eos::press_ent(const Scalar& ent, int nzet, int l_min, Param* par) const {
420 
421  Scalar resu(ent.get_mp()) ;
422 
423  calcule(ent, nzet, l_min, &Eos::press_ent_p, par, resu) ;
424 
425  return resu ;
426 
427 }
428 // Derivative of baryon density from enthalpy
429 //-------------------------------------------
430 
431 Cmp Eos::der_nbar_ent(const Cmp& ent, int nzet, int l_min, Param* par) const {
432 
433  Cmp resu(ent.get_mp()) ;
434 
435  calcule(ent, nzet, l_min, &Eos::der_nbar_ent_p, par, resu) ;
436 
437  return resu ;
438 
439 }
440 
441 Scalar Eos::der_nbar_ent(const Scalar& ent, int nzet, int l_min, Param* par) const {
442 
443  Scalar resu(ent.get_mp()) ;
444 
445  calcule(ent, nzet, l_min, &Eos::der_nbar_ent_p, par, resu) ;
446 
447  return resu ;
448 
449 }
450 
451 // Derivative of energy density from enthalpy
452 //-------------------------------------------
453 
454 Cmp Eos::der_ener_ent(const Cmp& ent, int nzet, int l_min, Param* par) const {
455 
456  Cmp resu(ent.get_mp()) ;
457 
458  calcule(ent, nzet, l_min, &Eos::der_ener_ent_p, par, resu) ;
459 
460  return resu ;
461 
462 }
463 
464 Scalar Eos::der_ener_ent(const Scalar& ent, int nzet, int l_min, Param* par) const {
465 
466  Scalar resu(ent.get_mp()) ;
467 
468  calcule(ent, nzet, l_min, &Eos::der_ener_ent_p, par, resu) ;
469 
470  return resu ;
471 
472 }
473 // Derivative of pressure from enthalpy
474 //-------------------------------------------
475 
476 Cmp Eos::der_press_ent(const Cmp& ent, int nzet, int l_min, Param* par) const {
477 
478  Cmp resu(ent.get_mp()) ;
479 
480  calcule(ent, nzet, l_min, &Eos::der_press_ent_p, par, resu) ;
481 
482  return resu ;
483 
484 }
485 
486 Scalar Eos::der_press_ent(const Scalar& ent, int nzet, int l_min, Param* par) const {
487 
488  Scalar resu(ent.get_mp()) ;
489 
490  calcule(ent, nzet, l_min, &Eos::der_press_ent_p, par, resu) ;
491 
492  return resu ;
493 
494 }
495 
496 // Sound speed from enthalpy
497 //---------------------------------
498 
499 Scalar Eos::csound_square_ent(const Scalar& ent, int nzet,
500  int l_min, Param* par) const {
501  Scalar resu(ent.get_mp()) ;
502 
503  calcule(ent, nzet, l_min, &Eos::csound_square_ent_p, par, resu) ;
504 
505  return resu ;
506  }
507 
508 }
509 
510 
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Cmp ener_ent(const Cmp &ent, int nzet, int l_min=0, Param *par=0x0) const
Computes the total energy density from the log-enthalpy and extra parameters.
Definition: eos.C:387
virtual double der_press_ent_p(double ent, const Param *par=0x0) const =0
Computes the logarithmic derivative from the log-enthalpy and extra parameters (virtual function imp...
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const =0
Computes the logarithmic derivative from the log-enthalpy and extra parameters (virtual function imp...
Cmp der_nbar_ent(const Cmp &ent, int nzet, int l_min=0, Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy and extra parameters.
Definition: eos.C:431
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
const char * get_name() const
Returns the EOS name.
Definition: eos.C:179
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:351
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
Lorene prototypes.
Definition: app_hor.h:67
Equation of state base class.
Definition: eos.h:209
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Scalar csound_square_ent(const Scalar &ent, int nzet, int l_min=0, Param *par=0x0) const
Computes the sound speed squared from the enthalpy with extra parameters.
Definition: eos.C:499
void coef_i() const
Computes the physical value of *this.
Base class for coordinate mappings.
Definition: map.h:682
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
int get_etat() const
Gives the logical state.
Definition: tbl.h:414
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:359
virtual double press_ent_p(double ent, const Param *par=0x0) const =0
Computes the pressure from the log-enthalpy and extra parameters (virtual function implemented in the...
EOS with domain dependency.
Definition: eos.h:2780
virtual void sauve(FILE *) const
Save in a file.
Definition: eos.C:189
Eos()
Standard constructor.
Definition: eos.C:118
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
virtual double ener_ent_p(double ent, const Param *par=0x0) const =0
Computes the total energy density from the log-enthalpy and extra parameters (virtual function implem...
Cmp nbar_ent(const Cmp &ent, int nzet, int l_min=0, Param *par=0x0) const
Computes the baryon density field from the log-enthalpy field and extra parameters.
Definition: eos.C:362
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:292
virtual double nbar_ent_p(double ent, const Param *par=0x0) const =0
Computes the baryon density from the log-enthalpy and extra parameters (virtual function implemented ...
double * t
The array of double.
Definition: tbl.h:176
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:309
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
virtual double csound_square_ent_p(double ent, const Param *par=0x0) const =0
Computes the sound speed squared from the enthapy with extra parameters (virtual function implemente...
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
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:307
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:302
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: tbl.C:350
Multi-domain grid.
Definition: grilles.h:279
virtual ~Eos()
Destructor.
Definition: eos.C:162
void set_name(const char *name_i)
Sets the EOS name.
Definition: eos.C:173
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:417
Cmp der_ener_ent(const Cmp &ent, int nzet, int l_min=0, Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy and extra parameters.
Definition: eos.C:454
char name[100]
EOS name.
Definition: eos.h:215
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c )...
Definition: valeur.C:704
Basic array class.
Definition: tbl.h:164
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Cmp der_press_ent(const Cmp &ent, int nzet, int l_min=0, Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy and extra parameters.
Definition: eos.C:476
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s.
Definition: mtbl.h:132
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:874
void calcule(const Cmp &thermo, int nzet, int l_min, double(Eos::*fait)(double, const Param *) const, Param *par, Cmp &resu) const
General computational method for Cmp &#39;s.
Definition: eos.C:213
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const =0
Computes the logarithmic derivative from the log-enthalpy with extra parameters (virtual function im...
Cmp press_ent(const Cmp &ent, int nzet, int l_min=0, Param *par=0x0) const
Computes the pressure from the log-enthalpy and extra parameters.
Definition: eos.C:409
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:388
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607