Improve ESTER chemistry fixing freeeos implementation

python/ester_wrap.cxx

@@ -14844,37 +14844,6 @@ SWIGINTERN PyObject *_wrap_atm_onelayer(PyObject *SWIGUNUSEDPARM(self), PyObject
 }
 
 
-SWIGINTERN PyObject *_wrap_initial_composition(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  double arg1 ;
-  double arg2 ;
-  double val1 ;
-  int ecode1 = 0 ;
-  double val2 ;
-  int ecode2 = 0 ;
-  PyObject * obj0 = 0 ;
-  PyObject * obj1 = 0 ;
-  double_map result;
-
-  if (!PyArg_ParseTuple(args,(char *)"OO:initial_composition",&obj0,&obj1)) SWIG_fail;
-  ecode1 = SWIG_AsVal_double(obj0, &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "initial_composition" "', argument " "1"" of type '" "double""'");
-  } 
-  arg1 = static_cast< double >(val1);
-  ecode2 = SWIG_AsVal_double(obj1, &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "initial_composition" "', argument " "2"" of type '" "double""'");
-  } 
-  arg2 = static_cast< double >(val2);
-  result = initial_composition(arg1,arg2);
-  resultobj = SWIG_NewPointerObj((new double_map(static_cast< const double_map& >(result))), SWIGTYPE_p_double_map, SWIG_POINTER_OWN |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_mapping_gl_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   mapping *arg1 = (mapping *) 0 ;
@@ -31325,7 +31294,6 @@ static PyMethodDef SwigMethods[] = {
 	 { (char *)"eos_idealrad", _wrap_eos_idealrad, METH_VARARGS, NULL},
 	 { (char *)"eos_opal", _wrap_eos_opal, METH_VARARGS, NULL},
 	 { (char *)"atm_onelayer", _wrap_atm_onelayer, METH_VARARGS, NULL},
-	 { (char *)"initial_composition", _wrap_initial_composition, METH_VARARGS, NULL},
 	 { (char *)"mapping_gl_set", _wrap_mapping_gl_set, METH_VARARGS, NULL},
 	 { (char *)"mapping_gl_get", _wrap_mapping_gl_get, METH_VARARGS, NULL},
 	 { (char *)"mapping_leg_set", _wrap_mapping_leg_set, METH_VARARGS, NULL},

python/ester_wrap.py

@@ -902,9 +902,6 @@ def atm_onelayer(X, Z, g, Teff, eos_name, opa_name, atm):
     return _ester_wrap.atm_onelayer(X, Z, g, Teff, eos_name, opa_name, atm)
 atm_onelayer = _ester_wrap.atm_onelayer
 
-def initial_composition(X, Z):
-    return _ester_wrap.initial_composition(X, Z)
-initial_composition = _ester_wrap.initial_composition
 MAP_BONAZZOLA = _ester_wrap.MAP_BONAZZOLA
 MAP_LINEAR = _ester_wrap.MAP_LINEAR
 class mapping(_object):

src/include/matrix.h

@@ -289,7 +289,9 @@ class matrix_map : public std::map<std::string,matrix> {
 	inline matrix &operator[](const std::string& key) {return std::map<std::string,matrix>::operator[](key);}
 	inline const matrix &operator[](const std::string& key) const {return find(key)->second;};
 	matrix_map_elem operator()(int nfil, int ncol);
+	matrix_map_elem operator()(int ielem);
 	const double_map operator()(int nfil, int ncol) const;
+	const double_map operator()(int ielem) const;
 	const matrix_map row(int irow) const;
 	const matrix_map col(int icol) const;
 	const matrix_map block(int irow1,int irow2,int icol1,int icol2) const;

src/include/physics.h

@@ -34,12 +34,12 @@ class composition_map : public matrix_map {
 
 int opa_calc(const matrix &X,double Z,const matrix &T,const matrix &rho,
 		opa_struct &opa);
-int eos_calc(const matrix &X,double Z,const matrix &T,const matrix &p,
-		matrix &rho,eos_struct &eos);
-int nuc_calc(const matrix_map &X,const matrix &T,const matrix &rho,
+int eos_calc(const composition_map &, const matrix &T, const matrix &p,
+		matrix &rho, eos_struct &eos);
+int nuc_calc(const composition_map &, const matrix &T, const matrix &rho,
 		nuc_struct &nuc);
-int atm_calc(const matrix &X,double Z,const matrix &g,const matrix &Teff,
-		const char *eos_name,const char *opa_name,atm_struct &atm);
+int atm_calc(const composition_map &, const matrix &g, const matrix &Teff,
+		const char *eos_name, const char *opa_name, atm_struct &atm);
 
 int opa_opal(const matrix &X,double Z,const matrix &T,const matrix &rho,
 		opa_struct &opa);
@@ -58,19 +58,17 @@ int nuc_cesam(const composition_map &comp,const matrix &T,const matrix &rho,
 int nuc_cesam_dcomp(composition_map &comp,const matrix &T,const matrix &rho,
 		nuc_struct &nuc);
 
-int eos_ideal(const matrix &X,double Z,const matrix &T,const matrix &p,
-		matrix &rho,eos_struct &eos);
-int eos_idealrad(const matrix &X,double Z,const matrix &T,const matrix &p,
-		matrix &rho,eos_struct &eos);
-int eos_opal(const matrix &X,double Z,const matrix &T,const matrix &p,
-		matrix &rho,eos_struct &eos);
-int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p,
+int eos_ideal(const composition_map &, const matrix &T, const matrix &p,
+		matrix &rho, eos_struct &eos);
+int eos_idealrad(const composition_map &, const matrix &T, const matrix &p,
+		matrix &rho, eos_struct &eos);
+int eos_opal(const composition_map &, const matrix &T, const matrix &p,
+		matrix &rho, eos_struct &eos);
+int eos_freeeos(const composition_map &, const matrix &T, const matrix &p,
 		matrix &rho, eos_struct &eos);
 
-int atm_onelayer(const matrix &X,double Z,const matrix &g,const matrix &Teff,
-		const char *eos_name,const char *opa_name,atm_struct &atm);
-
-double_map initial_composition(double X,double Z);
+int atm_onelayer(const composition_map &, const matrix &g, const matrix &Teff,
+		const char *eos_name, const char *opa_name, atm_struct &atm);
 
 #endif
 

src/include/star.h

@@ -37,6 +37,9 @@ class star2d {
         matrix rho, phi, p, T;
         matrix phiex;
         matrix vr, vt, G, w;
+        // the m_ prefix is used to explicit that it's a class' member https://stackoverflow.com/a/13018190/22158934 :
+        double m_He_isotopic_ratio; // He3/He4 hardcoded in star*d::init_comp // TODO: change it to be configurable
+        double_map m_metal_mix;
         composition_map comp;
         double X0, Y0, Z0;
         double R, M;
@@ -97,6 +100,7 @@ class star2d {
 
         virtual void dump_info();
 
+        virtual void init_metal_mix();
         virtual void init_comp();
 
         virtual solver *init_solver(int nvar_add=0);

src/matrix/matrix_map.cpp

@@ -51,6 +51,17 @@ matrix_map_elem matrix_map::operator()(int nfil, int ncol) {
 	return elem;
 }
 
+matrix_map_elem matrix_map::operator()(int ielem) {
+
+	matrix_map_elem elem;
+	matrix_map::iterator it;
+	for(it=begin();it!=end();it++) {
+		elem[it->first]=&(it->second)(ielem);
+	}
+
+	return elem;
+}
+
 const double_map matrix_map::operator()(int nfil, int ncol) const {
 
 	double_map elem;
@@ -62,6 +73,17 @@ const double_map matrix_map::operator()(int nfil, int ncol) const {
 	return elem;
 }
 
+const double_map matrix_map::operator()(int ielem) const {
+
+	double_map elem;
+	matrix_map::const_iterator it;
+	for(it=begin();it!=end();it++) {
+		elem[it->first]=(it->second)(ielem);
+	}
+
+	return elem;
+}
+
 const matrix_map matrix_map::row(int irow) const {
 
 	matrix_map res;

src/physics/EOS5_xtrin.F

@@ -121,7 +121,7 @@ subroutine esac (xh,ztab,t6,r,iorder,irad)
 c..... read the data files
         call readcoeos
         z=zz(1)
-        if (ztab .ne. z) then
+        if (abs(ztab-z) .ge. 1e-7) then
         write (*,'("requested z=",f10.6," EOS5_data is for z=",
      x             f10.6)')
      x    ztab,z      

src/physics/atm_onelayer.cpp

@@ -7,8 +7,8 @@
 #include <string.h>
 #include <cmath>
 
-int atm_onelayer(const matrix &X,double Z,const matrix &g,const matrix &Teff,
-		const char *eos_name,const char *opa_name,atm_struct &atm) {
+int atm_onelayer(const composition_map &chemical_comp, const matrix &g, const matrix &Teff,
+		const char *eos_name, const char *opa_name, atm_struct &atm) {
 
 	int n=g.nrows();
 	int m=g.ncols();
@@ -43,8 +43,8 @@ int atm_onelayer(const matrix &X,double Z,const matrix &g,const matrix &Teff,
 				}
 				double F,dF,dlogps;
 				p=pow(10,logps)*ones(1,1);
-				if(eos_calc(X(i,j)*ones(1,1),Z,T,p,rho,eos)) return 1;
-				if(opa_calc(X(i,j)*ones(1,1),Z,T,rho,opa)) return 1;
+				if(eos_calc(chemical_comp(i,j)*ones(1,1), T, p, rho, eos)) return 1;
+				if(opa_calc(chemical_comp.X()(i,j)*ones(1,1), chemical_comp.Z()(i,j), T, rho, opa)) return 1;
 				F=logps+log10(opa.k)(0)-logg-log10(2./3.);
 				dF=1.-(1+opa.dlnxi_lnrho(0))/eos.chi_rho(0);
 				dlogps=-F/dF;

src/physics/composition.cpp

@@ -1,30 +1,8 @@
 #ifndef WITH_CMAKE
 #include "ester-config.h"
 #endif
-#include "physics.h"
-
-double_map initial_composition(double X, double Z) {
-	double_map comp;
-
-// The mixture below has an unidentified origin but is close to the one
-// of Grevesse & Sauval 1998.
-	comp["H"]=X;
-	comp["He3"]=3.15247417638132e-04*(1.-X-Z);
-	comp["He4"]=(1.-X-Z)-comp["He3"];
-	comp["C12"]=Z*1.71243418737847e-01;
-	comp["C13"]=Z*2.06388003380057e-03;
-	comp["N14"]=Z*5.29501630871695e-02;
-	comp["N15"]=Z*2.08372414940812e-04;
-	comp["O16"]=Z*4.82006487350336e-01;
-	comp["O17"]=Z*1.95126448826986e-04;
-
-	double tot=comp.sum();
-
-	comp["Ex"] =1-tot;
-
-	return comp;
 
-}
+#include "physics.h"
 
 matrix composition_map::X() const {
 

src/physics/eos_freeeos.cpp

@@ -21,7 +21,33 @@ extern "C" {
             double *entropy, int *iteration_count);
 }
 
-int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p,
+void set_eps(double* eps, const double_map &local_chemical_mix) {
+        // TODO: maybe use atomic masses defined in src/global/global.cpp
+        // Only isotops that exists in non negligible quantities have been considered
+        // according to their wikipedia page on the 3rd of July 2023
+        eps[0] = local_chemical_mix["H"] / 1.008e0;                     // H
+        eps[1] = local_chemical_mix["He3"] / 3. + local_chemical_mix["He4"] / 4.;     // He3 + He4
+        eps[2] = local_chemical_mix["C12"] / 12. + local_chemical_mix["C13"] / 13.;   // C12 + C13
+        eps[3] = local_chemical_mix["N14"] / 14. + local_chemical_mix["N15"] / 15.;   // N14 + N15
+        eps[4] = local_chemical_mix["O16"] / 16. + local_chemical_mix["O17"] / 17.;   // O16 + O17
+        eps[5] = local_chemical_mix["Ne20"] / 20. + local_chemical_mix["Ne22"] / 22.; // Ne20 + Ne22
+        eps[6] = local_chemical_mix["Na23"] / 23.; // Na23
+        eps[7] = local_chemical_mix["Mg24"] / 24 + local_chemical_mix["Mg25"] / 25 + local_chemical_mix["Mg26"] / 26; // Mg24 + Mg25 + Mg26
+        eps[8] = local_chemical_mix["Al27"] / 27.0; // Al27
+        eps[9] = local_chemical_mix["Si28"] / 28.0; // Si28
+        eps[10] = local_chemical_mix["P31"] / 31.0; // P31
+        eps[11] = local_chemical_mix["S32"] / 32.0; // S32
+        eps[12] = local_chemical_mix["Cl35"] / 35.0 + local_chemical_mix["Cl37"] / 37.0; // Cl35 + Cl37
+        eps[13] = local_chemical_mix["A40"] / 40.0; // A40
+        eps[14] = local_chemical_mix["Ca40"] / 40.0; // Ca40
+        eps[15] = local_chemical_mix["Ti"] / 47.867; // Ti //there are too many isotops
+        eps[16] = local_chemical_mix["Cr"] / 52.0; // Cr //there are too many isotops
+        eps[17] = local_chemical_mix["Mn55"] / 55.0; // Mn
+        eps[18] = local_chemical_mix["Fe"] / 55.845; // Fe //there are too many isotops
+        eps[19] = local_chemical_mix["Ni"] / 58.693; // Ni //there are too many isotops
+}
+
+int eos_freeeos(const composition_map &chemical_comp, const matrix &T, const matrix &p,
         matrix &rho, eos_struct &eos) {
 
     double t;
@@ -84,29 +110,7 @@ int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p,
     eos.prad.dim(T.nrows(), T.ncols()); // added MR june 2023
 
     for (int i=0; i<N; i++) {
-
-        double_map comp = initial_composition(X(i), Z);
-
-        eps[0] = comp["H"] / 1.008e0;                       // H
-        eps[1] = (comp["He3"] + comp["He4"]) / 4.0026e0;   // He3 + He4
-        eps[2] = (comp["C12"] + comp["C13"]) / 12.0111e0; // C12 + C13
-        eps[3] = (comp["N14"] + comp["N15"]) / 14.0067e0; // N14 + N15
-        eps[4] = (comp["O16"] + comp["O17"]) / 15.9994e0; // O16 + O17
-        eps[5] = 0.0; // Ne
-        eps[6] = 0.0; // Na
-        eps[7] = 0.0; // Mg
-        eps[8] = 0.0; // AL
-        eps[9] = 0.0; // Si
-        eps[10] = .0; // P
-        eps[11] = .0; // S
-        eps[12] = .0; // Cl
-        eps[13] = .0; // A
-        eps[14] = .0; // Ca
-        eps[15] = .0; // Ti
-        eps[16] = .0; // Cr
-        eps[17] = .0; // Mn
-        eps[18] = .0; // Fe
-        eps[19] = .0; // Ni
+        set_eps(eps, chemical_comp(i));
 
         double pi = p(i);
         double match_variable = log(pi);
@@ -123,12 +127,13 @@ int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p,
                 &iteration_count);
         rho(i) = rhoi;
         if (iteration_count < 0) {
+            // TODO: this makes no sense freeEOS has no tables...
             ester_err(
                     "Values outside freeEOS eos table:\n"
                     "  X = %e\n"
                     "  Z = %e\n"
                     "  T = %e\n"
-                    "  p = %e", X(i), Z, t, p(i));
+                    "  p = %e", chemical_comp.X()(i), chemical_comp.Z()(i), t, p(i));
         }
         eos.s(i) = entropy[0];
         eos.G1(i) = gamma1;

src/physics/eos_ideal.cpp

@@ -5,12 +5,13 @@
 #include "physics.h"
 #include "constants.h"
 
-int eos_ideal(const matrix &X,double Z,const matrix &T,const matrix &p,
-		matrix &rho,eos_struct &eos) {
+int eos_ideal(const composition_map &chemical_comp, const matrix &T, const matrix &p,
+		matrix &rho, eos_struct &eos) {
 
 	matrix mu,b;
-
-	mu=4/(3+5*X-Z);
+
+	// TODO: this mix between X and Z implicitly assume a certain Z mixture
+	mu = 4/(3 + 5*chemical_comp.X() - chemical_comp.Z());
 	eos.prad=zeros(T.nrows(),T.ncols());
 	rho=p/T/K_BOL*mu*HYDROGEN_MASS;
 	b=ones(T.nrows(),T.ncols());

src/physics/eos_idealrad.cpp

@@ -4,12 +4,13 @@
 #include "physics.h"
 #include "constants.h"
 
-int eos_idealrad(const matrix &X,double Z,const matrix &T,const matrix &p,
-		matrix &rho,eos_struct &eos) {
+int eos_idealrad(const composition_map &chemical_comp, const matrix &T, const matrix &p,
+		matrix &rho, eos_struct &eos) {
 
 	matrix mu,b;
 
-	mu=4/(3+5*X-Z);
+	// TODO: this mix between X and Z implicitly assume a certain Z mixture
+	mu = 4/(3 + 5*chemical_comp.X() - chemical_comp.Z());
 	eos.prad=A_RAD/3*pow(T,4);
 	rho=(p-eos.prad)/T/K_BOL*mu*HYDROGEN_MASS;
 	b=(p-eos.prad)/p;