Add an ISMIP7 melt-module calibration test group

.flake8.cfg

@@ -14,3 +14,7 @@ exclude =
     .idea,
     .mypy_cache,
     .pytest_cache,
+    # a verbatim copy of the upstream ISMIP7 parameter-selection toolbox; it
+    # must not be modified, so it is not linted.  See PROVENANCE.md alongside
+    # it.
+    compass/landice/tests/ismip7_calibration/toolbox/parameter_selection_toolbox.py,

.github/workflows/build_workflow.yml

@@ -133,3 +133,4 @@ jobs:
           compass clean --help
           create_compass_load_script --help
           pip check
+          pytest -v

compass/landice/__init__.py

@@ -11,6 +11,7 @@
 from compass.landice.tests.hydro_radial import HydroRadial
 from compass.landice.tests.ismip6_forcing import Ismip6Forcing
 from compass.landice.tests.ismip6_run import Ismip6Run
+from compass.landice.tests.ismip7_calibration import Ismip7Calibration
 from compass.landice.tests.ismip7_forcing import Ismip7Forcing
 from compass.landice.tests.ismip7_run import Ismip7Run
 from compass.landice.tests.isunnguata_sermia import IsunnguataSermia
@@ -48,6 +49,7 @@ def __init__(self):
         self.add_test_group(HydroRadial(mpas_core=self))
         self.add_test_group(Ismip6Forcing(mpas_core=self))
         self.add_test_group(Ismip6Run(mpas_core=self))
+        self.add_test_group(Ismip7Calibration(mpas_core=self))
         self.add_test_group(Ismip7Forcing(mpas_core=self))
         self.add_test_group(Ismip7Run(mpas_core=self))
         self.add_test_group(IsunnguataSermia(mpas_core=self))

compass/landice/ismip7/__init__.py

@@ -0,0 +1,9 @@
+"""
+Shared framework code for the ISMIP7 test groups.
+
+The ISMIP7 test groups -- ``ismip7_forcing``, ``ismip7_run`` and
+``ismip7_calibration`` -- all remap data from the ISMIP7 polar stereographic
+grids onto a MALI mesh, and all need the same handful of helpers to do it.
+Those helpers live here rather than in any one test group, so that using them
+from another does not mean importing across test-group boundaries.
+"""

compass/landice/tests/ismip7_forcing/create_mapfile.py → compass/landice/ismip7/mapping.py

@@ -7,7 +7,7 @@
 
 def build_mapping_file(config, logger, ismip7_grid_file,
                        mapping_file, mali_mesh_file=None,
-                       method_remap=None, projection=None):
+                       method_remap=None, projection=None, ntasks=None):
     """
     Build a mapping file for regridding from an ISMIP7 polar
     stereographic grid to the MALI unstructured mesh.
@@ -35,6 +35,11 @@ def build_mapping_file(config, logger, ismip7_grid_file,
     projection : str, optional
         Projection flag for SCRIP generation (e.g., 'ais-bedmap2',
         'gis-bamber'). If not provided, reads from ice_sheet_params.
+
+    ntasks : int, optional
+        Number of MPI tasks to use for ESMF_RegridWeightGen. If not
+        provided, reads ``esmf_ntasks`` from the ``[ismip7]`` config
+        section.
     """
 
     if os.path.exists(mapping_file):
