LORENE
eos_mag.C
1  /* Methods of class Eos_mag
2  *
3  * (see file eos_mag.h for documentation).
4  *
5  */
6 
7 /*
8  * Copyright (c) 2011 Thomas Elghozi & Jerome Novak
9  * Copyright (c) 2013 Debarati Chatterjee
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_mag_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/eos_mag.C,v 1.13 2014/10/13 08:52:53 j_novak Exp $" ;
31 
32 /*
33  * $Id: eos_mag.C,v 1.13 2014/10/13 08:52:53 j_novak Exp $
34  * $Log: eos_mag.C,v $
35  * Revision 1.13 2014/10/13 08:52:53 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.12 2014/09/22 16:13:10 j_novak
39  * Minor modif.
40  *
41  * Revision 1.11 2014/05/27 12:32:28 j_novak
42  * Added possibility to converge to a given magnetic moment.
43  *
44  * Revision 1.10 2014/05/13 15:37:12 j_novak
45  * Updated to new magnetic units.
46  *
47  * Revision 1.9 2014/04/28 12:48:13 j_novak
48  * Minor modifications.
49  *
50  * Revision 1.8 2014/03/11 14:27:26 j_novak
51  * Corrected a missing 4pi term.
52  *
53  * Revision 1.7 2014/03/03 16:23:08 j_novak
54  * Updated error message
55  *
56  * Revision 1.6 2013/12/12 16:07:30 j_novak
57  * interpol_herm_2d outputs df/dx, used to get the magnetization.
58  *
59  * Revision 1.5 2013/11/25 15:00:52 j_novak
60  * Correction of memory error.
61  *
62  * Revision 1.4 2013/11/14 16:12:55 j_novak
63  * Corrected a mistake in the units.
64  *
65  * Revision 1.2 2011/10/04 16:05:19 j_novak
66  * Update of Eos_mag class. Suppression of loge, re-definition of the derivatives
67  * and use of interpol_herm_2d.
68  *
69  * Revision 1.1 2011/06/16 10:49:18 j_novak
70  * New class Eos_mag for EOSs depending on density and magnetic field.
71  *
72  *
73  * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_mag.C,v 1.13 2014/10/13 08:52:53 j_novak Exp $
74  *
75  */
76 
77 // headers C
78 #include <cmath>
79 
80 // Headers Lorene
81 #include "eos.h"
82 #include "cmp.h"
83 #include "param.h"
84 #include "utilitaires.h"
85 #include "unites.h"
86 
87 
88 namespace Lorene {
89 void interpol_herm_2d(const Tbl&, const Tbl&, const Tbl&, const Tbl&, const Tbl&, const Tbl&, double, double, double&, double&, double&) ;
90 
91 
92  //----------------------------//
93  // Constructors //
94  //----------------------------//
95 
96 // Standard constructor
97 // --------------------
98 Eos_mag::Eos_mag(const char* name_i, const char* table,
99  const char* path) : Eos(name_i), tablename(path) {
100 
101  tablename += '/' ;
102  tablename += table ;
103 
104  read_table() ;
105 
106 }
107 
108 // Standard constructor with full filename
109 // ---------------------------------------
110 Eos_mag::Eos_mag(const char* name_i, const char* file_name)
111  : Eos(name_i), tablename(file_name) {
112 
113  read_table() ;
114 
115 }
116 
117 
118 // Constructor from binary file
119 // ----------------------------
120 Eos_mag::Eos_mag(FILE* fich) : Eos(fich) {
121 
122  char tmp_name[160] ;
123 
124  fread(tmp_name, sizeof(char), 160, fich) ;
125  tablename += tmp_name ;
126 
127  read_table() ;
128 
129 }
130 
131 
132 
133 // Constructor from a formatted file
134 // ---------------------------------
135 Eos_mag::Eos_mag(ifstream& fich, const char* table) : Eos(fich) {
136 
137  fich >> tablename ;
138  tablename += '/' ;
139  tablename += table ;
140 
141  read_table() ;
142 
143 }
144 
145 Eos_mag::Eos_mag(ifstream& fich) : Eos(fich) {
146 
147  fich >> tablename ;
148 
149  read_table() ;
150 
151 }
152 
153 
154  //--------------//
155  // Destructor //
156  //--------------//
157 
159  delete d2lp ;
160  delete dlpsdB ;
161  delete dlpsdlh ;
162  delete Bfield ;
163  delete logh ;
164  delete logp ;
165 }
166 
167  //------------//
168  // Outputs //
169  //------------//
170 
171 void Eos_mag::sauve(FILE* fich) const {
172 
173  Eos::sauve(fich) ;
174 
175  fwrite(tablename.c_str(), sizeof(char), tablename.size(), fich) ;
176 
177 }
178  //------------------------//
179  // Comparison operators //
180  //------------------------//
181 
182 
183 bool Eos_mag::operator==(const Eos& eos_i) const {
184 
185  bool resu = true ;
186 
187  if ( eos_i.identify() != identify() ) {
188  cerr << "The second EOS is not of type Eos_mag !" << endl ;
189  resu = false ;
190  }
191 
192  return resu ;
193 
194 }
195 
196 bool Eos_mag::operator!=(const Eos& eos_i) const {
197 
198  return !(operator==(eos_i)) ;
199 
200 }
201 
202  //------------//
203  // Outputs //
204  //------------//
205 
206 
207 ostream& Eos_mag::operator>>(ostream & ost) const {
208 
209  ost <<
210  "EOS of class Eos_mag : tabulated EOS depending on two variables: enthalpy and magnetic field."
211  << '\n' ;
212  ost << "Table read from " << tablename << endl ;
213 
214  return ost ;
215 
216 }
217 
218 
219 
220  //------------------------//
221  // Reading of the table //
222  //------------------------//
223 
225 
226  using namespace Unites_mag ;
227 
228  ifstream fich(tablename.data()) ;
229 
230  if (!fich) {
231  cerr << "Eos_mag::read_table(): " << endl ;
232  cerr << "Problem in opening the EOS file!" << endl ;
233  cerr << "Aborting..." << endl ;
234  abort() ;
235  }
236 
237  for (int i=0; i<5; i++) { // jump over the file
238  fich.ignore(1000, '\n') ; // header
239  } //
240 
241  int nbp1, nbp2 ;
242  fich >> nbp1 >> nbp2 ; fich.ignore(1000, '\n') ; // number of data
243  if ((nbp1<=0) || (nbp2<=0)) {
244  cerr << "Eos_mag::read_table(): " << endl ;
245  cerr << "Wrong value for the number of lines!" << endl ;
246  cerr << "nbp1 = " << nbp1 << ", nbp2 = " << nbp2 << endl ;
247  cerr << "Aborting..." << endl ;
248  abort() ;
249  }
250 
251  for (int i=0; i<3; i++) { // jump over the table
252  fich.ignore(1000, '\n') ;
253  }
254 
255  logp = new Tbl(nbp2, nbp1) ;
256  logh = new Tbl(nbp2, nbp1) ;
257  Bfield = new Tbl(nbp2, nbp1) ;
258  dlpsdlh = new Tbl(nbp2, nbp1) ;
259  dlpsdB = new Tbl(nbp2, nbp1) ;
260  d2lp = new Tbl(nbp2, nbp1) ;
261 
262  logp->set_etat_qcq() ;
263  logh->set_etat_qcq() ;
264  Bfield->set_etat_qcq() ;
265  dlpsdlh->set_etat_qcq() ;
266  dlpsdB->set_etat_qcq() ;
267  d2lp->set_etat_qcq() ;
268 
269  double c2 = c_si * c_si ;
270  double B_unit = mag_unit / 1.e11 ;
271  double M_unit = B_unit*mu0/(4*M_PI) ;
272 
273  int no1 ;
274  double nb_fm3, rho_cgs, p_cgs, mu_MeV, magB_PG, magM_PG, chi_PGpMeV ;
275 
276  double ww = 0. ;
277 
278  for (int j=0; j<nbp2; j++) {
279 
280  for (int i=0; i<nbp1; i++) {
281  fich >> no1 >> nb_fm3 >> rho_cgs >> p_cgs >> mu_MeV
282  >> magB_PG >> magM_PG >> chi_PGpMeV ;
283  fich.ignore(1000,'\n') ;
284  if ( (nb_fm3<0) || (rho_cgs<0)) { // || (p_cgs < 0) ){
285  cerr << "Eos_mag::read_table(): " << endl ;
286  cerr << "Negative value in table!" << endl ;
287  cerr << "nb = " << nb_fm3 << ", rho = " << rho_cgs <<
288  ", p = " << p_cgs << endl ;
289  cerr << "Aborting..." << endl ;
290  abort() ;
291  }
292 
293  double psc2 = 0.1*p_cgs/c2 ; // in kg m^-3
294  double rho_si = rho_cgs*1000. ; // in kg m^-3
295 
296  double h_read = log(mu_MeV) ;
297  if ( (i==0) && (j==0) ) ww = h_read ;
298  double h_new = h_read - ww ;
299 
300  logp->set(j, i) = psc2/rhonuc_si ;
301  logh->set(j, i) = h_new ;
302  Bfield->set(j, i) = magB_PG / B_unit ; // in Lorene units
303  dlpsdlh->set(j, i) = (rho_si + psc2)/rhonuc_si ;
304  dlpsdB->set(j, i) = magM_PG / M_unit ;
305  d2lp->set(j, i) = mu_MeV*chi_PGpMeV / (M_unit) ;
306 
307  }
308  }
309 
310  hmin = (*logh)(0, 0) ;
311  hmax = (*logh)(0, nbp1-1) ;
312 
313  Bmax = (*Bfield)(nbp2-1, 0) ;
314 
315  // cout << "hmin: " << hmin << ", hmax: " << hmax << endl ;
316  // cout << "Bmax: " << Bmax << endl ;
317 
318  fich.close();
319 
320 }
321 
322 
323  //------------------------------//
324  // Computational routines //
325  //------------------------------//
326 
327 // Baryon density from enthalpy
328 //------------------------------
329 
330 double Eos_mag::nbar_ent_p(double ent, const Param* par ) const {
331 
332  using namespace Unites_mag ;
333 
334  if ((ent > hmin - 1.e-12) && (ent < hmin))
335  ent = hmin ;
336 
337  if ( ent >= hmin) {
338  if (ent > hmax) {
339  cerr << "Eos_tabul::nbar_ent_p : ent > hmax !" << endl ;
340  abort() ;
341  }
342  // recuperer magB0 (input)
343  double magB0 = 0. ;
344  if (par != 0x0)
345  if (par->get_n_double_mod() > 0) {
346  magB0 = par->get_double_mod() ;
347  }
348 
349  if ( magB0 > Bmax) {
350  cerr << "Eos_tabul::nbar_ent_p : B > Bmax !" << endl ;
351  cerr << "B = " << magB0*mag_unit << ", Bmax = " << Bmax*mag_unit << endl ;
352  abort() ;
353  }
354 
355  double p_int, dpdb_int, dpdh_int ;
356 
357  interpol_herm_2d(*Bfield, *logh, *logp, *dlpsdB, *dlpsdlh, *d2lp, magB0, ent,
358  p_int, dpdb_int, dpdh_int) ;
359 
360  double nbar_int = dpdh_int * exp(-ent) ;
361 
362  return nbar_int ;
363  }
364  else{
365  return 0 ;
366  }
367 }
368 
369 
370 // Energy density from enthalpy
371 //------------------------------
372 
373 double Eos_mag::ener_ent_p(double ent, const Param* par ) const {
374 
375  using namespace Unites_mag ;
376 
377  if ((ent > hmin - 1.e-12) && (ent < hmin))
378  ent = hmin ;
379 
380  if ( ent >= hmin ) {
381  if (ent > hmax) {
382  cerr << "Eos_tabul::ener_ent_p : ent > hmax !" << endl ;
383  abort() ;
384  }
385 
386  // recuperer magB0 (input)
387  double magB0 = 0. ;
388  if (par != 0x0)
389  if (par->get_n_double_mod() > 0) {
390  magB0 = par->get_double_mod() ;
391  }
392 
393  if ( magB0 > Bmax) {
394  cerr << "Eos_tabul::ener_ent_p : B > Bmax !" << endl ;
395  cerr << "B = " << magB0*mag_unit << ", Bmax = " << Bmax*mag_unit << endl ;
396  abort() ;
397  }
398 
399 
400  double p_int, dpdb_int, dpdh_int ;
401 
402  interpol_herm_2d(*Bfield, *logh, *logp, *dlpsdB, *dlpsdlh, *d2lp, magB0, ent,
403  p_int, dpdb_int, dpdh_int) ;
404 
405  double nbar_int = dpdh_int * exp(-ent) ;
406 
407  double f_int = - p_int + exp(ent) * nbar_int;
408 
409  return f_int ;
410  }
411  else{
412  return 0 ;
413  }
414 }
415 
416 // Pressure from enthalpy
417 //------------------------
418 
419 double Eos_mag::press_ent_p(double ent, const Param* par ) const {
420 
421  using namespace Unites_mag ;
422 
423  if ((ent > hmin - 1.e-12) && (ent < hmin))
424  ent = hmin ;
425 
426  if ( ent >= hmin ) {
427  if (ent > hmax) {
428  cout << "Eos_mag::press_ent_p : ent > hmax !" << endl ;
429  abort() ;
430  }
431 
432  // recuperer magB0 (input)
433  double magB0 = 0. ;
434  if (par != 0x0)
435  if (par->get_n_double_mod() > 0) {
436  magB0 = par->get_double_mod() ;
437  }
438 
439  if ( magB0 > Bmax) {
440  cerr << "Eos_tabul::press_ent_p : B > Bmax !" << endl ;
441  cerr << "B = " << magB0*mag_unit << ", Bmax = " << Bmax*mag_unit << endl ;
442  abort() ;
443  }
444 
445  double p_int, dpdb_int, dpdh_int ;
446 
447  interpol_herm_2d(*Bfield, *logh, *logp, *dlpsdB, *dlpsdlh, *d2lp, magB0, ent,
448  p_int, dpdb_int, dpdh_int) ;
449 
450  return p_int;
451  }
452  else{
453  return 0 ;
454  }
455 }
456 
457 // Magnetization from enthalpy
458 //----------------------------
459 
460 double Eos_mag::mag_ent_p(double ent, const Param* par) const {
461 
462  using namespace Unites_mag ;
463 
464  if ((ent > hmin - 1.e-12) && (ent < hmin))
465  ent = hmin ;
466 
467  if ( ent >= hmin ) {
468  if (ent > hmax) {
469  cout << "Eos_mag::mag_ent_p : ent > hmax !" << endl ;
470  abort() ;
471  }
472 
473  // recuperer magB0 (input)
474  double magB0 = 0. ;
475  if (par != 0x0)
476  if (par->get_n_double_mod() > 0) {
477  magB0 = par->get_double_mod() ;
478  }
479 
480  if ( magB0 > Bmax) {
481  cerr << "Eos_tabul::mag_ent_p : B > Bmax !" << endl ;
482  cerr << "B = " << magB0*mag_unit << ", Bmax = " << Bmax*mag_unit << endl ;
483  abort() ;
484  }
485 
486  double p_int, dpdb_int, dpdh_int ;
487 
488  interpol_herm_2d(*Bfield, *logh, *logp, *dlpsdB, *dlpsdlh, *d2lp, magB0, ent,
489  p_int, dpdb_int, dpdh_int) ;
490 
491  double magnetization = dpdb_int ;
492 
493  if (magB0 == 0.)
494  return 0 ;
495  else
496  return mu0*magnetization / magB0 ;
497  }
498 
499  else
500  return 0. ;
501 
502 }
503 
504 
505 // dln(n)/ln(H) from enthalpy
506 //---------------------------
507 
508 double Eos_mag::der_nbar_ent_p(double ent, const Param* ) const {
509 
510  c_est_pas_fait("Eos_mag::der_nbar_ent_p" ) ;
511 
512  return ent ;
513 
514 }
515 
516 
517 // dln(e)/ln(H) from enthalpy
518 //---------------------------
519 
520 double Eos_mag::der_ener_ent_p(double ent, const Param* ) const {
521 
522 
523  c_est_pas_fait("Eos_mag::der_ener_enr_p" ) ;
524 
525  return ent ;
526 
527 }
528 
529 
530 // dln(p)/ln(H) from enthalpy
531 //---------------------------
532 
533 double Eos_mag::der_press_ent_p(double ent, const Param* ) const {
534 
535  c_est_pas_fait("Eos_mag::" ) ;
536 
537  return ent ;
538 
539 }
540 
541 
542 // dln(p)/dln(n) from enthalpy
543 //---------------------------
544 
545 double Eos_mag::der_press_nbar_p(double ent, const Param*) const {
546 
547  c_est_pas_fait("Eos_mag::der_press_nbar_p" ) ;
548 
549  return ent ;
550 
551 }
552 }
Tbl * dlpsdB
Table of .
Definition: eos_mag.h:109
double mag_ent_p(double ent, const Param *par=0x0) const
Computes the magnetisation.
Definition: eos_mag.C:460
Tbl * d2lp
Table of .
Definition: eos_mag.h:112
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
Definition: eos_mag.C:196
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
Tbl * logp
Table of .
Definition: eos_mag.h:97
Tbl * logh
Table of .
Definition: eos_mag.h:100
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_mag.C:533
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
Definition: eos_mag.C:419
virtual void sauve(FILE *) const
Save in a file.
Definition: eos_mag.C:171
double Bmax
Upper boundary of the magnetic field interval.
Definition: eos_mag.h:94
virtual ~Eos_mag()
Destructor.
Definition: eos_mag.C:158
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
Definition: eos_mag.C:330
Tbl * Bfield
Table of .
Definition: eos_mag.h:103
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_mag.C:508
virtual bool operator==(const Eos &) const
Comparison operator (egality)
Definition: eos_mag.C:183
double hmin
Lower boundary of the log-enthalpy interval.
Definition: eos_mag.h:88
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_mag.C:520
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
Definition: eos_mag.C:373
string tablename
Name of the file containing the tabulated data.
Definition: eos_mag.h:85
void read_table()
Reads the file containing the table and initializes in the arrays logh , logp and dlpsdlh .
Definition: eos_mag.C:224
Tbl * dlpsdlh
Table of .
Definition: eos_mag.h:106
virtual double der_press_nbar_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_mag.C:545
Eos_mag(const char *name_i, const char *table, const char *path)
Standard constructor.
Definition: eos_mag.C:98
double hmax
Upper boundary of the log-enthalpy interval.
Definition: eos_mag.h:91
virtual ostream & operator>>(ostream &) const
Operator >>
Definition: eos_mag.C:207
Equation of state base class.
Definition: eos.h:190
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual void sauve(FILE *) const
Save in a file.
Definition: eos.C:179
Parameter storage.
Definition: param.h:125
int get_n_double_mod() const
Returns the number of stored modifiable double 's addresses.
Definition: param.C:446
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
Definition: param.C:498
Basic array class.
Definition: tbl.h:161
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Lorene prototypes.
Definition: app_hor.h:64
Standard electro-magnetic units.