From fe42ca7da193f3f81bc1843dcb4c4a3529650991 Mon Sep 17 00:00:00 2001 From: ctmbl Date: Thu, 29 Jun 2023 09:54:18 +0200 Subject: [PATCH 01/20] Update intial_composition to read an hardcoded metal-mix.cfg file --- src/physics/composition.cpp | 60 ++++++++++++++++++++++++++----------- 1 file changed, 43 insertions(+), 17 deletions(-) diff --git a/src/physics/composition.cpp b/src/physics/composition.cpp index 1deacba7..ec3574a8 100644 --- a/src/physics/composition.cpp +++ b/src/physics/composition.cpp @@ -2,28 +2,54 @@ #include "ester-config.h" #endif #include "physics.h" +#include "parser.h" + double_map initial_composition(double X, double Z) { double_map comp; + file_parser fp; -// 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; + // 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 Y = 1. - (X + Z); + comp["H"] = X; + comp["He3"] = 3.15247417638132e-04 * Y; + comp["He4"] = Y - comp["He3"]; + + // TODO: change this hardcoded + char file[] = "metal-mix.cfg"; + + char* arg = NULL; + char* val = NULL; + + if(!fp.open(file)){ + printf("Can't open configuration file %s\n", file); + perror("Error:"); + exit(1); + } else { + int line; + while(line = fp.get(arg,val)) { + if(val == NULL){ + printf("Syntax error in configuration file %s, line %d\n", file, line); + exit(1); + } + comp[arg] = Z * atof(val); + } + } + fp.close(); + comp["Ex"] = 1 - comp.sum(); + + return comp; } matrix composition_map::X() const { From 8276b8ceb9a9479dac6ec2f3eb9bcbcd6081c98b Mon Sep 17 00:00:00 2001 From: ctmbl Date: Mon, 3 Jul 2023 10:54:53 +0200 Subject: [PATCH 02/20] Remove old metal comp initialization (commented) --- src/physics/composition.cpp | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/physics/composition.cpp b/src/physics/composition.cpp index ec3574a8..9a0744a6 100644 --- a/src/physics/composition.cpp +++ b/src/physics/composition.cpp @@ -9,17 +9,6 @@ double_map initial_composition(double X, double Z) { double_map comp; file_parser fp; - // 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 Y = 1. - (X + Z); comp["H"] = X; From ff472e66abda99507d25924e17c2f0360cc307a2 Mon Sep 17 00:00:00 2001 From: ctmbl Date: Mon, 3 Jul 2023 10:58:03 +0200 Subject: [PATCH 03/20] Initialize comp object with all metals (isotopic or not) to 0 --- src/physics/composition.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/physics/composition.cpp b/src/physics/composition.cpp index 9a0744a6..8121ba5b 100644 --- a/src/physics/composition.cpp +++ b/src/physics/composition.cpp @@ -1,6 +1,8 @@ #ifndef WITH_CMAKE #include "ester-config.h" #endif +#include +#include #include "physics.h" #include "parser.h" @@ -20,6 +22,11 @@ double_map initial_composition(double X, double Z) { char* arg = NULL; char* val = NULL; + std::set metals = {"C12","C13","N14","N15","O16","O17","Ne","Na","Mg","AL","Si","P","S","Cl","A","Ca","Ti","Cr","Mn","Fe","Ni"}; + // Initialization of comp + for(std::string metal: metals){ + comp[metal] = .0; + } if(!fp.open(file)){ printf("Can't open configuration file %s\n", file); From 73f89d6343d32dfb73c59b40af2704b37972cd4f Mon Sep 17 00:00:00 2001 From: ctmbl Date: Mon, 3 Jul 2023 11:14:48 +0200 Subject: [PATCH 04/20] Add a TODO --- src/physics/composition.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/physics/composition.cpp b/src/physics/composition.cpp index 8121ba5b..63ea61e0 100644 --- a/src/physics/composition.cpp +++ b/src/physics/composition.cpp @@ -14,6 +14,7 @@ double_map initial_composition(double X, double Z) { double Y = 1. - (X + Z); comp["H"] = X; + // TODO: should this one be hardcoded too? comp["He3"] = 3.15247417638132e-04 * Y; comp["He4"] = Y - comp["He3"]; From c0c2b859466370509997fd5112b7eb3e9c10353b Mon Sep 17 00:00:00 2001 From: ctmbl Date: Mon, 3 Jul 2023 11:38:52 +0200 Subject: [PATCH 05/20] Check that metal-mix.cfg keys are metals in a specific set and print the remaining of the sum of fractions --- src/physics/composition.cpp | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/physics/composition.cpp b/src/physics/composition.cpp index 63ea61e0..3792383b 100644 --- a/src/physics/composition.cpp +++ b/src/physics/composition.cpp @@ -6,6 +6,7 @@ #include "physics.h" #include "parser.h" +bool _init_metals = true; double_map initial_composition(double X, double Z) { double_map comp; @@ -17,6 +18,7 @@ double_map initial_composition(double X, double Z) { // TODO: should this one be hardcoded too? comp["He3"] = 3.15247417638132e-04 * Y; comp["He4"] = Y - comp["He3"]; + // Note that the previous way of setting He3 He4 ensure that He3+He4 = Y // TODO: change this hardcoded char file[] = "metal-mix.cfg"; @@ -28,6 +30,7 @@ double_map initial_composition(double X, double Z) { for(std::string metal: metals){ comp[metal] = .0; } + double metals_fraction = .0; if(!fp.open(file)){ printf("Can't open configuration file %s\n", file); @@ -40,12 +43,27 @@ double_map initial_composition(double X, double Z) { printf("Syntax error in configuration file %s, line %d\n", file, line); exit(1); } + if(metals.find(arg) == metals.end()){ + printf("%s is unknown, possible metals composition are: ", val); + for(std::string metal: metals){ + printf("%s ", metal); + } + puts("\n"); + exit(1); + } + metals_fraction += atof(val); comp[arg] = Z * atof(val); } } fp.close(); comp["Ex"] = 1 - comp.sum(); + if(_init_metals){ + printf("A config file %s has been used to config metal composition\n", file); + printf("The sum of metals fractions (of Z sum of metal mass fraction) is %f, and should be as close to 1 as possible.\n", metals_fraction); + _init_metals = false; + } + return comp; } From 21ae1a7d7e324fb2c06ee36d17422a5f11a7d4f1 Mon Sep 17 00:00:00 2001 From: ctmbl Date: Mon, 3 Jul 2023 13:20:46 +0200 Subject: [PATCH 06/20] Update eos_freeeos to extract metals fractions from comp and choose most import isotops --- src/physics/composition.cpp | 4 +++- src/physics/eos_freeeos.cpp | 43 ++++++++++++++++++++----------------- 2 files changed, 26 insertions(+), 21 deletions(-) diff --git a/src/physics/composition.cpp b/src/physics/composition.cpp index 3792383b..10b20ef3 100644 --- a/src/physics/composition.cpp +++ b/src/physics/composition.cpp @@ -25,7 +25,9 @@ double_map initial_composition(double X, double Z) { char* arg = NULL; char* val = NULL; - std::set metals = {"C12","C13","N14","N15","O16","O17","Ne","Na","Mg","AL","Si","P","S","Cl","A","Ca","Ti","Cr","Mn","Fe","Ni"}; + std::set metals = {"C12","C13","N14","N15","O16","O17","Ne20","Ne22","Na23", + "Mg24","Mg25","Mg26","Al27","Si28","P31","S32","Cl35","Cl37", + "A40","Ca40","Ti","Cr","Mn55","Fe","Ni"}; // Initialization of comp for(std::string metal: metals){ comp[metal] = .0; diff --git a/src/physics/eos_freeeos.cpp b/src/physics/eos_freeeos.cpp index 87da2cbc..fcc7ec4d 100644 --- a/src/physics/eos_freeeos.cpp +++ b/src/physics/eos_freeeos.cpp @@ -87,26 +87,29 @@ int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p, 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 + + // Only isotops that exists in non negligible quantities have been considered + // according to their wikipedia page on the 3rd of July 2023 + eps[0] = comp["H"] / 1.008e0; // H + eps[1] = comp["He3"] / 3. + comp["He4"] / 4.; // He3 + He4 + eps[2] = comp["C12"] / 12. + comp["C13"] / 13.; // C12 + C13 + eps[3] = comp["N14"] / 14. + comp["N15"] / 15.; // N14 + N15 + eps[4] = comp["O16"] / 16. + comp["O17"] / 17.; // O16 + O17 + eps[5] = comp["Ne20"] / 20. + comp["Ne22"] / 22.; // Ne20 + Ne22 + eps[6] = comp["Na23"] / 23.; // Na23 + eps[7] = comp["Mg24"] / 24 + comp["Mg25"] / 25 + comp["Mg26"] / 26; // Mg24 + Mg25 + Mg26 + eps[8] = comp["Al27"] / 27.0; // Al27 + eps[9] = comp["Si28"] / 28.0; // Si28 + eps[10] = comp["P31"] / 31.0; // P31 + eps[11] = comp["S32"] / 32.0; // S32 + eps[12] = comp["Cl35"] / 35.0 + comp["Cl37"] / 37.0; // Cl35 + Cl37 + eps[13] = comp["A40"] / 40.0; // A40 + eps[14] = comp["Ca40"] / 40.0; // Ca40 + eps[15] = comp["Ti"] / 47.867; // Ti //there are too many isotops + eps[16] = comp["Cr"] / 52.0; // Cr //there are too many isotops + eps[17] = comp["Mn55"] / 55.0; // Mn + eps[18] = comp["Fe"] / 55.845; // Fe //there are too many isotops + eps[19] = comp["Ni"] / 58.693; // Ni //there are too many isotops double pi = p(i); double match_variable = log(pi); From fdaac5dcb1e1dbb2737e35e85254dcac3606e99a Mon Sep 17 00:00:00 2001 From: ctmbl Date: Thu, 13 Jul 2023 16:40:29 +0200 Subject: [PATCH 07/20] (REVERTED LATER) Add printf to debug --- src/star/star1d_class.cpp | 1 + src/star/star1d_solvers.cpp | 4 ++++ src/star/star2d_solvers.cpp | 3 +++ 3 files changed, 8 insertions(+) diff --git a/src/star/star1d_class.cpp b/src/star/star1d_class.cpp index 86ec087c..cb8253e0 100644 --- a/src/star/star1d_class.cpp +++ b/src/star/star1d_class.cpp @@ -156,6 +156,7 @@ int star1d::init(const char *input_file, const char *param_file, int argc, char } init_comp(); + printf("from star1d::init ->"); fill(); phi = solve_phi(); diff --git a/src/star/star1d_solvers.cpp b/src/star/star1d_solvers.cpp index c83b8ff1..dcbb8cec 100644 --- a/src/star/star1d_solvers.cpp +++ b/src/star/star1d_solvers.cpp @@ -9,6 +9,7 @@ #include "matplotlib.h" //--------------------------------------------------------------------- void star1d::fill() { + printf("star1d::fill "); Y0=1.-X0-Z0; init_comp(); @@ -93,6 +94,7 @@ double star1d::solve(solver *op) { } //--------------------------------------------------------------------- double star1d::solve(solver *op, matrix_map& error_map, int nit) { + printf("star1d::solve "); int info[5]; matrix rho0; double err,err2; @@ -177,6 +179,7 @@ double star1d::solve(solver *op, matrix_map& error_map, int nit) { rho0=rho; + printf("from star1d::solve ->"); fill(); err2=max(abs(rho-rho0));err=err2>err?err2:err; @@ -728,6 +731,7 @@ void star1d::check_jacobian(solver *op,const char *eqn) { B.map.remap(); } + printf("from star1d::check_jacobian ->"); B.fill(); i=op->get_id("rho"); diff --git a/src/star/star2d_solvers.cpp b/src/star/star2d_solvers.cpp index 832aeab5..16557656 100644 --- a/src/star/star2d_solvers.cpp +++ b/src/star/star2d_solvers.cpp @@ -11,6 +11,7 @@ /// nuclear reaction, mass definition, pi_c, Lambda, velocity, units, atmosphere, /// flatness, scaled keplerian angular velocity void star2d::fill() { + printf("star2d::fill "); Y0=1.-X0-Z0; init_comp(); @@ -37,6 +38,7 @@ void star2d::fill() { } void star2d::init_comp() { + printf("star2d::init_comp! "); // Update the object comp comp=initial_composition(X0,Z0)*ones(nr,nth); @@ -166,6 +168,7 @@ double star2d::solve(solver *op) { /// \brief Performs one step of the Newton algorithm to compute the star's /// internal structure. double star2d::solve(solver *op, matrix_map& error_map, int nit) { + printf("star2d::solve "); int info[5]; matrix rho0; double err,err2,h,dmax; From 5923d19dcd1b3fcd23ca4474554baf690e2824c3 Mon Sep 17 00:00:00 2001 From: ctmbl Date: Thu, 13 Jul 2023 16:41:39 +0200 Subject: [PATCH 08/20] (BREAKS ESTER INFO ON PURPOSE) Call init_comp only in star*d::init not in fill/read, remove useless redondant init_comp --- src/star/star1d_solvers.cpp | 1 - src/star/star2d_solvers.cpp | 3 +-- utils/starR/star1dR_class.cpp | 2 -- utils/starR/star2dR_class.cpp | 2 -- 4 files changed, 1 insertion(+), 7 deletions(-) diff --git a/src/star/star1d_solvers.cpp b/src/star/star1d_solvers.cpp index dcbb8cec..fb0a2507 100644 --- a/src/star/star1d_solvers.cpp +++ b/src/star/star1d_solvers.cpp @@ -11,7 +11,6 @@ void star1d::fill() { printf("star1d::fill "); Y0=1.-X0-Z0; - init_comp(); eq_state(); diff --git a/src/star/star2d_solvers.cpp b/src/star/star2d_solvers.cpp index 16557656..2a9fd44b 100644 --- a/src/star/star2d_solvers.cpp +++ b/src/star/star2d_solvers.cpp @@ -7,13 +7,12 @@ #include #include "symbolic.h" -/// \brief Initialize star's chemical composition, equation of state, opacity, +/// \brief Initialize star's equation of state, opacity, /// nuclear reaction, mass definition, pi_c, Lambda, velocity, units, atmosphere, /// flatness, scaled keplerian angular velocity void star2d::fill() { printf("star2d::fill "); Y0=1.-X0-Z0; - init_comp(); eq_state(); opacity(); diff --git a/utils/starR/star1dR_class.cpp b/utils/starR/star1dR_class.cpp index 91962f8c..9b1243ed 100644 --- a/utils/starR/star1dR_class.cpp +++ b/utils/starR/star1dR_class.cpp @@ -58,7 +58,6 @@ double star1dR::solve(solver *op) { void star1dR::fill() { - init_comp(); eq_state(); m=4*PI*(map.gl.I,rho*r*r)(0); M=m*rhoc*R*R*R; @@ -112,7 +111,6 @@ void star1dR::solve_definitions(solver *op) { double dXc=1e-8; Xc+=dXc; - init_comp(); nuclear(); opacity(); eq_state(); diff --git a/utils/starR/star2dR_class.cpp b/utils/starR/star2dR_class.cpp index d3c551ac..db5fa3b7 100644 --- a/utils/starR/star2dR_class.cpp +++ b/utils/starR/star2dR_class.cpp @@ -65,7 +65,6 @@ double star2dR::solve(solver *op) { void star2dR::fill() { - init_comp(); eq_state(); m=2*PI*(map.gl.I,rho*r*r*map.rz,map.leg.I_00)(0); M=m*rhoc*R*R*R; @@ -124,7 +123,6 @@ void star2dR::solve_definitions(solver *op) { double dXc=1e-8; Xc+=dXc; - init_comp(); nuclear(); opacity(); eq_state(); From 82c3e25deca21b154f181e0140213a0043779643 Mon Sep 17 00:00:00 2001 From: ctmbl Date: Thu, 13 Jul 2023 17:10:44 +0200 Subject: [PATCH 09/20] TO BE EDITED Crash ester if initial_composition is called multiple times --- src/physics/composition.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/physics/composition.cpp b/src/physics/composition.cpp index 10b20ef3..145d3a61 100644 --- a/src/physics/composition.cpp +++ b/src/physics/composition.cpp @@ -3,12 +3,20 @@ #endif #include #include +#include "utils.h" #include "physics.h" #include "parser.h" bool _init_metals = true; +bool first = true; + double_map initial_composition(double X, double Z) { + if(!first){ + //ester_err("initial_composition must be called only once"); + } + first = false; + double_map comp; file_parser fp; From f2181dbd5b9e1eb754ba37b59789a81af3ad32ab Mon Sep 17 00:00:00 2001 From: ctmbl Date: Mon, 17 Jul 2023 13:30:00 +0200 Subject: [PATCH 10/20] Remove fat old comment about Daniel's and Delphine's special star layers comp in stard::init_comp --- src/star/star2d_solvers.cpp | 34 ---------------------------------- 1 file changed, 34 deletions(-) diff --git a/src/star/star2d_solvers.cpp b/src/star/star2d_solvers.cpp index 2a9fd44b..670adb16 100644 --- a/src/star/star2d_solvers.cpp +++ b/src/star/star2d_solvers.cpp @@ -51,40 +51,6 @@ void star2d::init_comp() { } // and put some fraction Xc of X_envelope in the core comp.setblock(0,n-1,0,-1,initial_composition(Xc*X0,Z0)*ones(n,nth)); - -// Put a special composition in the last domain: Daniel's request -/* printf("WARNING : Daniel's superstar computed!\n"); - - n = 0; // Count the number of points to the before before last domain - for (int i=0; i Date: Mon, 17 Jul 2023 14:08:15 +0200 Subject: [PATCH 11/20] Add an important TODO about code moving to star*d_class.cpp --- src/star/star1d_solvers.cpp | 51 ++++++++++++++++++++++++++----------- src/star/star2d_solvers.cpp | 21 ++++++++++++--- 2 files changed, 53 insertions(+), 19 deletions(-) diff --git a/src/star/star1d_solvers.cpp b/src/star/star1d_solvers.cpp index fb0a2507..083f3799 100644 --- a/src/star/star1d_solvers.cpp +++ b/src/star/star1d_solvers.cpp @@ -7,7 +7,13 @@ #include "symbolic.h" #include "matplotlib.h" -//--------------------------------------------------------------------- + +// TODO: apply the following move (at the end of the work) +// The code between [BEGIN MOVE] and [END MODE] should be placed in star2d_class.cpp +// it hasn't anything to do with solving, but is required by reading, and ester info for exmaple +// whereas ester info doesn't require any other methods from star2d_solvers.cpp +// [BEGIN MOVE] + void star1d::fill() { printf("star1d::fill "); Y0=1.-X0-Z0; @@ -33,7 +39,9 @@ void star1d::fill() { Omega=0;Omega_bk=0;Ekman=0;Omegac=0; } -//--------------------------------------------------------------------- + +// [END MOVE] + solver *star1d::init_solver(int nvar_add) { int nvar; solver *op; @@ -49,7 +57,8 @@ solver *star1d::init_solver(int nvar_add) { return op; } -//--------------------------------------------------------------------- + + void star1d::register_variables(solver *op) { int i,var_nr[ndomains]; @@ -86,12 +95,14 @@ void star1d::register_variables(solver *op) { op->regvar_dep("nuc.eps"); } -//--------------------------------------------------------------------- + + double star1d::solve(solver *op) { matrix_map error_map; return solve(op, error_map, 0); } -//--------------------------------------------------------------------- + + double star1d::solve(solver *op, matrix_map& error_map, int nit) { printf("star1d::solve "); int info[5]; @@ -186,7 +197,8 @@ double star1d::solve(solver *op, matrix_map& error_map, int nit) { return err; } -//--------------------------------------------------------------------- + + void star1d::update_map(matrix dR) { if(ndomains==1) return; @@ -203,7 +215,8 @@ void star1d::update_map(matrix dR) { map.R=R0+h*dR; } } -//--------------------------------------------------------------------- + + void star1d::solve_definitions(solver *op) { // substitution of the dep variables by their expression // with the basic variables @@ -257,7 +270,8 @@ void star1d::solve_definitions(solver *op) { op->add_d("nuc.eps","log_Tc",nuc.dlneps_lnT*nuc.eps); } -//--------------------------------------------------------------------- + + void star1d::solve_poisson(solver *op) { int n,j0; matrix rhs; @@ -305,7 +319,8 @@ void star1d::solve_poisson(solver *op) { } op->set_rhs("Phi",rhs); } -//--------------------------------------------------------------------- + + void star1d::solve_pressure(solver *op) { int n,j0; matrix rhs_p,rhs_pi_c; @@ -341,7 +356,8 @@ void star1d::solve_pressure(solver *op) { op->set_rhs(eqn,rhs_p); op->set_rhs("pi_c",rhs_pi_c); } -//--------------------------------------------------------------------- + + void star1d::solve_temp(solver *op) { int n,j0,j1; matrix q; @@ -550,7 +566,8 @@ void star1d::solve_temp(solver *op) { op->set_rhs("Lambda",rhs_Lambda); } -//--------------------------------------------------------------------- + + void star1d::solve_dim(solver *op) { // write the equations defining the global variables: // Radius, R, Tc, Pc and non-dim mass. @@ -615,7 +632,8 @@ void star1d::solve_dim(solver *op) { op->set_rhs("log_R",rhs); } -//--------------------------------------------------------------------- + + void star1d::solve_map(solver *op) { int n,j0,j1; matrix rhs=zeros(ndomains,1); @@ -663,7 +681,8 @@ void star1d::solve_map(solver *op) { } op->set_rhs("Ri",rhs); } -//--------------------------------------------------------------------- + + void star1d::solve_gsup(solver *op) { // Comfort equation // Code the equation gsup=4*pi*G*rho_c*R*m @@ -679,7 +698,8 @@ void star1d::solve_gsup(solver *op) { op->set_rhs("gsup",zeros(1,1)); } -//--------------------------------------------------------------------- + + void star1d::solve_Teff(solver *op) { matrix q,Te,F; int n=ndomains-1; @@ -705,7 +725,8 @@ void star1d::solve_Teff(solver *op) { */ op->set_rhs("Teff",zeros(1,1)); } -//--------------------------------------------------------------------- + + void star1d::check_jacobian(solver *op,const char *eqn) { star1d B; matrix rhs,drhs,drhs2,qq; diff --git a/src/star/star2d_solvers.cpp b/src/star/star2d_solvers.cpp index 670adb16..f8bd6d99 100644 --- a/src/star/star2d_solvers.cpp +++ b/src/star/star2d_solvers.cpp @@ -7,6 +7,12 @@ #include #include "symbolic.h" +// TODO: apply the following move (at the end of the work) +// The code between [BEGIN MOVE] and [END MODE] should be placed in star2d_class.cpp +// it hasn't anything to do with solving, but is required by reading, and ester info for exmaple +// whereas ester info doesn't require any other methods from star2d_solvers.cpp +// [BEGIN MOVE] + /// \brief Initialize star's equation of state, opacity, /// nuclear reaction, mass definition, pi_c, Lambda, velocity, units, atmosphere, /// flatness, scaled keplerian angular velocity @@ -53,6 +59,8 @@ void star2d::init_comp() { comp.setblock(0,n-1,0,-1,initial_composition(Xc*X0,Z0)*ones(n,nth)); } +// [END MOVE] + void star2d::calc_veloc() { // vr=rz*V^zeta+rt*V^theta, vt=r*V^theta //vr=(G,map.leg.D_11)/r+(map.rt/r+cos(th)/sin(th))/r*G; @@ -65,6 +73,7 @@ void star2d::calc_veloc() { vt/=rho; } + solver *star2d::init_solver(int nvar_add) { int nvar; solver *op; @@ -80,6 +89,7 @@ solver *star2d::init_solver(int nvar_add) { return op; } + void star2d::register_variables(solver *op) { int i,var_nr[ndomains+1]; @@ -125,11 +135,13 @@ void star2d::register_variables(solver *op) { } + double star2d::solve(solver *op) { matrix_map error_map; return solve(op, error_map, 0); } + /// \brief Performs one step of the Newton algorithm to compute the star's /// internal structure. double star2d::solve(solver *op, matrix_map& error_map, int nit) { @@ -249,9 +261,11 @@ double star2d::solve(solver *op, matrix_map& error_map, int nit) { } + // Special treatment for updating R_i because we need // a different relaxation parameter "h". + void star2d::update_map(matrix dR) { double h=1,dmax=config.newton_dmax; @@ -267,6 +281,7 @@ void star2d::update_map(matrix dR) { } + /// \brief insert the definitions depending on opacity and eos tables into the solver, /// and the definitions used by the mapping (eta,deta,Ri,dRi,...), and the entropy void star2d::solve_definitions(solver *op) { @@ -887,7 +902,6 @@ void star2d::solve_temp(solver *op) { /// \brief Writes the equations for the dimensional quantities (T_c, rho_c, R, etc.) - void star2d::solve_dim(solver *op) { int n,j0,j1; matrix q,rhs; @@ -1077,8 +1091,6 @@ void star2d::solve_map(solver *op) { } - - /// \brief Equation setting the equatorial angular velocity void star2d::solve_Omega(solver *op) { int n; @@ -1103,7 +1115,6 @@ void star2d::solve_Omega(solver *op) { } - /// \brief Equation giving the effective surface gravity gsup /// gsup=(-\vn\cdot\grad P)/rho void star2d::solve_gsup(solver *op) { @@ -1140,6 +1151,7 @@ void star2d::solve_gsup(solver *op) { } + /// \brief Equation setting the surface effective temperature /// Derived from sigma T_e^4 = -xi\vn\cdot\gradT void star2d::solve_Teff(solver *op) { @@ -1177,6 +1189,7 @@ void star2d::solve_Teff(solver *op) { } + /// \brief Equation setting the 'simple' atmosphere model equations /// 'simple' == the polytropic model for the atmosphere /// checked but approximate on ... "p",-1/p.row(-1)/(n_atm+1) ... From 85d8c3d75756c7fac059edd5c054f84023e51c9f Mon Sep 17 00:00:00 2001 From: ctmbl Date: Mon, 17 Jul 2023 14:33:09 +0200 Subject: [PATCH 12/20] Introduce star2d::init_metal_mix method to replace old initial_composition --- src/include/star.h | 2 ++ src/star/star2d_solvers.cpp | 58 +++++++++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+) diff --git a/src/include/star.h b/src/include/star.h index 04e34b00..8fcbf975 100644 --- a/src/include/star.h +++ b/src/include/star.h @@ -37,6 +37,7 @@ class star2d { matrix rho, phi, p, T; matrix phiex; matrix vr, vt, G, w; + double_map m_metal_mix; // the m_ prefix is used to explicit that it's a class' member https://stackoverflow.com/a/13018190/22158934 composition_map comp; double X0, Y0, Z0; double R, M; @@ -97,6 +98,7 @@ class star2d { virtual void dump_info(); + virtual void init_metal_mix(); virtual void init_comp(); virtual solver *init_solver(int nvar_add=0); diff --git a/src/star/star2d_solvers.cpp b/src/star/star2d_solvers.cpp index f8bd6d99..e51252c0 100644 --- a/src/star/star2d_solvers.cpp +++ b/src/star/star2d_solvers.cpp @@ -4,8 +4,10 @@ #include "star.h" #include #include +#include #include #include "symbolic.h" +#include "utils.h" // TODO: apply the following move (at the end of the work) // The code between [BEGIN MOVE] and [END MODE] should be placed in star2d_class.cpp @@ -42,6 +44,62 @@ void star2d::fill() { } + +void star2d::init_metal_mix() { + // a security to avoid initializing multiple time the metal mix + // it should be done only once, then m_metal_mix should be used + if(!m_metal_mix.empty()){ + ester_err("init_metal_mix must be called only once"); + } + + file_parser fp; + + // TODO: change this hardcoded path + char file[] = "metal-mix.cfg"; + + char* arg = NULL; + char* val = NULL; + std::set metals = {"C12","C13","N14","N15","O16","O17","Ne20","Ne22","Na23", + "Mg24","Mg25","Mg26","Al27","Si28","P31","S32","Cl35","Cl37", + "A40","Ca40","Ti","Cr","Mn55","Fe","Ni"}; + // Initialization of m_metal_mix + for(std::string metal: metals){ + m_metal_mix[metal] = .0; + } + double metals_fraction = .0; + + if(!fp.open(file)){ + printf("Can't open configuration file %s\n", file); + perror("Error:"); + exit(1); + } else { + int line; + while(line = fp.get(arg,val)) { + if(val == NULL){ + printf("Syntax error in configuration file %s, line %d\n", file, line); + exit(1); + } + if(metals.find(arg) == metals.end()){ + // the metal specified in the config file isn't supported + printf("%s is unknown, possible metals composition are: ", val); + for(std::string metal: metals){ + printf("%s ", metal); + } + puts("\n"); + exit(1); + } + metals_fraction += atof(val); + m_metal_mix[arg] = Z0 * atof(val); + } + } + fp.close(); + + // will be removed: + printf("A config file %s has been used to config metal composition\n", file); + printf("The sum of metals fractions (of Z sum of metal mass fraction) is %f, and should be as close to 1 as possible.\n", metals_fraction); +} + + void star2d::init_comp() { printf("star2d::init_comp! "); // Update the object comp From 089cad0ef83d15978b0017fdd6d558b81668bf2c Mon Sep 17 00:00:00 2001 From: ctmbl Date: Mon, 17 Jul 2023 15:22:47 +0200 Subject: [PATCH 13/20] Update star2d::init_comp to use new star2d::init_metal_mix and remove old initial_composition --- python/ester_wrap.cxx | 32 ---------------- python/ester_wrap.py | 3 -- src/include/physics.h | 2 - src/physics/composition.cpp | 75 +------------------------------------ src/physics/eos_freeeos.cpp | 53 +++++++++++++------------- src/star/star1d_solvers.cpp | 1 - src/star/star2d_solvers.cpp | 35 ++++++++++++----- 7 files changed, 53 insertions(+), 148 deletions(-) diff --git a/python/ester_wrap.cxx b/python/ester_wrap.cxx index 1cdb2411..c0b9a35a 100644 --- a/python/ester_wrap.cxx +++ b/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}, diff --git a/python/ester_wrap.py b/python/ester_wrap.py index 73436d59..862715a0 100644 --- a/python/ester_wrap.py +++ b/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): diff --git a/src/include/physics.h b/src/include/physics.h index 223305cd..8e3cd487 100644 --- a/src/include/physics.h +++ b/src/include/physics.h @@ -70,7 +70,5 @@ int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p, 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); - #endif diff --git a/src/physics/composition.cpp b/src/physics/composition.cpp index 145d3a61..4a6fa36c 100644 --- a/src/physics/composition.cpp +++ b/src/physics/composition.cpp @@ -1,81 +1,8 @@ #ifndef WITH_CMAKE #include "ester-config.h" #endif -#include -#include -#include "utils.h" -#include "physics.h" -#include "parser.h" - -bool _init_metals = true; - -bool first = true; - -double_map initial_composition(double X, double Z) { - if(!first){ - //ester_err("initial_composition must be called only once"); - } - first = false; - - double_map comp; - file_parser fp; - - double Y = 1. - (X + Z); - - comp["H"] = X; - // TODO: should this one be hardcoded too? - comp["He3"] = 3.15247417638132e-04 * Y; - comp["He4"] = Y - comp["He3"]; - // Note that the previous way of setting He3 He4 ensure that He3+He4 = Y - // TODO: change this hardcoded - char file[] = "metal-mix.cfg"; - - char* arg = NULL; - char* val = NULL; - std::set metals = {"C12","C13","N14","N15","O16","O17","Ne20","Ne22","Na23", - "Mg24","Mg25","Mg26","Al27","Si28","P31","S32","Cl35","Cl37", - "A40","Ca40","Ti","Cr","Mn55","Fe","Ni"}; - // Initialization of comp - for(std::string metal: metals){ - comp[metal] = .0; - } - double metals_fraction = .0; - - if(!fp.open(file)){ - printf("Can't open configuration file %s\n", file); - perror("Error:"); - exit(1); - } else { - int line; - while(line = fp.get(arg,val)) { - if(val == NULL){ - printf("Syntax error in configuration file %s, line %d\n", file, line); - exit(1); - } - if(metals.find(arg) == metals.end()){ - printf("%s is unknown, possible metals composition are: ", val); - for(std::string metal: metals){ - printf("%s ", metal); - } - puts("\n"); - exit(1); - } - metals_fraction += atof(val); - comp[arg] = Z * atof(val); - } - } - fp.close(); - comp["Ex"] = 1 - comp.sum(); - - if(_init_metals){ - printf("A config file %s has been used to config metal composition\n", file); - printf("The sum of metals fractions (of Z sum of metal mass fraction) is %f, and should be as close to 1 as possible.\n", metals_fraction); - _init_metals = false; - } - - return comp; -} +#include "physics.h" matrix composition_map::X() const { diff --git a/src/physics/eos_freeeos.cpp b/src/physics/eos_freeeos.cpp index fcc7ec4d..0b0de635 100644 --- a/src/physics/eos_freeeos.cpp +++ b/src/physics/eos_freeeos.cpp @@ -21,6 +21,32 @@ extern "C" { double *entropy, int *iteration_count); } +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 matrix &X, double Z, const matrix &T, const matrix &p, matrix &rho, eos_struct &eos) { @@ -84,32 +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 Date: Tue, 18 Jul 2023 17:18:24 +0200 Subject: [PATCH 14/20] Make the He3/He4 ratio a star2d attr and hardcode it in star2d::init_comp --- src/include/star.h | 4 +++- src/star/star2d_solvers.cpp | 3 ++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/include/star.h b/src/include/star.h index 8fcbf975..a0f7c313 100644 --- a/src/include/star.h +++ b/src/include/star.h @@ -37,7 +37,9 @@ class star2d { matrix rho, phi, p, T; matrix phiex; matrix vr, vt, G, w; - double_map m_metal_mix; // the m_ prefix is used to explicit that it's a class' member https://stackoverflow.com/a/13018190/22158934 + // 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; diff --git a/src/star/star2d_solvers.cpp b/src/star/star2d_solvers.cpp index 7a296ed8..291d01d2 100644 --- a/src/star/star2d_solvers.cpp +++ b/src/star/star2d_solvers.cpp @@ -102,6 +102,7 @@ void star2d::init_metal_mix() { void star2d::init_comp() { printf("star2d::init_comp! "); Y0 = 1. - (X0 + Z0); + m_He_isotopic_ratio = 3.15247417638132e-04; init_metal_mix(); @@ -113,7 +114,7 @@ void star2d::init_comp() { temp_chemical_mix["H"] = X0; // TODO: should the following be hardcoded? - temp_chemical_mix["He3"] = 3.15247417638132e-04 * Y0; + temp_chemical_mix["He3"] = m_He_isotopic_ratio * Y0; temp_chemical_mix["He4"] = Y0 - temp_chemical_mix["He3"]; // Note that the previous way of setting He3 He4 enforce He3+He4 = Y From d7f56a367e7e71b2cd004827c01895c56cc16596 Mon Sep 17 00:00:00 2001 From: ctmbl Date: Tue, 18 Jul 2023 17:20:43 +0200 Subject: [PATCH 15/20] Add 2 matrix_map::operator() methods that uses matrix::operator() --- src/include/matrix.h | 2 ++ src/matrix/matrix_map.cpp | 22 ++++++++++++++++++++++ 2 files changed, 24 insertions(+) diff --git a/src/include/matrix.h b/src/include/matrix.h index f2e1303f..5f61b753 100644 --- a/src/include/matrix.h +++ b/src/include/matrix.h @@ -289,7 +289,9 @@ class matrix_map : public std::map { inline matrix &operator[](const std::string& key) {return std::map::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; diff --git a/src/matrix/matrix_map.cpp b/src/matrix/matrix_map.cpp index 58336024..d756ac70 100644 --- a/src/matrix/matrix_map.cpp +++ b/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; From 4bb473a07ee084038f1ee811fc5f5bc28934c9f3 Mon Sep 17 00:00:00 2001 From: ctmbl Date: Tue, 18 Jul 2023 17:29:08 +0200 Subject: [PATCH 16/20] Fix wrongly named nuc_calc signature --- src/include/physics.h | 2 +- src/physics/physics.cpp | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/include/physics.h b/src/include/physics.h index 8e3cd487..f5cb4e17 100644 --- a/src/include/physics.h +++ b/src/include/physics.h @@ -36,7 +36,7 @@ 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 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); diff --git a/src/physics/physics.cpp b/src/physics/physics.cpp index 457a4246..a0409d5d 100644 --- a/src/physics/physics.cpp +++ b/src/physics/physics.cpp @@ -48,15 +48,15 @@ int eos_calc(const matrix &X,double Z,const matrix &T,const matrix &p, } -int nuc_calc(const matrix_map &X,const matrix &T,const matrix &rho, +int nuc_calc(const composition_map &comp, const matrix &T, const matrix &rho, nuc_struct &nuc) { int error=0; if(!strcmp(nuc.name,"simple")) { - error=nuc_simple(X,T,rho,nuc); + error = nuc_simple(comp, T, rho, nuc); } else if(!strcmp(nuc.name,"cesam")) { - error=nuc_cesam(X,T,rho,nuc); + error = nuc_cesam(comp, T, rho, nuc); } else { ester_err("Unknown nuc. reac. type: %s",nuc.name); return 1; From a07ed233d87a8e3394feda666009938e7b522077 Mon Sep 17 00:00:00 2001 From: ctmbl Date: Thu, 20 Jul 2023 11:54:33 +0200 Subject: [PATCH 17/20] Clean src/physics/eos_opal.cpp, removing useless code --- src/physics/eos_opal.cpp | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/src/physics/eos_opal.cpp b/src/physics/eos_opal.cpp index c36b3978..1f1a873c 100644 --- a/src/physics/eos_opal.cpp +++ b/src/physics/eos_opal.cpp @@ -14,16 +14,13 @@ extern"C" { extern struct{ double esact,eos[10]; } eeos_; - extern struct{ - int itime; - } lreadco_; } int eos_opal(const matrix &X,double Z,const matrix &T,const matrix &p, matrix &rho,eos_struct &eos) { matrix t6,p_mb; - int i,N,error=0; + int N; static double Z_table=-99; @@ -31,16 +28,10 @@ int eos_opal(const matrix &X,double Z,const matrix &T,const matrix &p, p_mb=p*1e-12; if(Z!=Z_table) { - lreadco_.itime=0; zfs_interp_eos5_(&Z); Z_table=Z; } - if(error) { - printf("Can't initialize OPAL EOS table\n"); - return error; - } - rho.dim(T.nrows(),T.ncols()); eos.s.dim(T.nrows(),T.ncols()); eos.G1.dim(T.nrows(),T.ncols()); @@ -57,7 +48,7 @@ int eos_opal(const matrix &X,double Z,const matrix &T,const matrix &p, N=T.nrows()*T.ncols(); double Xi,Zi,t6i,p_mbi,rhoi; - for(i=0;i Date: Thu, 20 Jul 2023 11:55:05 +0200 Subject: [PATCH 18/20] Change eos_* signature to pass the whole star composition --- src/include/physics.h | 26 +++++++++++++------------- src/physics/atm_onelayer.cpp | 8 ++++---- src/physics/eos_freeeos.cpp | 5 +++-- src/physics/eos_ideal.cpp | 9 +++++---- src/physics/eos_idealrad.cpp | 7 ++++--- src/physics/eos_opal.cpp | 20 +++++++++----------- src/physics/physics.cpp | 18 +++++++++--------- src/star/star_phys.cpp | 10 ++++++---- 8 files changed, 53 insertions(+), 50 deletions(-) diff --git a/src/include/physics.h b/src/include/physics.h index f5cb4e17..190f5caa 100644 --- a/src/include/physics.h +++ b/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 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,17 +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); +int atm_onelayer(const composition_map &, const matrix &g, const matrix &Teff, + const char *eos_name, const char *opa_name, atm_struct &atm); #endif diff --git a/src/physics/atm_onelayer.cpp b/src/physics/atm_onelayer.cpp index db366430..9f6dbd18 100644 --- a/src/physics/atm_onelayer.cpp +++ b/src/physics/atm_onelayer.cpp @@ -7,8 +7,8 @@ #include #include -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; diff --git a/src/physics/eos_freeeos.cpp b/src/physics/eos_freeeos.cpp index 0b0de635..253c693a 100644 --- a/src/physics/eos_freeeos.cpp +++ b/src/physics/eos_freeeos.cpp @@ -47,7 +47,7 @@ void set_eps(double* eps, const double_map &local_chemical_mix) { eps[19] = local_chemical_mix["Ni"] / 58.693; // Ni //there are too many isotops } -int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p, +int eos_freeeos(const composition_map &chemical_comp, const matrix &T, const matrix &p, matrix &rho, eos_struct &eos) { double t; @@ -127,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; diff --git a/src/physics/eos_ideal.cpp b/src/physics/eos_ideal.cpp index 48169903..d2976d74 100644 --- a/src/physics/eos_ideal.cpp +++ b/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()); diff --git a/src/physics/eos_idealrad.cpp b/src/physics/eos_idealrad.cpp index 439021c5..be163bcb 100644 --- a/src/physics/eos_idealrad.cpp +++ b/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; diff --git a/src/physics/eos_opal.cpp b/src/physics/eos_opal.cpp index 1f1a873c..0a7d6e6f 100644 --- a/src/physics/eos_opal.cpp +++ b/src/physics/eos_opal.cpp @@ -16,22 +16,15 @@ extern"C" { } eeos_; } -int eos_opal(const matrix &X,double Z,const matrix &T,const matrix &p, - matrix &rho,eos_struct &eos) { +int eos_opal(const composition_map &chemical_comp, const matrix &T, const matrix &p, + matrix &rho, eos_struct &eos) { matrix t6,p_mb; int N; - static double Z_table=-99; - t6=T*1e-6; p_mb=p*1e-12; - if(Z!=Z_table) { - zfs_interp_eos5_(&Z); - Z_table=Z; - } - rho.dim(T.nrows(),T.ncols()); eos.s.dim(T.nrows(),T.ncols()); eos.G1.dim(T.nrows(),T.ncols()); @@ -49,8 +42,13 @@ int eos_opal(const matrix &X,double Z,const matrix &T,const matrix &p, double Xi,Zi,t6i,p_mbi,rhoi; for(int i = 0; i < N; i++) { - Xi=X(i);Zi=Z;t6i=t6(i);p_mbi=p_mb(i); - eos5_xtrin_(&Xi,&Zi,&t6i,&p_mbi,&rhoi); + Xi = chemical_comp.X()(i); + Zi = chemical_comp.Z()(i); + t6i = t6(i); + p_mbi = p_mb(i); + + zfs_interp_eos5_(&Zi); + eos5_xtrin_(&Xi,&Zi,&t6i,&p_mbi,&rhoi); rho(i)=rhoi; eos.s(i)=1e6*(*(eeos_.eos+2)); eos.G1(i)=*(eeos_.eos+7); diff --git a/src/physics/physics.cpp b/src/physics/physics.cpp index a0409d5d..3f4bc5fa 100644 --- a/src/physics/physics.cpp +++ b/src/physics/physics.cpp @@ -26,19 +26,19 @@ int opa_calc(const matrix &X,double Z,const matrix &T,const matrix &rho, return error; } -int eos_calc(const matrix &X,double Z,const matrix &T,const matrix &p, - matrix &rho,eos_struct &eos) { +int eos_calc(const composition_map &chemical_comp, const matrix &T, const matrix &p, + matrix &rho, eos_struct &eos) { int error=0; if(!strcmp(eos.name,"ideal")) - error=eos_ideal(X,Z,T,p,rho,eos); + error = eos_ideal(chemical_comp, T, p, rho, eos); else if(!strcmp(eos.name,"ideal+rad")) - error=eos_idealrad(X,Z,T,p,rho,eos); + error = eos_idealrad(chemical_comp, T, p, rho, eos); else if(!strcmp(eos.name,"opal")) - error=eos_opal(X,Z,T,p,rho,eos); + error = eos_opal(chemical_comp, T, p, rho, eos); else if(!strcmp(eos.name,"freeeos")) - error = eos_freeeos(X, Z, T, p, rho, eos); + error = eos_freeeos(chemical_comp, T, p, rho, eos); else { ester_err("Unknown equation of state: %s",eos.name); return 1; @@ -66,13 +66,13 @@ int nuc_calc(const composition_map &comp, const matrix &T, const matrix &rho, } -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 &chemical_comp, const matrix &g, const matrix &Teff, + const char *eos_name, const char *opa_name, atm_struct &atm) { int error=0; if(!strcmp(atm.name,"onelayer")) { - error=atm_onelayer(X,Z,g,Teff,eos_name,opa_name,atm); + error = atm_onelayer(chemical_comp, g, Teff, eos_name, opa_name, atm); } else { ester_err("Unknown atmosphere type: %s", atm.name); return 1; diff --git a/src/star/star_phys.cpp b/src/star/star_phys.cpp index 0a9153fc..7d1435c8 100644 --- a/src/star/star_phys.cpp +++ b/src/star/star_phys.cpp @@ -35,10 +35,12 @@ void star2d::eq_state() { int error; matrix rhoc_m(1,1); - eos_calc(comp.X()(0,0)*ones(1,1),Z0,ones(1,1)*Tc,ones(1,1)*pc,rhoc_m,eos); - rhoc=rhoc_m(0); + // TODO: is the following really usefull??? + // couldn't we just get rhoc from rho(0) ??? instead of this weird call only computing central value + eos_calc(comp(0,0)*ones(1,1), Tc*ones(1,1), pc*ones(1,1), rhoc_m, eos); + rhoc = rhoc_m(0); - error=eos_calc(comp.X(),Z0,T*Tc,p*pc,rho,eos); + error = eos_calc(comp, T*Tc, p*pc, rho, eos); // rhoc=rho(0); @@ -61,7 +63,7 @@ void star2d::atmosphere() { } - error=atm_calc(comp.X(),Z0,gsup(),Teff(),eos.name,opa.name,atm); + error = atm_calc(comp, gsup(), Teff(), eos.name, opa.name, atm); if(error) exit(1); From 09553f6e9c83c7a2c04bf7478a472a7d396f4692 Mon Sep 17 00:00:00 2001 From: ctmbl Date: Wed, 19 Jul 2023 15:46:45 +0200 Subject: [PATCH 19/20] Fix OPAL interpolation and double comparison pb --- src/physics/EOS5_xtrin.F | 2 +- src/physics/eos_opal.cpp | 27 ++++++++++++++++++++------- 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/src/physics/EOS5_xtrin.F b/src/physics/EOS5_xtrin.F index fe6d1043..d0af5778 100644 --- a/src/physics/EOS5_xtrin.F +++ b/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 diff --git a/src/physics/eos_opal.cpp b/src/physics/eos_opal.cpp index 0a7d6e6f..1da56aaa 100644 --- a/src/physics/eos_opal.cpp +++ b/src/physics/eos_opal.cpp @@ -16,12 +16,25 @@ extern"C" { } eeos_; } +static double Z_OPAL = -1.; + int eos_opal(const composition_map &chemical_comp, const matrix &T, const matrix &p, matrix &rho, eos_struct &eos) { - + matrix t6,p_mb; int N; - + + if(Z_OPAL == -1.){ + Z_OPAL = chemical_comp.Z()(0,0); + ester_warn( + "OPAL EOS doesn't support non-Z-Uniform stars (interpolation table computing is too long).\n" + "For **EOS ONLY** Z will be uniform.\n" + "Taking Z = Z(0,0) = %.6f", + Z_OPAL + ); + zfs_interp_eos5_(&Z_OPAL); + } + t6=T*1e-6; p_mb=p*1e-12; @@ -40,15 +53,15 @@ int eos_opal(const composition_map &chemical_comp, const matrix &T, const matrix N=T.nrows()*T.ncols(); - double Xi,Zi,t6i,p_mbi,rhoi; + double Xi, t6i, p_mbi, rhoi; for(int i = 0; i < N; i++) { Xi = chemical_comp.X()(i); - Zi = chemical_comp.Z()(i); + t6i = t6(i); p_mbi = p_mb(i); - zfs_interp_eos5_(&Zi); - eos5_xtrin_(&Xi,&Zi,&t6i,&p_mbi,&rhoi); + // NOTE: Z_OPAL is used, not Z(i) see ester_warn above + eos5_xtrin_(&Xi,&Z_OPAL,&t6i,&p_mbi,&rhoi); rho(i)=rhoi; eos.s(i)=1e6*(*(eeos_.eos+2)); eos.G1(i)=*(eeos_.eos+7); @@ -65,7 +78,7 @@ int eos_opal(const composition_map &chemical_comp, const matrix &T, const matrix " X = %e\n" " Z = %e\n" " T = %e\n" - " p = %e", Xi, Zi, t6i, p_mbi); + " p = %e", Xi, Z_OPAL, t6i, p_mbi); } } if(exist(rho==-9e99)) { From 0a8c5c9bfd89668f545d9f1df8b2d242d67f5fa5 Mon Sep 17 00:00:00 2001 From: ctmbl Date: Wed, 19 Jul 2023 15:57:02 +0200 Subject: [PATCH 20/20] Revert "Add printf to debug" This reverts commit 208454be780d1d82365c285ef5c5fe6b8c91f0af. --- src/star/star1d_class.cpp | 1 - src/star/star1d_solvers.cpp | 5 ----- src/star/star2d_solvers.cpp | 4 ---- 3 files changed, 10 deletions(-) diff --git a/src/star/star1d_class.cpp b/src/star/star1d_class.cpp index cb8253e0..86ec087c 100644 --- a/src/star/star1d_class.cpp +++ b/src/star/star1d_class.cpp @@ -156,7 +156,6 @@ int star1d::init(const char *input_file, const char *param_file, int argc, char } init_comp(); - printf("from star1d::init ->"); fill(); phi = solve_phi(); diff --git a/src/star/star1d_solvers.cpp b/src/star/star1d_solvers.cpp index 46acbe22..8ca02754 100644 --- a/src/star/star1d_solvers.cpp +++ b/src/star/star1d_solvers.cpp @@ -15,8 +15,6 @@ // [BEGIN MOVE] void star1d::fill() { - printf("star1d::fill "); - eq_state(); opacity(); @@ -103,7 +101,6 @@ double star1d::solve(solver *op) { double star1d::solve(solver *op, matrix_map& error_map, int nit) { - printf("star1d::solve "); int info[5]; matrix rho0; double err,err2; @@ -188,7 +185,6 @@ double star1d::solve(solver *op, matrix_map& error_map, int nit) { rho0=rho; - printf("from star1d::solve ->"); fill(); err2=max(abs(rho-rho0));err=err2>err?err2:err; @@ -750,7 +746,6 @@ void star1d::check_jacobian(solver *op,const char *eqn) { B.map.remap(); } - printf("from star1d::check_jacobian ->"); B.fill(); i=op->get_id("rho"); diff --git a/src/star/star2d_solvers.cpp b/src/star/star2d_solvers.cpp index 291d01d2..e4efde0b 100644 --- a/src/star/star2d_solvers.cpp +++ b/src/star/star2d_solvers.cpp @@ -19,8 +19,6 @@ /// nuclear reaction, mass definition, pi_c, Lambda, velocity, units, atmosphere, /// flatness, scaled keplerian angular velocity void star2d::fill() { - printf("star2d::fill "); - eq_state(); opacity(); nuclear(); @@ -100,7 +98,6 @@ void star2d::init_metal_mix() { void star2d::init_comp() { - printf("star2d::init_comp! "); Y0 = 1. - (X0 + Z0); m_He_isotopic_ratio = 3.15247417638132e-04; @@ -219,7 +216,6 @@ double star2d::solve(solver *op) { /// \brief Performs one step of the Newton algorithm to compute the star's /// internal structure. double star2d::solve(solver *op, matrix_map& error_map, int nit) { - printf("star2d::solve "); int info[5]; matrix rho0; double err,err2,h,dmax;