From ada8a1fb8117c7caa65be23ff19112d92fd83f33 Mon Sep 17 00:00:00 2001 From: Cesar Diaz Date: Thu, 21 Oct 2021 12:43:32 -0500 Subject: [PATCH 1/4] Read harm files with multiple electron models through keln parameter --- model/iharm/example.par | 3 +++ model/iharm/model.c | 12 +++++++++--- model/iharm/model_params.h | 4 ++-- src/par.c | 4 ++++ src/par.h | 2 ++ 5 files changed, 20 insertions(+), 5 deletions(-) diff --git a/model/iharm/example.par b/model/iharm/example.par index e223c51..235aa19 100644 --- a/model/iharm/example.par +++ b/model/iharm/example.par @@ -1,3 +1,6 @@ +# 9 Kawazura, 10 Werner, 11 Rowan, 12 Sharma, default is 9 +keln 10 + # example.par -- ipole example parameter file for iharm problem # Common diff --git a/model/iharm/model.c b/model/iharm/model.c index 33ada1c..878c9f8 100644 --- a/model/iharm/model.c +++ b/model/iharm/model.c @@ -19,7 +19,7 @@ #define NVAR (10) #define USE_FIXED_TPTE (0) -#define USE_MIXED_TPTE (1) +#define USE_MIXED_TPTE (0) // KEL models are actually used #define NSUP (3) #define USE_GEODESIC_SIGMACUT (1) @@ -37,6 +37,7 @@ double RHO_unit; double U_unit; double B_unit; double Te_unit; +int keln; // MOLECULAR WEIGHTS static double Ne_factor = 1.; // e.g., used for He with 2 protons+neutrons per 2 electrons @@ -149,6 +150,9 @@ void try_set_model_parameter(const char *word, const char *value) set_by_word_val(word, value, "dump_max", &dumpmax, TYPE_INT); set_by_word_val(word, value, "dump_skip", &dumpskip, TYPE_INT); dumpidx = dumpmin; + + // for selecting electron model + set_by_word_val(word, value, "keln", &keln, TYPE_INT); } // Advance through dumps until we are closer to the next set @@ -524,6 +528,7 @@ void set_units() fprintf(stderr,"MBH: %g [Msun]\n",MBH/MSUN); fprintf(stderr,"L,T,M units: %g [cm] %g [s] %g [g]\n",L_unit,T_unit,M_unit) ; fprintf(stderr,"rho,u,B units: %g [g cm^-3] %g [g cm^-1 s^-2] %g [G] \n",RHO_unit,U_unit,B_unit) ; + fprintf(stderr,"Model is: %d \n", keln) ; } void init_physical_quantities(int n) @@ -1205,6 +1210,7 @@ void output_hdf5() hdf5_write_single_val(&M_unit, "M_unit", H5T_IEEE_F64LE); hdf5_write_single_val(&T_unit, "T_unit", H5T_IEEE_F64LE); hdf5_write_single_val(&Te_unit, "Thetae_unit", H5T_IEEE_F64LE); + hdf5_write_single_val(&keln, "keln", H5T_IEEE_F64LE); hdf5_set_directory("/"); @@ -1797,9 +1803,9 @@ void load_iharm_data(int n, char *fnam, int dumpidx, int verbose) hdf5_read_array(data[n]->p[B3][0][0], "prims", 4, fdims, fstart, fcount, mdims, mstart, H5T_IEEE_F64LE); if (ELECTRONS == 1) { - fstart[3] = 8; + fstart[3] = keln; // This is what actually selects what is taken from h5 file. hdf5_read_array(data[n]->p[KEL][0][0], "prims", 4, fdims, fstart, fcount, mdims, mstart, H5T_IEEE_F64LE); - fstart[3] = 9; + fstart[3] = 8; hdf5_read_array(data[n]->p[KTOT][0][0], "prims", 4, fdims, fstart, fcount, mdims, mstart, H5T_IEEE_F64LE); } diff --git a/model/iharm/model_params.h b/model/iharm/model_params.h index 0167fa7..3092878 100644 --- a/model/iharm/model_params.h +++ b/model/iharm/model_params.h @@ -12,8 +12,8 @@ #define B1 5 #define B2 6 #define B3 7 -#define KEL 8 -#define KTOT 9 +#define KTOT 8 +#define KEL 9 // These two will never be used simultaneously, // and never with KEL. #define TFLK 8 // temperature of fluid in Kelvin diff --git a/src/par.c b/src/par.c index a953e58..d9842d5 100644 --- a/src/par.c +++ b/src/par.c @@ -90,6 +90,8 @@ void load_par_from_argv(int argc, char *argv[], Params *params) { params->trace_i = -1; params->trace_j = -1; + params->keln = 9; // Kawazura + // I'm not sure there's still any advantage to "const" if we do this, // but hey, no warnings sscanf("trace.h5", "%s", (char *) (void *) params->trace_outf); @@ -187,6 +189,8 @@ void try_set_parameter(const char *word, const char *value, Params *params) { set_by_word_val(word, value, "trace_j", &(params->trace_j), TYPE_INT); set_by_word_val(word, value, "trace_outf", (void *)(params->trace_outf), TYPE_STR); + set_by_word_val(word, value, "keln", &(params->keln), TYPE_INT); + // Let models add/parse their own parameters we don't understand try_set_model_parameter(word, value); diff --git a/src/par.h b/src/par.h index b2c0da8..0d1ffa1 100644 --- a/src/par.h +++ b/src/par.h @@ -62,6 +62,8 @@ typedef struct params_t { int trace_stride; int trace_i, trace_j; const char trace_outf[STRLEN]; + + int keln; } Params; // modify this to set default values From 8fe4523fd64b3ff91ee6f02d4cd60b5ce8998f41 Mon Sep 17 00:00:00 2001 From: Cesar Diaz Date: Wed, 3 Nov 2021 14:59:05 -0500 Subject: [PATCH 2/4] keln dataype is now an integer and is under header/electrons --- model/iharm/model.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/model/iharm/model.c b/model/iharm/model.c index 878c9f8..9065eef 100644 --- a/model/iharm/model.c +++ b/model/iharm/model.c @@ -37,7 +37,6 @@ double RHO_unit; double U_unit; double B_unit; double Te_unit; -int keln; // MOLECULAR WEIGHTS static double Ne_factor = 1.; // e.g., used for He with 2 protons+neutrons per 2 electrons @@ -83,6 +82,7 @@ double tf; // TODO the way this is selected is horrid. Make it a parameter. #define ELECTRONS_TFLUID (3) static int RADIATION, ELECTRONS; +static int keln; static double gam = 1.444444, game = 1.333333, gamp = 1.666667; static double Thetae_unit, Mdotedd; @@ -1202,6 +1202,7 @@ void output_hdf5() hdf5_write_single_val(&mu_tot, "mu_tot", H5T_IEEE_F64LE); } hdf5_write_single_val(&ELECTRONS, "type", H5T_STD_I32LE); + hdf5_write_single_val(&keln, "keln", H5T_STD_I32LE); hdf5_set_directory("/header/"); hdf5_make_directory("units"); @@ -1210,7 +1211,6 @@ void output_hdf5() hdf5_write_single_val(&M_unit, "M_unit", H5T_IEEE_F64LE); hdf5_write_single_val(&T_unit, "T_unit", H5T_IEEE_F64LE); hdf5_write_single_val(&Te_unit, "Thetae_unit", H5T_IEEE_F64LE); - hdf5_write_single_val(&keln, "keln", H5T_IEEE_F64LE); hdf5_set_directory("/"); From 4d530f79dd118f9a252a36526a9e0c6475d34d31 Mon Sep 17 00:00:00 2001 From: Cesar Diaz Date: Fri, 10 Dec 2021 16:07:35 -0600 Subject: [PATCH 3/4] rhigh model is now default --- model/iharm/model.c | 25 ++++++++++++++++--------- src/par.c | 2 ++ src/par.h | 1 + 3 files changed, 19 insertions(+), 9 deletions(-) diff --git a/model/iharm/model.c b/model/iharm/model.c index 9065eef..14dca4c 100644 --- a/model/iharm/model.c +++ b/model/iharm/model.c @@ -19,7 +19,7 @@ #define NVAR (10) #define USE_FIXED_TPTE (0) -#define USE_MIXED_TPTE (0) // KEL models are actually used +// #define USE_MIXED_TPTE (1) // Now rhigh_model #define NSUP (3) #define USE_GEODESIC_SIGMACUT (1) @@ -83,6 +83,7 @@ double tf; #define ELECTRONS_TFLUID (3) static int RADIATION, ELECTRONS; static int keln; +static int rhigh_model; static double gam = 1.444444, game = 1.333333, gamp = 1.666667; static double Thetae_unit, Mdotedd; @@ -153,6 +154,7 @@ void try_set_model_parameter(const char *word, const char *value) // for selecting electron model set_by_word_val(word, value, "keln", &keln, TYPE_INT); + set_by_word_val(word, value, "rhigh_model", &rhigh_model, TYPE_INT); } // Advance through dumps until we are closer to the next set @@ -528,7 +530,11 @@ void set_units() fprintf(stderr,"MBH: %g [Msun]\n",MBH/MSUN); fprintf(stderr,"L,T,M units: %g [cm] %g [s] %g [g]\n",L_unit,T_unit,M_unit) ; fprintf(stderr,"rho,u,B units: %g [g cm^-3] %g [g cm^-1 s^-2] %g [G] \n",RHO_unit,U_unit,B_unit) ; - fprintf(stderr,"Model is: %d \n", keln) ; + if (rhigh_model == 1) { + fprintf(stderr,"Rhigh model is used \n") ; + } else { + fprintf(stderr,"Electron Model #%d is used\n", keln) ; + } } void init_physical_quantities(int n) @@ -687,17 +693,17 @@ void init_hamr_grid(char *fnam, int dumpidx) // we can override which electron model to use here. print results if we're // overriding anything. ELECTRONS should only be nonzero if we need to make // use of extra variables (instead of just UU and RHO) for thetae - if (!USE_FIXED_TPTE && !USE_MIXED_TPTE) { + if (!USE_FIXED_TPTE && !rhigh_model) { fprintf(stderr, "! no electron temperature model specified in model/iharm.c\n"); exit(-3); - } else if (USE_FIXED_TPTE && !USE_MIXED_TPTE) { + } else if (USE_FIXED_TPTE && !rhigh_model) { ELECTRONS = 0; // force TP_OVER_TE to overwrite bad electrons fprintf(stderr, "using fixed tp_over_te ratio = %g\n", tp_over_te); //Thetae_unit = MP/ME*(gam-1.)*1./(1. + tp_over_te); // see, e.g., Eq. 8 of the EHT GRRT formula list. // this formula assumes game = 4./3 and gamp = 5./3 Thetae_unit = 2./3. * MP/ME / (2. + tp_over_te); - } else if (USE_MIXED_TPTE && !USE_FIXED_TPTE) { + } else if (rhigh_model && !USE_FIXED_TPTE) { ELECTRONS = 2; fprintf(stderr, "using mixed tp_over_te with trat_small = %g, trat_large = %g, and beta_crit = %g\n", trat_small, trat_large, beta_crit); @@ -849,7 +855,7 @@ void init_iharm_grid(char *fnam, int dumpidx) // we can override which electron model to use here. print results if we're // overriding anything. ELECTRONS should only be nonzero if we need to make // use of extra variables (instead of just UU and RHO) for thetae - if (!USE_FIXED_TPTE && !USE_MIXED_TPTE) { + if (!USE_FIXED_TPTE && !rhigh_model) { if (ELECTRONS != 1) { fprintf(stderr, "! no electron temperature model specified! Cannot continue\n"); exit(-3); @@ -859,14 +865,14 @@ void init_iharm_grid(char *fnam, int dumpidx) } else if (ELECTRONS == ELECTRONS_TFLUID) { fprintf(stderr, "Using Ressler/Athena electrons with mixed tp_over_te and\n"); fprintf(stderr, "trat_small = %g, trat_large = %g, and beta_crit = %g\n", trat_small, trat_large, beta_crit); - } else if (USE_FIXED_TPTE && !USE_MIXED_TPTE) { + } else if (USE_FIXED_TPTE && !rhigh_model) { ELECTRONS = 0; // force TP_OVER_TE to overwrite bad electrons fprintf(stderr, "Using fixed tp_over_te ratio = %g\n", tp_over_te); //Thetae_unit = MP/ME*(gam-1.)*1./(1. + tp_over_te); // see, e.g., Eq. 8 of the EHT GRRT formula list. // this formula assumes game = 4./3 and gamp = 5./3 Thetae_unit = 2./3. * MP/ME / (2. + tp_over_te); - } else if (USE_MIXED_TPTE && !USE_FIXED_TPTE) { + } else if (rhigh_model && !USE_FIXED_TPTE) { ELECTRONS = 2; fprintf(stderr, "Using mixed tp_over_te with trat_small = %g, trat_large = %g, and beta_crit = %g\n", trat_small, trat_large, beta_crit); @@ -1144,7 +1150,7 @@ void init_koral_grid(char *fnam, int dumpidx) fprintf(stderr, "koral dump has native electron temperature. forcing Thetae...\n"); ELECTRONS = 9; } else { - if (USE_MIXED_TPTE && !USE_FIXED_TPTE) { + if (rhigh_model && !USE_FIXED_TPTE) { fprintf(stderr, "Using mixed tp_over_te with trat_small = %g, trat_large = %g, and beta_crit = %g\n", trat_small, trat_large, beta_crit); ELECTRONS = 2; } else { @@ -1203,6 +1209,7 @@ void output_hdf5() } hdf5_write_single_val(&ELECTRONS, "type", H5T_STD_I32LE); hdf5_write_single_val(&keln, "keln", H5T_STD_I32LE); + hdf5_write_single_val(&rhigh_model, "rhigh_model", H5T_STD_I32LE); hdf5_set_directory("/header/"); hdf5_make_directory("units"); diff --git a/src/par.c b/src/par.c index d9842d5..78ba352 100644 --- a/src/par.c +++ b/src/par.c @@ -91,6 +91,7 @@ void load_par_from_argv(int argc, char *argv[], Params *params) { params->trace_j = -1; params->keln = 9; // Kawazura + params->rhigh_model = 1; // R-high model // I'm not sure there's still any advantage to "const" if we do this, // but hey, no warnings @@ -190,6 +191,7 @@ void try_set_parameter(const char *word, const char *value, Params *params) { set_by_word_val(word, value, "trace_outf", (void *)(params->trace_outf), TYPE_STR); set_by_word_val(word, value, "keln", &(params->keln), TYPE_INT); + set_by_word_val(word, value, "rhigh_model", &(params->rhigh_model), TYPE_INT); // Let models add/parse their own parameters we don't understand try_set_model_parameter(word, value); diff --git a/src/par.h b/src/par.h index 0d1ffa1..c3bea69 100644 --- a/src/par.h +++ b/src/par.h @@ -64,6 +64,7 @@ typedef struct params_t { const char trace_outf[STRLEN]; int keln; + int rhigh_model; } Params; // modify this to set default values From 2d20d1fdd4621ff73a237ee5f70add00ea12ff20 Mon Sep 17 00:00:00 2001 From: Cesar Diaz Date: Fri, 10 Dec 2021 16:13:36 -0600 Subject: [PATCH 4/4] parameter rhigh_model different to 1 to use electron models instead of rhigh --- model/iharm/example.par | 1 + 1 file changed, 1 insertion(+) diff --git a/model/iharm/example.par b/model/iharm/example.par index 235aa19..c07ea0c 100644 --- a/model/iharm/example.par +++ b/model/iharm/example.par @@ -1,5 +1,6 @@ # 9 Kawazura, 10 Werner, 11 Rowan, 12 Sharma, default is 9 keln 10 +rhigh_model 0 # example.par -- ipole example parameter file for iharm problem