Merge pull request #200 from mabruzzo/HydrogenFractionByMass

Updating behavior of `HydrogenFractionByMass` parameter

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
 

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 */

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,59 @@ 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. 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",
+              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;