Modified default compilation mode to not use -ffast-math

It turns out that this was the reason certain summary statistics were varying by more than 0.1% when I was using 2 different machines with different compiler versions (gcc versions 9.4.0 vs 7.5.0). This was not completely obvious because Grackle itself was never compiled with these flags (although it took me longer to determine than it should have). As far as I can tell most of the differences arose in EnzoInitialGrackleTest.

We were already making this change in the CMake PR (#139), anyways, so this shouldn't be controversial. As a result of this change, the following modifications were required for existing tests:
- the expected values for the VL+CT dual_energy_cloud test needed to be updated. This change was previously discussed in the context of PR #139.
- the final time for the PPML tests had to be updated for binaries compiled with single precision
- the cosmology test had to be disabled for binaries compiled with single precision. As is described in Issue #154, it appears the `-ffast-math` had been suppressing an existing bug for this case (the test now crashes without the flag). Note that the test ran perfectly fine when compiled in double-precision.

The double-precision results for the Grackle test also had to be updated (but the single-precision results were basically fine). Note that the tolerances for the double-precision results were calibrated so that the tests would pass if you used gcc versions 7.5.0 or 9.4.0 (they may need to be loosened more in the future).

This commit also added the ability for the field_summary to accept different tolerances for different summary statistics (this is useful since we might want 0 tolerance for the location of a min/max while still applying a tolerance for the min/max itself). It also a runtime warnings raised by the script.

config/linux_gnu.py

@@ -6,7 +6,7 @@
 
 #
 #flags_arch = '-g -fprofile-arcs -ftest-coverage' # gcov
-flags_arch = '-O3 -g -ffast-math -funroll-loops -fPIC -pedantic'
+flags_arch = '-O3 -g -funroll-loops -fPIC -pedantic' # + ' -ffast-math'
 #flags_arch = '-Wall -O3 -g'
 #flags_arch = '-Wall -O0 -g'
 

input/PPML/ppml.incl

@@ -103,8 +103,8 @@
 
    Testing
    {
-      time_final  = [0.00294580065398299,
-                     0.00294612924335524 ];
+      time_final  = [0.0029466382402461,    # single-precision
+                     0.00294612924335524 ]; # double-precision
       time_tolerance = 1.0e-4;
    }
 

input/vlct/run_dual_energy_cloud_test.py

@@ -81,7 +81,7 @@ def check_cloud_asym(fname, name, max_asym):
 def analyze_tests():
     r = []
     r += check_cloud_asym('hlld_cloud_0.0625/hlld_cloud_0.0625.block_list',
-                          'hlld_cloud', 5.5e-13)
+                          'hlld_cloud', 5.9e-13)
     r += check_cloud_asym('hllc_cloud_0.0625/hllc_cloud_0.0625.block_list',
                           'hllc_cloud', 3.8e-13)
     r += check_cloud_asym('hlle_cloud_0.0625/hlle_cloud_0.0625.block_list',

test/MethodCosmology/SConscript

@@ -5,6 +5,7 @@ Import('ip_charm')
 
 Import('bin_path')
 Import('test_path')
+Import('prec')
 
 #----------------------------------------------------------
 #defines
@@ -37,35 +38,38 @@ env_mv_cosmology_8 = env.Clone(COPY = 'mkdir -p ' + test_path + '/MethodCosmolog
 #-------------------------------------------------------------
 #load balancing
 #-------------------------------------------------------------
+if prec == 'double':
+    # as documented in ISSUE #154, there are some issues with running cosmology
+    # tests in single-precision mode. This test has been disabled for now.
 
-#serial
-balance_cosmology_1 = env_mv_cosmology_1.RunCosmology_1 (
-     'test_method_cosmology-1.unit',
-     bin_path + '/enzo-e',
-     ARGS='input/Cosmology/method_cosmology-1.in')
+    #serial
+    balance_cosmology_1 = env_mv_cosmology_1.RunCosmology_1 (
+         'test_method_cosmology-1.unit',
+         bin_path + '/enzo-e',
+         ARGS='input/Cosmology/method_cosmology-1.in')
 
-Clean(balance_cosmology_1,
-     [Glob('#/' + test_path + '/Dir_COSMO-1'),
-      Glob('#/' + test_path + '/Dir_COSMO-1')])
+    Clean(balance_cosmology_1,
+         [Glob('#/' + test_path + '/Dir_COSMO-1'),
+          Glob('#/' + test_path + '/Dir_COSMO-1')])
 
-#env.MakeMovie("method_cosmology-1.swf", "test_method_cosmology-1.unit", \
-#              ARGS = test_path + '/method_cosmology-1*.png");
-#env.PngToGif("method_cosmology-1.gif", "test_method_cosmology-1.unit", \
-#              ARGS = test_path + '/method_cosmology-1*.png");
+    #env.MakeMovie("method_cosmology-1.swf", "test_method_cosmology-1.unit", \
+    #              ARGS = test_path + '/method_cosmology-1*.png");
+    #env.PngToGif("method_cosmology-1.gif", "test_method_cosmology-1.unit", \
+    #              ARGS = test_path + '/method_cosmology-1*.png");
 
 
-#parallel
+    #parallel
 
-balance_cosmology_8 = env_mv_cosmology_8.RunCosmology_8 (
-     'test_method_cosmology-8.unit',
-     bin_path + '/enzo-e',
-     ARGS='input/Cosmology/method_cosmology-8.in')
+    balance_cosmology_8 = env_mv_cosmology_8.RunCosmology_8 (
+         'test_method_cosmology-8.unit',
+         bin_path + '/enzo-e',
+         ARGS='input/Cosmology/method_cosmology-8.in')
 
-Clean(balance_cosmology_8,
-     [Glob('#/' + test_path + '/Dir_COSMO-8'),
-      Glob('#/' + test_path + '/Dir_COSMO-8')])
+    Clean(balance_cosmology_8,
+         [Glob('#/' + test_path + '/Dir_COSMO-8'),
+          Glob('#/' + test_path + '/Dir_COSMO-8')])
 
-#env.MakeMovie("method_cosmology-8.swf", "test_method_cosmology-8.unit", \
-#              ARGS = test_path + '/method_cosmology-8*.png");
-#env.PngToGif("method_cosmology-8.gif", "test_method_cosmology-8.unit", \
-#              ARGS = test_path + '/method_cosmology-8*.png");
+    #env.MakeMovie("method_cosmology-8.swf", "test_method_cosmology-8.unit", \
+    #              ARGS = test_path + '/method_cosmology-8*.png");
+    #env.PngToGif("method_cosmology-8.gif", "test_method_cosmology-8.unit", \
+    #              ARGS = test_path + '/method_cosmology-8*.png");

test/MethodGrackle/SConscript

@@ -137,12 +137,21 @@ if grackle_data_dir != '' and use_grackle:
 
     # Part 3 of 3: Measure the field summary statistics and compare against
     # a table of reference values
+    _dict_to_json_cmd_arg = lambda d: "'" + str(d).replace('\'', '"') + "'"
     if prec == 'double':
         _ref_tab = test_path + '/MethodGrackle/ref_general_grackle-double.csv'
-        atol, rtol = 1e-15, 1e-14
+        atol = 0
+        # these should be flexible for different compiler versions
+        rtol = _dict_to_json_cmd_arg({"min" : 5e-15, "max" : 5e-6,
+                                      "mean" : 5e-8,
+                                      "standard_deviation" : 5e-8})
     else:
         _ref_tab = test_path + '/MethodGrackle/ref_general_grackle-single.csv'
-        atol, rtol = 1e-7, 1e-6
+        atol = 0
+        # the following may need to be relaxed for different compiler versions
+        rtol = _dict_to_json_cmd_arg(
+            dict((k, 1e-7) for k in ["min","max","mean","standard_deviation"])
+        )
 
     _compare_log_path = test_path + '/MethodGrackle/grackle_general_cmp.log'
     general_grackle = env.CompareFieldSummary(

test/MethodGrackle/ref_general_grackle-double.csv

@@ -1,6 +1,6 @@
 # {"code unit definitions": {"code_length": [3.086e+16, "cm"], "code_mass": [4.9157019637146824e+25, "g"], "code_density": [1.6726219e-24, "g/cm**3"], "code_specific_energy": [956277.4378779504, "erg/g"], "code_time": [31557600000000.0, "s"], "code_magnetic": [4.483279099741026e-09, "G"], "code_temperature": [1.0, "K"], "code_pressure": [1.5994905850705492e-18, "dyn/cm**2"], "code_velocity": [977.8943899409335, "cm/s"], "code_metallicity": [1.0, "dimensionless"]}, "field units": {"pressure": "dimensionless", "temperature": "code_temperature", "cooling_time": "dimensionless"}}
 #
 # name,min,min_xloc,min_yloc,min_zloc,max,max_xloc,max_yloc,max_zloc,mean,standard_deviation
-pressure,1.520195901060065e+02,7.812500000000000e-02,1.718750000000000e-01,1.562500000000000e-02,2.753947492701502e+06,7.812500000000000e-02,2.343750000000000e-01,1.562500000000000e-02,6.675724774669249e+04,3.557197556605695e+05
-temperature,3.549497639185964e+00,2.343750000000000e-01,1.718750000000000e-01,1.562500000000000e-02,6.376232820587005e+06,1.562500000000000e-02,2.343750000000000e-01,1.562500000000000e-02,2.417771901546382e+05,1.115430124658054e+06
-cooling_time,-6.450617403727277e+04,1.562500000000000e-02,2.343750000000000e-01,1.562500000000000e-02,3.866436169424824e+03,1.093750000000000e-01,1.093750000000000e-01,1.562500000000000e-02,-1.578969463027807e+03,8.485129947191494e+03
+pressure,1.520195900675427e+02,7.812500000000000e-02,1.718750000000000e-01,1.562500000000000e-02,2.753947492701290e+06,7.812500000000000e-02,2.343750000000000e-01,1.562500000000000e-02,6.660831661417476e+04,3.557307908963668e+05
+temperature,3.330949646514823e+00,2.343750000000000e-01,1.718750000000000e-01,1.562500000000000e-02,6.376232820586985e+06,1.562500000000000e-02,2.343750000000000e-01,1.562500000000000e-02,2.417771759620777e+05,1.115430127734497e+06
+cooling_time,-6.450617403727239e+04,1.562500000000000e-02,2.343750000000000e-01,1.562500000000000e-02,3.866436114544078e+03,1.093750000000000e-01,1.093750000000000e-01,1.562500000000000e-02,-1.576307136350389e+03,8.485636211531526e+03

tools/field_summary.py

@@ -27,6 +27,8 @@
     ('standard_deviation', 'f8')
 ]
 
+_TABLE_COLS = frozenset(col for col, _ in _TABLE_DTYPE_ARGS)
+
 _FIXED_COL_UNITS = dict((col, 'code_length') for col,_ in _TABLE_DTYPE_ARGS \
                         if col[-5:] in ['_xloc', '_yloc', '_zloc'])
 
@@ -168,9 +170,7 @@ def test_equivalent_field_tables(cur_field_table, ref_field_table,
     else:
         tr = test_report
 
-    tr.write(
-        'Comparing Field Tables (atol = {}, rtol = {})\n'.format(atol,rtol)
-    )
+    tr.write('Comparing Field Tables\n')
 
     # first check that table properties are consistent
     consistent_table_prop = _check_consistent_cols_and_names(cur_field_table,
@@ -184,18 +184,24 @@ def test_equivalent_field_tables(cur_field_table, ref_field_table,
         'General table properties are consistent. Proceeding with comparison\n'
     )
 
-    # Now actually verify equality of properties
     colnames = tuple(ref_field_table.dtype.fields.keys())
+
+    filtered_colnames = list(filter(lambda col: col != 'name', colnames))
+    for tol, tolname in [(atol, 'atol'), (rtol, 'rtol')]:
+        tr.write(f'{tolname} = {tol.represent_as_string(filtered_colnames)}\n')
+
+    # Now actually verify equality of properties
     comparison_arr = np.empty((len(cur_field_table), len(colnames)),
-                              dtype = np.bool)
+                              dtype = np.bool_)
     for j,col in enumerate(colnames):
         if col == "name":
             comparison_arr[:,j] = True
         else:
             cur_vals = cur_field_table[col]
             ref_vals = ref_field_table[col]
-            comparison_arr[:,j] = np.isclose(cur_vals, ref_vals, rtol = rtol,
-                                             atol = atol, equal_nan = False)
+            comparison_arr[:,j] = np.isclose(cur_vals, ref_vals,
+                                             rtol = rtol[col], atol = atol[col],
+                                             equal_nan = False)
 
     all_consistent_vals = True
     # Finally check comparison results on a row-by-row basis
@@ -249,6 +255,53 @@ def _main_measure(args):
         print(output.getvalue())
     output.close()
 
+class ToleranceConfig:
+    def __init__(self, fallback_tol, col_specific_vals):
+        self._fallback_tol = fallback_tol
+        self._col_specific_vals = col_specific_vals
+
+    def __getitem__(self, key):
+        return self._col_specific_vals.get(key, self._fallback_tol)
+
+    def represent_as_string(self, colnames):
+        if len(self._col_specific_vals) == 0:
+            return str(self._fallback_tol)
+        else:
+            return str(dict((col,self[col]) for col in colnames))
+
+def _process_tol_args(arg_namespace):
+    # Process the tolerance arguments. They can be a string encoding an
+    # - int or float (this applies to all columns)
+    # - JSON object that assoicates tolerances with columns of the summary
+    #   table (columns without entries have a tolerance of 0)
+
+    def _process(tol_name):
+        try:
+            tmp = json.loads(getattr(arg_namespace, tol_name))
+        except json.JSONDecodeError:
+            tmp = None
+
+        if isinstance(tmp, (float, int)):
+            return ToleranceConfig(fallback_tol = tmp, col_specific_vals = {})
+
+        elif isinstance(tmp, dict):
+            for col, val in tmp.items(): # check contents of dict
+                if col not in _TABLE_COLS:
+                    raise ValueError(f"{tol_name} specified for unknown table "
+                                     f"column: {col}")
+                elif not isinstance(val, (int,float)):
+                    raise ValueError(f"{tol_name} for '{col}', {val}, isn't "
+                                     "an int or float")
+            return ToleranceConfig(fallback_tol = 0., col_specific_vals = tmp)
+
+        else:
+            raise ValueError(
+                f"{tol_name} option expects an int/float or a JSON object "
+                "that pairs summary table column-names with int/floats. \n"
+                f"Received: '{getattr(arg_namespace, tol_name)}'"
+            )
+    return _process('atol'), _process('rtol')
+
 
 def _main_cmp(args):
     # Program to use by the cmp subcommand. Just construct the field summary
@@ -257,6 +310,8 @@ def _main_cmp(args):
           "summarized")
     ref_field_table = read_field_summary(args.ref)
 
+    atol,rtol = _process_tol_args(args)
+
     field_names = ref_field_table['name'].tolist()
     print("Measuring the Field Summary Properties")
     cur_field_table, sim_props = measure_field_summary(
@@ -272,8 +327,7 @@ def _main_cmp(args):
         print("Comparing field summary tables")
         test_rslt = test_equivalent_field_tables(cur_field_table,
                                                  ref_field_table,
-                                                 atol = args.atol,
-                                                 rtol = args.rtol,
+                                                 atol = atol, rtol = rtol,
                                                  test_report = tr)
         if test_rslt:
             print("Field summary tables are consistent")
@@ -327,13 +381,21 @@ def _add_target_path_arg(parser):
     help = ('Path to file where a properly formatted test report describing '
             'the sucess of the comparison should be written.')
 )
+
+_extended_explanation = (
+    'This expects a single tolerance value (an integer/floating point value) '
+    'that is used for all comparing values from any column of the '
+    'summary-table. Alternatively, this can accept a JSON object that pairs '
+    'column-names with the individual tolerance values. When a tolerance '
+    'valueIf a tolerance value isn\'t provided, it defaults to 0.'
+)
 cmp_parser.add_argument(
-    '--rtol', required = True, action = 'store', type = float,
-    help = "Relative tolerance"
+    '--rtol', action = 'store', default = '0',
+    help = f"Relative tolerance. {_extended_explanation}"
 )
 cmp_parser.add_argument(
-    '--atol', required = True, action = 'store', type = float,
-    help = "Absolute tolerance"
+    '--atol', action = 'store', default = '0',
+    help = f"Absolute tolerance. {_extended_explanation}"
 )
 
 cmp_parser.set_defaults(func = _main_cmp)