LORENE
eos_compose_fit_build.C
1 /*
2  * Methods of class Eos_compose_fit building.
3  *
4  * (see file eos_compose_fit.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2022 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 
31 
32 /*
33  * $Id: eos_compose_fit_build.C,v 1.6 2023/07/12 15:39:30 j_novak Exp $
34  * $Log: eos_compose_fit_build.C,v $
35  * Revision 1.6 2023/07/12 15:39:30 j_novak
36  * Reversed order of matching the EoS fit: it now starts from the crust (polytrope) and adjusts the kappa of this crust so as to have the right pressure at the highest density in the original table.
37  *
38  * Revision 1.4 2023/03/17 16:10:46 j_novak
39  * Better computation of gamma_low, to avoid jumps in sound speed.
40  *
41  * Revision 1.3 2023/01/27 16:10:35 j_novak
42  * A polytrope (class Eos_poly) is used for low and high enthalpies.
43  *
44  * Revision 1.2 2022/07/21 12:34:18 j_novak
45  * Corrected units
46  *
47  * Revision 1.1 2022/04/15 13:39:24 j_novak
48  * New class Eos_compose_fit to generate fitted EoSs from CompOSE tables.
49  *
50  *
51  * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_compose_fit_build.C,v 1.6 2023/07/12 15:39:30 j_novak Exp $
52  *
53  */
54 
55 // Headers Lorene
56 #include "eos_compose_fit.h"
57 #include "scalar.h"
58 #include "utilitaires.h"
59 #include "unites.h"
60 #include "graphique.h"
61 
62 namespace Lorene {
63  void Eos_compose_fit::read_compose_data(int& nbp, Tbl*& logh_read,
64  Tbl*& logp_read, Tbl*& loge_read,
65  Tbl*& lognb_read, Tbl*& gam1_read) {
66  // Files containing data and a test
67  //---------------------------------
68 
69  string file_nb = tablename + "eos.nb" ;
70  string file_thermo = tablename + "eos.thermo" ;
71 
72  cout << "opening " << file_nb << " and " << file_thermo << " ... " << flush ;
73  ifstream in_nb(file_nb.data()) ;
74  if (!in_nb) {
75  cerr << "Eos_compose_fit::read_compose_data:" << endl ;
76  cerr << "Problem in opening the EOS file!" << endl ;
77  cerr << "While trying to open " << file_nb << endl ;
78  cerr << "Aborting..." << endl ;
79  abort() ;
80  }
81 
82  // obtaining the size of the tables for memory allocation
83  //-------------------------------------------------------
84 
85  int index1, index2 ;
86  in_nb >> index1 >> index2 ;
87  nbp = index2 - index1 + 1 ;
88  assert(nbp > 0) ;
89 
90  Tbl press(nbp) ; press.set_etat_qcq() ;
91  Tbl nb(nbp) ; nb.set_etat_qcq() ;
92  Tbl ener(nbp) ; ener.set_etat_qcq() ;
93 
94  if (logh_read != nullptr) delete logh_read ;
95  logh_read = new Tbl(nbp) ;
96  logh_read->set_etat_qcq() ;
97  if (logp_read != nullptr) delete logp_read ;
98  logp_read = new Tbl(nbp) ;
99  logp_read->set_etat_qcq() ;
100  if (loge_read != nullptr) delete loge_read ;
101  loge_read = new Tbl(nbp) ;
102  loge_read->set_etat_qcq() ;
103  if (lognb_read != nullptr) delete lognb_read ;
104  lognb_read = new Tbl(nbp) ;
105  lognb_read->set_etat_qcq() ;
106  if (gam1_read != nullptr) delete gam1_read ;
107  gam1_read = new Tbl(nbp) ;
108  gam1_read->set_etat_qcq() ;
109 
110  // Variables and conversion
111  //-------------------------
112  double nb_fm3, rho_cgs, p_cgs, p_over_nb_comp, eps_comp ;
113  double dummy_x ;
114  int dummy_n ;
115 
116  //# Dealing with units could be optimized...
117  double rhonuc_cgs = Unites::rhonuc_si * 1e-3 ;
118  double c2_cgs = Unites::c_si * Unites::c_si * 1e4 ;
119  double m_neutron_MeV, m_proton_MeV ;
120 
121  ifstream in_p_rho (file_thermo.data()) ;
122  if (!in_p_rho) {
123  cerr << "Reading CompOSE data: " << endl ;
124  cerr << "Problem in opening the EOS file!" << endl ;
125  cerr << "While trying to open " << file_thermo << endl ;
126  cerr << "Aborting..." << endl ;
127  abort() ;
128  }
129  in_p_rho >> m_neutron_MeV >> m_proton_MeV ; //Neutron and proton masses
130  in_p_rho.ignore(1000, '\n') ;
131 
132  // Conversion from MeV/fm^3 to cgs
133  double p_convert = Unites::mev_si * 1.e45 * 10. ;
134  // From meV/fm^3 to g/cm^3
135  double eps_convert = Unites::mev_si * 1.e42 / (Unites::c_si*Unites::c_si) ;
136 
137  cout << " done. " << endl ;
138  cout << "Reading the table ... " << flush ;
139 
140  //--------------------------------------
141  // Main loop reading the CompOSE tables
142  //--------------------------------------
143  for (int i=0; i<nbp; i++) {
144  in_nb >> nb_fm3 ;
145  in_p_rho >> dummy_n >> dummy_n >> dummy_n >> p_over_nb_comp ;
146  in_p_rho >> dummy_x >> dummy_x >> dummy_x >> dummy_x >> dummy_x
147  >> eps_comp ;
148  in_p_rho.ignore(1000, '\n') ;
149  p_cgs = p_over_nb_comp * nb_fm3 * p_convert ;
150  rho_cgs = ( eps_comp + 1. ) * m_neutron_MeV * nb_fm3 * eps_convert ;
151 
152  if ( (nb_fm3 < 0) || (rho_cgs < 0) || (p_cgs < 0) ){
153  cerr << "Eos_compose_fit::read_compose_data: " << endl ;
154  cerr << "Negative value in table!" << endl ;
155  cerr << "nb = " << nb_fm3 << ", rho = " << rho_cgs <<
156  ", p = " << p_cgs << endl ;
157  cerr << "Aborting..." << endl ;
158  abort() ;
159  }
160 
161  press.set(i) = p_cgs / c2_cgs ;
162  nb.set(i) = nb_fm3 ;
163  ener.set(i) = rho_cgs ;
164 
165  // log-quantities in cgs units
166  logp_read->set(i) = log( press(i) / rhonuc_cgs ) ;
167  loge_read->set(i) = log( ener(i) / rhonuc_cgs ) ;
168  lognb_read->set(i) = log( 10.* nb(i) ) ;
169  }
170  nb_max = nb(nbp-1) ;
171  cout << " done." << endl ;
172  // -----------------------------------------
173  // End of the reading of the CompOSE tables
174  // -----------------------------------------
175 
176  cout << "Computation of derived quantities and file outputs.. " << flush ;
177 
178  // Computation of various derived quantities
179  //-------------------------------------------
180 
181  // log-enthallpy, with a shift ensuring h=10^{-14} for the lowest density
182  double ww = 0. ;
183  for (int i=0; i<nbp; i++) {
184  double h = log( (ener(i) + press(i)) / (10 * nb(i) * rhonuc_cgs) ) ;
185 
186  if (i==0) { ww = h ; }
187  h = h - ww + 1.e-14 ;
188 
189  logh_read->set(i) = log( h ) ;
190  }
191 
192  // d log(p) / d log(n) -> adiabatic index \Gamma_1
193  compute_derivative((*lognb_read), (*logp_read), (*gam1_read) ) ;
194 
195  cout << "done" << endl ;
196 
197  } // End (read_compose_data)
198 
199 
200  void Eos_compose_fit::adiabatic_index_fit(int i_min, int i_mid, const Tbl&
201  logh_read, const Tbl& logp_read,
202  const Tbl&, const Tbl&
203  lognb_read, const Tbl& gam1_read) {
204 
205  int nbp = logh_read.get_dim(0) ;
206  double nb0_max = exp(lognb_read(nbp-1)) ;
207 
208  assert((mp != nullptr) && (mg != nullptr)) ;
209  const Map_af& mpd = *mp ;
210  int nz = mg->get_nzone() ;
211  assert(nz == 2) ;
212  int nr = mg->get_nr(0) ;
213 #ifndef NDEBUG
214  for (int l=1; l<nz; l++)
215  assert(mg->get_nr(l) == nr) ;
216 #endif
217  const Coord& xx = mpd.r ; // xx = log(h) = log(log(e+p/m_B n_B))
218  double x_max = (+xx)(nz-1, 0, 0, nr-1) ;
219  double x_mid = (+xx)(1, 0, 0, 0) ;
220  double x_min = (+xx)(0, 0, 0, 0) ;
221 
222  if (log_p != nullptr) delete log_p ;
223  if (log_e != nullptr) delete log_e ;
224  if (log_nb != nullptr) delete log_nb ;
225  if (log_cs2 != nullptr) delete log_cs2 ;
226 
227  int np_gam = 0 ;
228  double mean_gam = 0. ;
229  for (int i=0; i<i_min; i++) {
230  mean_gam += gam1_read(i) ;
231  np_gam++ ;
232  }
233 
234  Scalar expexpx(mpd) ; // e^(e^x) , or simply (e+p)/(m_B n_B)
235  expexpx = exp(exp(xx)) ;
236  expexpx.std_spectral_base() ;
237  Scalar expx(mpd) ; // e^x or simply h = log((e+p)/(m_B n_B))
238  expx = exp(xx) ;
239  expx.std_spectral_base() ;
240 
241  mean_gam /= double(np_gam) ;
242  cout << "Mean Gamma_1 = " << mean_gam << endl ;
243 
244  // Fitting of the adiabatic index \Gamma_1
245  Scalar gam1_spec(mpd) ;
246  gam1_spec = mean_gam ;
247  gam1_spec.std_spectral_base() ;
248 
249  // High density part : polynomial fit (polynomial regression)
250  //--------------------------------------------------------------
251  int ndata = nbp - i_mid ;
252  Tbl xdata(ndata) ; xdata.set_etat_qcq() ;
253  Tbl ydata(ndata) ; ydata.set_etat_qcq() ;
254  for (int i=i_mid; i<nbp; i++) {
255  xdata.set(i-i_mid) = 2.*(logh_read(i) - x_mid)
256  / (x_max - x_mid) - 1. ;
257  ydata.set(i-i_mid) = gam1_read(i) ;
258  }
259 
260  Tbl tcoefs = poly_regression(xdata, ydata, n_coefs) ;
261 
262  // Putting coefficients to the 'Scalar' object
263  gam1_spec.set_spectral_va().coef() ;
264  for (int i=0; i<=n_coefs; i++) {
265  (*gam1_spec.set_spectral_va().c_cf).set(1, 0, 0, i) = tcoefs(i) ;
266  }
267 
268  // Not nice, but only simple way to update values on grid points
269  if (gam1_spec.set_spectral_va().c != nullptr) {
270  delete gam1_spec.set_spectral_va().c ;
271  gam1_spec.set_spectral_va().c = nullptr ;
272  }
273  gam1_spec.set_spectral_va().coef_i() ;
274 
275  // Intermediate part : linear interpolation
276  //------------------------------------------
277  Scalar interm(mpd) ;
278  // Linear formula
279  interm = mean_gam
280  + (gam1_spec.val_grid_point(1, 0, 0, 0) - mean_gam)
281  * (xx - x_min) / (x_mid - x_min) ;
282 
283  gam1_spec.set_domain(0) = interm.domain(0) ;
284 
285  // Low density part : polytrope, \Gamma_1 = const.
286  //-------------------------------------------------
287 
288  Scalar YYY(mpd) ;
289  // Determining kappa_low
290  double p_min = exp(logp_read(i_min)) ;
291  double ccc = (mean_gam - 1.) / mean_gam ;
292  double kappa_low0 = pow(ccc*expm1(exp(x_min)), mean_gam)
293  / pow(p_min, mean_gam*ccc) ;
294  double kappa_low1 = 0.1*kappa_low0 ;
295  double logp_max0 = integrate_equations(kappa_low0, mean_gam, gam1_spec, expx,
296  expexpx, YYY) ;
297  double logp_max = logp_read(nbp-1) ;
298 
299  // 'kappa_low' is determined to get the correct pressure at the higher end
300  //-------------------------------------------------------------------------
301  while(fabs(logp_max0 - logp_max)/logp_max > 1.e-3) {
302  double logp_max1 = integrate_equations(kappa_low1, mean_gam, gam1_spec, expx,
303  expexpx, YYY) ;
304  double kappa_low_new
305  = (kappa_low0*(logp_max1 -logp_max)- kappa_low1*(logp_max0 - logp_max))
306  / (logp_max1 - logp_max0) ;
307  cout << "Kappa(new) = " << kappa_low_new << ", logp = " << logp_max1 << endl ;
308  if (kappa_low_new < 0.) kappa_low_new = 0.01*kappa_low0 ;
309  kappa_low0 = kappa_low1 ;
310  logp_max0 = logp_max1 ;
311  kappa_low1 = kappa_low_new ;
312  }
313 
314  //des_profile(YYY, x_min, x_max, 0, 0) ;
315 
316  // sound speed : log(cs^2) = log(Y*\Gamma_1)
317  log_cs2 = new Scalar(log(YYY*gam1_spec)) ;
319 
320  // Pressure
321  Scalar spec_press = exp((*log_p)) ;
322  spec_press.std_spectral_base( );
323 
324  // Energy density
325  Scalar spec_ener = spec_press*(1./YYY - 1.) ; spec_ener.std_spectral_base() ;
326  log_e = new Scalar(log(spec_ener)) ;
328 
329  // Baryon density
330  Scalar spec_nbar = spec_press / (YYY * expexpx) ;
331  cout << "Relative difference in baryon density at the last point" << endl ;
332  cout << "(fitted computed / original tabulated) : "
333  << fabs(1. - spec_nbar.val_grid_point(nz-1, 0, 0, nr-1) / nb0_max)
334  << endl ;
335  log_nb = new Scalar(log(spec_nbar)) ;
336  log_nb->std_spectral_base() ; // log_nb_spec = log(n_B)
337 
341  spec_nbar.set_spectral_va().coef_i() ;
343 
344  // Matching to the high-density part polytrope
345  //--------------------------------------------
346 
347  double gam_high = gam1_spec.val_grid_point(1, 0, 0, nr-1) ;
348  double kappa_high = spec_press.val_grid_point(1, 0, 0, nr-1)
349  / pow(spec_nbar.val_grid_point(1, 0, 0, nr-1), gam_high) ;
350 
351  double mu_high = (spec_ener.val_grid_point(1, 0, 0, nr-1)
352  - kappa_high/(gam_high - 1.)
353  * pow(spec_nbar.val_grid_point(1, 0, 0, nr-1), gam_high))
354  / spec_nbar.val_grid_point(1, 0, 0, nr-1);
355 
356  if (p_eos_high != nullptr) delete p_eos_high ;
357  p_eos_high = new Eos_poly(gam_high, kappa_high, mu_high, mu_high) ;
358 
359  cout << "... done (analytic representation)" << endl ;
360 
361  } // End (adiabatic_index_fit)
362 
364  (double kappa_low, double mean_gam, const Scalar& gam1_spec,
365  const Scalar& expx, const Scalar& expexpx, Scalar& YYY) {
366 
367  // Enthalpy at lower end of the grid
368  double ent_min = expx.val_grid_point(0, 0, 0, 0) ;
369 
370  // The low-density part is defined through an "Eos_poly" object
371  if (p_eos_low != nullptr) delete p_eos_low ;
372  p_eos_low = new Eos_poly(mean_gam, kappa_low) ;
373 
374  // pressure & energy from th epolytropic EoS
375  double p_min = p_eos_low->press_ent_p(ent_min) ;
376  double e_min = p_eos_low->ener_ent_p(ent_min) ;
377 
378  // Solution of the differential equations to compute pressure, energy, etc
379  //--------------------------------------------------------------------------
380  // rhs: f = (\Gamma_1 - 1.)/\Gamma_1 * e^x * e^(e^x)
381  Scalar fff = ((gam1_spec-1.)/gam1_spec)*expx*expexpx ;
382  fff.std_spectral_base() ;
383 
384  // First integration
385  //-------------------
386  // F = primitive(f)
387  Scalar FFF = fff.primr() ;
388  FFF.std_spectral_base() ;
389  // Integration constant
390  double Y_min = p_min / (e_min + p_min) ;
391  double A_min = expexpx.val_grid_point(0, 0, 0, 0)*Y_min
392  - FFF.val_grid_point(0, 0, 0, 0) ;
393 
394  // Y is solution of Y' + Y e^x = f
395  YYY = (A_min + FFF)/expexpx ;
396 
397  // Second integration
398  //--------------------
399  // \Pi' = e^x / Y
400  Scalar Pidot = expx / YYY ;
401  if (log_p != nullptr) delete log_p ;
402  log_p = new Scalar(Pidot.primr()) ; // log_p_spec = log(p)
403  // Integration constant
404  double lnp0 = log(p_min) - log_p->val_grid_point(0, 0, 0, 0) ;
405  (*log_p) = (*log_p) + lnp0 ;
406 
407  // returning the log of pressure at the higher end of the grid
408  return log_p->val_grid_point(1, 0, 0, mg->get_nr(1)-1) ;
409  }
410 
411 }
void read_compose_data(int &nbp, Tbl *&logh, Tbl *&logp, Tbl *&loge, Tbl *&lognb, Tbl *&gam1)
Reads Compose data and stores the values of thermodynamic quantities.
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
Tbl poly_regression(const Tbl &, const Tbl &, int)
Polynomial regression, giving Chebyshev coefficients.
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
const Mg3d * mg
Multi-grid defining the number of domains and spectral coefficients.
const Tbl & domain(int l) const
Read-only of the value in a given domain.
Definition: scalar.h:631
Lorene prototypes.
Definition: app_hor.h:67
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
void coef_i() const
Computes the physical value of *this.
string tablename
Name of the file containing the tabulated data.
Scalar primr(bool null_infty=true) const
Computes the radial primitive which vanishes for .
Definition: scalar_integ.C:104
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition: scalar.h:621
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:643
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
void adiabatic_index_fit(int i_min, int i_max, const Tbl &logh_read, const Tbl &logp_read, const Tbl &loge_read, const Tbl &lognb_read, const Tbl &gam1_read)
From the read values, makes the fit on the adiabatic index and deduces the other quantities from ther...
Scalar * log_cs2
Table of .
double integrate_equations(double kappa, double mean_gam, const Scalar &gamma1, const Scalar &log_ent, const Scalar &ent, Scalar &YYY)
Integrates the differential system giving all thermo quantities from the adiabatic index...
Scalar * log_nb
Table of .
Scalar * log_e
Table of .
int n_coefs
Number of coeficients for polynomial regression.
const Map_af * mp
Mapping in .
double nb_max
Higher bound on the density, above which the EoS is assumed to be a polytrope.
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:309
Polytropic equation of state (relativistic case).
Definition: eos.h:812
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
void compute_derivative(const Tbl &xx, const Tbl &ff, Tbl &dfdx)
Derives a function defined on an unequally-spaced grid, approximating it by piecewise parabolae...
Definition: misc.C:64
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:423
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
const Eos_poly * p_eos_high
Pointer on a polytropic EoS for the description of high densities (nb>nb_max)
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Affine radial mapping.
Definition: map.h:2048
Scalar * log_p
Table of .
Basic array class.
Definition: tbl.h:164
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Coord r
r coordinate centered on the grid
Definition: map.h:736