From 096489e1a33916070a047f9a44d0a2da0c237c95 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Sun, 28 Apr 2024 15:17:01 -0400 Subject: [PATCH 1/2] updating behavior of HydrogenFractionByMass parameter --- doc/source/Parameters.rst | 7 ++- src/clib/grackle_chemistry_data_fields.def | 12 ++++- src/clib/initialize_chemistry_data.c | 60 ++++++++++++++++++++-- 3 files changed, 71 insertions(+), 8 deletions(-) diff --git a/doc/source/Parameters.rst b/doc/source/Parameters.rst index 0876b0a9..7905442a 100644 --- a/doc/source/Parameters.rst +++ b/doc/source/Parameters.rst @@ -297,8 +297,11 @@ For all on/off integer flags, 0 is off and 1 is on. When using the non-equilibrium solver, a sensible default is 0.76. However, the tables for tabulated mode were created assuming n\ :sub:`He`/n\ :sub:`H` = 0.1, which corresponds to an H mass fraction of - about 0.716. When running in tabulated mode, this parameter will automatically - be changed to this value. Default: 0.76. + about 0.716. + + The default value stored in this variable is a negative value that denotes the value is unset. + If the user doesn't modify the value, the value is overwritten with a value of about 0.716 in tabulated mode and a value of 0.76 in non-equilibrium mode. + While users are allowed to set arbitrary values for the non-equilibrium solver, tabulated mode reports an error if the user initializes this to a value that does not exactly match the default. .. c:var:: float DeuteriumToHydrogenRatio diff --git a/src/clib/grackle_chemistry_data_fields.def b/src/clib/grackle_chemistry_data_fields.def index 61ae2011..bcc8aa4a 100644 --- a/src/clib/grackle_chemistry_data_fields.def +++ b/src/clib/grackle_chemistry_data_fields.def @@ -105,7 +105,17 @@ ENTRY(ih2co, INT, 1) /* flag for photoionization cooling */ ENTRY(ipiht, INT, 1) -ENTRY(HydrogenFractionByMass, DOUBLE, 0.76) +/* This is initially set to an undefined value. If the user doesn't + * specify a value, this is subsequently assigned a value based on + * the value of primordial chemistry: + * -> primordial_chemistry == 0: based on the values used for generating + * the table + * -> primordial_chemistry >= 1: 0.76 + * + * Historically, this was set to a default value of 0.76 and was ALWAYS + * silently overwritten when primordial_chemistry is equal to 0. + */ +ENTRY(HydrogenFractionByMass, DOUBLE, FLOAT_UNDEFINED) /* The DToHRatio is by mass in the code, so multiply by 2. - the default value comes from Burles & Tytler 1998 */ diff --git a/src/clib/initialize_chemistry_data.c b/src/clib/initialize_chemistry_data.c index c3b1126a..3f4ac172 100644 --- a/src/clib/initialize_chemistry_data.c +++ b/src/clib/initialize_chemistry_data.c @@ -213,7 +213,7 @@ int local_initialize_chemistry_data(chemistry_data *my_chemistry, "OPENMP\n"); return FAIL; } -# else _OPENMP +# else /* _OPENMP */ if (my_chemistry->omp_nthreads < 1) { // this is the default behavior (unless the user intervenes) my_chemistry->omp_nthreads = omp_get_max_threads(); @@ -238,13 +238,63 @@ int local_initialize_chemistry_data(chemistry_data *my_chemistry, return FAIL; } - if (my_chemistry->primordial_chemistry == 0) { + // deal with my_chemistry->HydrogenFractionByMass + // ============================================== + // - set_default_chemistry_parameters & local_initialize_chemistry_parameters + // historically configured my_chemistry such that this member had a value + // of 0.76 (the default when primordial_chemistry >= 1). If + // primordial_chemistry == 0, this would always silently overwrite the + // value in fully tabulated mode. This led to cases where downstream users + // might mistakenly think they modified this parameter + // + // - New behavior: the member is now initialized to an undefined value + // - if the user hasn't changed it, we now set it to the appropriate + // default value + // - if the user has changed it, the precise handling depends on the choice + // of primordial_chemistry + + if (my_chemistry->HydrogenFractionByMass > 1) { + fprintf(stderr, "ERROR: HydrogenFractionByMass cannot exceed 1.0\n"); + return FAIL; + } else if (my_chemistry->primordial_chemistry == 0) { /* In fully tabulated mode, set H mass fraction according to the abundances in Cloudy, which assumes n_He / n_H = 0.1. - This gives a value of about 0.716. Using the default value + This gives a value of about 0.716. Using the other default value of 0.76 will result in negative electron densities at low - temperature. Below, we set X = 1 / (1 + m_He * n_He / n_H). */ - my_chemistry->HydrogenFractionByMass = 1. / (1. + 0.1 * 3.971); + temperature. Below, use X = 1 / (1 + m_He * n_He / n_H). */ + + // we precomputed this value (in enhanced precision) rather than compute + // it on the fly to allow users to directly specify this value (without + // worrying about round-off issues) + const double default_Hfrac = 0.715768377353088514; + // TODO: at some point in the future, we should store the above value as + // an attribute of the HDF5 table and read it in directly. + + if (my_chemistry->HydrogenFractionByMass < 0) { + my_chemistry->HydrogenFractionByMass = default_Hfrac; + } else if (my_chemistry->HydrogenFractionByMass != default_Hfrac) { + fprintf(stderr, + "ERROR: Invalid HydrogenFractionByMass value is specified.\n" + " -> when primordial_chemistry == 0, the allowed values are\n" + " are strictly enforced.\n" + " -> the specified value is: %.18f\n" + " -> it is allowed to have a negative value. In this case, the\n" + " value is overwritten to match the tabulated default\n" + " value of %.18f that is consistent with assumptions made\n" + " while constructing the table (this is the default \n" + " behavior if you don't touch the value of this member.)\n" + " -> alternatively you can EXACTLY specify the default value\n" + " -> NOTE: for primordial_chemistry == 0, prior versions of\n" + " Grackle would silently overwrite the value of\n" + " HydrogenFractionByMass instead of reporting this error\n", + my_chemistry->HydrogenFractionByMass, default_Hfrac); + return FAIL; + } + } else { + const double default_Hfrac = 0.76; + if (my_chemistry->HydrogenFractionByMass < 0) { + my_chemistry->HydrogenFractionByMass = default_Hfrac; + } } double co_length_units, co_density_units; From b7c7d71933ddd722676323ca9af5cad603caa9b1 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 19 Jun 2024 10:54:07 -0400 Subject: [PATCH 2/2] address Britton's concerns --- src/clib/initialize_chemistry_data.c | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/src/clib/initialize_chemistry_data.c b/src/clib/initialize_chemistry_data.c index 3f4ac172..cf349a8c 100644 --- a/src/clib/initialize_chemistry_data.c +++ b/src/clib/initialize_chemistry_data.c @@ -276,18 +276,14 @@ int local_initialize_chemistry_data(chemistry_data *my_chemistry, fprintf(stderr, "ERROR: Invalid HydrogenFractionByMass value is specified.\n" " -> when primordial_chemistry == 0, the allowed values are\n" - " are strictly enforced.\n" - " -> the specified value is: %.18f\n" - " -> it is allowed to have a negative value. In this case, the\n" - " value is overwritten to match the tabulated default\n" - " value of %.18f that is consistent with assumptions made\n" - " while constructing the table (this is the default \n" - " behavior if you don't touch the value of this member.)\n" - " -> alternatively you can EXACTLY specify the default value\n" + " are strictly enforced. You have 2 options: \n" + " 1. leave the value unset or set it to a negative number\n" + " to have it calculated internally as %.18f\n" + " 2. set the value exactly to the above number\n" " -> NOTE: for primordial_chemistry == 0, prior versions of\n" " Grackle would silently overwrite the value of\n" " HydrogenFractionByMass instead of reporting this error\n", - my_chemistry->HydrogenFractionByMass, default_Hfrac); + default_Hfrac); return FAIL; } } else {