@@ -54,9 +59,7 @@ def build_mapping_file(config, logger, ismip7_grid_file,
 
     # Determine projection from parameter or config
     if projection is None:
-        from compass.landice.tests.ismip7_forcing.ice_sheet_params import (
-            get_params,
-        )
+        from compass.landice.ismip7.ice_sheet_params import get_params
         projection = get_params(config)['projection']
 
     ismip7_projection = projection
@@ -77,8 +80,13 @@ def build_mapping_file(config, logger, ismip7_grid_file,
 
     # create a MALI mesh scrip file
     logger.info("Creating SCRIP file for MALI mesh...")
+    # copy the contents but not the mode bits: the copy has lat/lon fields
+    # written into it below, and a mesh distributed read-only would
+    # otherwise produce a read-only copy
     mali_mesh_copy = f"{mali_mesh_file}_copy"
-    shutil.copy(mali_mesh_file, mali_mesh_copy)
+    if os.path.exists(mali_mesh_copy):
+        os.remove(mali_mesh_copy)
+    shutil.copyfile(mali_mesh_file, mali_mesh_copy)
 
     args = ["set_lat_lon_fields_in_planar_grid",
             "--file", mali_mesh_copy,
@@ -91,12 +99,12 @@ def build_mapping_file(config, logger, ismip7_grid_file,
     # create a mapping file using ESMF_RegridWeightGen
     logger.info(f"Creating mapping file with method: {method_remap}")
 
-    section = config["ismip7"]
-    cores = section.getint("esmf_ntasks")
+    if ntasks is None:
+        ntasks = config.getint("ismip7", "esmf_ntasks")
 
     parallel_executable = config.get("parallel", "parallel_executable")
     args = parallel_executable.split(" ")
-    args.extend(["-n", f"{cores}",
+    args.extend(["-n", f"{ntasks}",
                  "ESMF_RegridWeightGen",
                  "-s", source_grid_scripfile,
                  "-d", mali_scripfile,

compass/landice/tests/ismip7_forcing/fracture/remap_utils.py → compass/landice/ismip7/remap.py

@@ -1,5 +1,5 @@
 """
-Shared helpers for remapping ISMIP7 fracture forcing data to the MALI mesh.
+Shared helpers for remapping ISMIP7 forcing data onto the MALI mesh.
 """
 import os
 
@@ -9,7 +9,8 @@
 from scipy.ndimage import distance_transform_edt
 
 
-def extrapolate_source(input_file, output_file, varnames, logger):
+def extrapolate_source(input_file, output_file, varnames, logger,
+                       decode_times=True):
     """
     Extrapolate fill/missing values on the source polar stereographic grid
     using nearest-neighbor via ``distance_transform_edt``. This must be done
@@ -29,14 +30,20 @@ def extrapolate_source(input_file, output_file, varnames, logger):
 
     logger : logging.Logger
         Logger for status messages
+
+    decode_times : bool, optional
+        Whether to let xarray decode the time coordinate.  The fracture
+        forcing files use ``units="year"`` (integer years), which is not
+        CF-compliant, so those callers must pass ``False``.
     """
     if isinstance(varnames, str):
         varnames = [varnames]
 
     logger.info(f"    Extrapolating fill values on source grid: "
                 f"{os.path.basename(input_file)}")
 
-    ds = xr.open_dataset(input_file, decode_times=False)
+    ds = xr.open_dataset(input_file, engine="netcdf4",
+                         decode_times=decode_times)
 
     for varname in varnames:
         data = ds[varname]

compass/landice/tests/ismip7_calibration/__init__.py

@@ -0,0 +1,30 @@
+from compass.landice.tests.ismip7_calibration.ais import Ais
+from compass.landice.tests.ismip7_calibration.replication import Replication
+from compass.testgroup import TestGroup
+
+
+class Ismip7Calibration(TestGroup):
+    """
+    A test group for calibrating MALI's sub-shelf melt parameterization
+    against the ISMIP7 Antarctic ice-ocean protocol (Reese et al., Sect. 4.2)
+
+    The protocol asks each ice-sheet model to calibrate the free parameter of
+    its melt module against four objective-function terms -- basin-integrated
+    present-day melt, melt by buttressing bin, the warm-minus-cold sensitivity
+    of ocean models, and observed Amundsen ice-shelf melt -- and to report the
+    5th, 50th and 95th percentiles of the resulting parameter distribution.
+    """
+
+    def __init__(self, mpas_core):
+        """
+        Create the test group
+
+        Parameters
+        ----------
+        mpas_core : compass.landice.Landice
+            the MPAS core that this test group belongs to
+        """
+        super().__init__(mpas_core=mpas_core, name='ismip7_calibration')
+
+        self.add_test_case(Replication(test_group=self))
+        self.add_test_case(Ais(test_group=self))

compass/landice/tests/ismip7_calibration/ais/__init__.py

@@ -0,0 +1,158 @@
+from compass.landice.tests.ismip7_calibration import datasets
+from compass.landice.tests.ismip7_calibration.ais.aggregate import Aggregate
+from compass.landice.tests.ismip7_calibration.ais.calibrate import Calibrate
+from compass.landice.tests.ismip7_calibration.ais.fit_delta_t import FitDeltaT
+from compass.landice.tests.ismip7_calibration.ais.make_graph import MakeGraph
+from compass.landice.tests.ismip7_calibration.ais.remap_forcing import (
+    RemapForcing,
+)
+from compass.landice.tests.ismip7_calibration.ais.remap_masks import RemapMasks
+from compass.landice.tests.ismip7_calibration.ais.report import Report
+from compass.landice.tests.ismip7_calibration.ais.run_state import RunState
+from compass.landice.tests.ismip7_calibration.ais.verify_melt import VerifyMelt
+from compass.landice.tests.ismip7_calibration.configure import (
+    check_options,
+    melt_forms,
+)
+from compass.testcase import TestCase
+from compass.validate import compare_variables
+
+#: the ocean state the verification and the dT_b fit use
+REFERENCE_STATE = 'climatology'
+
+#: multiples of the reference melt parameter used to measure the linearity
+#: that the one-run-per-ocean-state ensemble relies on.  Only the 'ismip7'
+#: form is scaled: its parameter is a namelist option, while the ISMIP6
+#: gamma0 is read from an input file.
+LINEARITY_SCALES = (0.5, 2.0)
+
+
+class Ais(TestCase):
+    """
+    A test case that calibrates MALI's sub-shelf melt parameterization on an
+    Antarctic MALI mesh, following the ISMIP7 protocol.
+
+    The steps, in dependency order:
+
+    ``make_graph``
+        The graph partition the melt diagnostics run on, built once from the
+        mesh rather than required as a pre-built file alongside it.
+
+    ``remap_masks``
+        ISMIP7 basins, buttressing bins, the floating mask and the
+        PIG/Dotson regions, onto the MALI mesh.
+
+    ``remap_forcing``
+        The calibration thermal forcing for each ocean state, onto the MALI
+        mesh.
+
+    ``<melt_form>_<state>``
+        One single-timestep MALI melt diagnostic per ocean state and melt
+        form.  **Not** one per parameter value: melt is exactly proportional
+        to the melt parameter, so the parameter sweep is a scaling of one
+        run.  That is what makes this 28 runs per form rather than about
+        1300.
+
+    ``verify_melt``
+        MALI's melt against an independent Python implementation, its
+        vertical interpolation against an independent one, and the linearity
+        that the previous point relies on.
+
+    ``aggregate``
+        Melt to basins, buttressing bins and shelf regions, area-weighted.
+
+    ``calibrate``
+        The 100,000-sample parameter selection, giving the percentiles the
+        ISMIP7 projections need.
+
+    ``fit_delta_t``
+        The per-basin correction dT_b, fitted **after** parameter selection
+        per protocol Sect. 4.2.1 option 2.
+
+    ``report``
+        Plots and a summary table.
+
+    Attributes
+    ----------
+    melt_forms : list of str
+        The melt forms being calibrated
+
+    states : list of compass.landice.tests.ismip7_calibration.datasets.OceanState
+        The ocean states in the ensemble
+    """  # noqa: E501
+
+    def __init__(self, test_group):
+        """
+        Create the test case
+
+        Parameters
+        ----------
+        test_group : compass.landice.tests.ismip7_calibration.Ismip7Calibration
+            The test group that this test case belongs to
+        """  # noqa: E501
+        name = 'ais'
+        super().__init__(test_group=test_group, name=name, subdir=name)
+        self.melt_forms = []
+        self.states = []
+
+    def configure(self):
+        """
+        Add a step per ocean state and melt form, once the config is known
+        """
+        config = self.config
+        check_options(config, ['base_path_ismip7', 'base_path_mali',
+                               'mali_mesh_file', 'mali_mesh_name'])
+
+        section = config['ismip7_calibration']
+        base_path = section.get('base_path_ismip7')
+        subset = section.get('ocean_state_subset')
+
+        self.melt_forms = melt_forms(config)
+        self.states = datasets.ocean_states(base_path, subset=subset)
+
+        self.add_step(MakeGraph(test_case=self))
+        self.add_step(RemapMasks(test_case=self))
+        self.add_step(RemapForcing(test_case=self))
+
+        for melt_form in self.melt_forms:
+            for state in self.states:
+                self.add_step(RunState(
+                    test_case=self, state_name=state.name,
+                    melt_form=melt_form,
+                    subdir=f'{melt_form}_{state.name}'))
+
+        # extra runs of one ocean state at other melt parameters, so that the
+        # linearity the ensemble design relies on is measured in MALI rather
+        # than argued from the code
+        verify_form = self.melt_forms[0]
+        scales = LINEARITY_SCALES if verify_form == 'ismip7' else ()
+        for scale in scales:
+            self.add_step(RunState(
+                test_case=self, state_name=REFERENCE_STATE,
+                melt_form=verify_form,
+                subdir=f'{verify_form}_{REFERENCE_STATE}_x{scale:g}',
+                parameter_scale=scale))
+
+        self.add_step(VerifyMelt(test_case=self, melt_form=verify_form,
+                                 state_name=REFERENCE_STATE,
+                                 linearity_scales=scales))
+        self.add_step(Aggregate(test_case=self, melt_forms=self.melt_forms,
+                                states=self.states))
+        self.add_step(Calibrate(test_case=self, melt_forms=self.melt_forms))
+        self.add_step(FitDeltaT(test_case=self, melt_forms=self.melt_forms,
+                                state_name=REFERENCE_STATE))
+        self.add_step(Report(test_case=self, melt_forms=self.melt_forms))
+
+    def validate(self):
+        """
+        Compare the calibration against a baseline, if one was provided
+        """
+        variables = ['p5', 'median', 'p95', 'mode']
+        for melt_form in self.melt_forms:
+            compare_variables(
+                test_case=self, variables=variables,
+                filename1=f'calibrate/calibration_{melt_form}.nc')
+        compare_variables(
+            test_case=self,
+            variables=['modelled_shelf_area', 'observed_shelf_area'],
+            filename1='aggregate/shelf_area.nc')