diff --git a/.flake8.cfg b/.flake8.cfg index f7ad654c10..39769d45a6 100644 --- a/.flake8.cfg +++ b/.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, diff --git a/.github/workflows/build_workflow.yml b/.github/workflows/build_workflow.yml index 309b19fe97..cafc30c184 100644 --- a/.github/workflows/build_workflow.yml +++ b/.github/workflows/build_workflow.yml @@ -133,3 +133,4 @@ jobs: compass clean --help create_compass_load_script --help pip check + pytest -v diff --git a/compass/landice/__init__.py b/compass/landice/__init__.py index aafe35fdd6..0cb34ca707 100644 --- a/compass/landice/__init__.py +++ b/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)) diff --git a/compass/landice/ismip7/__init__.py b/compass/landice/ismip7/__init__.py new file mode 100644 index 0000000000..dd68fbf14a --- /dev/null +++ b/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. +""" diff --git a/compass/landice/tests/ismip7_forcing/ice_sheet_params.py b/compass/landice/ismip7/ice_sheet_params.py similarity index 100% rename from compass/landice/tests/ismip7_forcing/ice_sheet_params.py rename to compass/landice/ismip7/ice_sheet_params.py diff --git a/compass/landice/tests/ismip7_forcing/create_mapfile.py b/compass/landice/ismip7/mapping.py similarity index 81% rename from compass/landice/tests/ismip7_forcing/create_mapfile.py rename to compass/landice/ismip7/mapping.py index de10f32b34..d99c79f01e 100644 --- a/compass/landice/tests/ismip7_forcing/create_mapfile.py +++ b/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, diff --git a/compass/landice/tests/ismip7_forcing/fracture/remap_utils.py b/compass/landice/ismip7/remap.py similarity index 90% rename from compass/landice/tests/ismip7_forcing/fracture/remap_utils.py rename to compass/landice/ismip7/remap.py index 8da5de4ae4..7a9eb36748 100644 --- a/compass/landice/tests/ismip7_forcing/fracture/remap_utils.py +++ b/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,6 +30,11 @@ 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] @@ -36,7 +42,8 @@ def extrapolate_source(input_file, output_file, varnames, logger): 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] diff --git a/compass/landice/tests/ismip7_calibration/__init__.py b/compass/landice/tests/ismip7_calibration/__init__.py new file mode 100644 index 0000000000..e356776695 --- /dev/null +++ b/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)) diff --git a/compass/landice/tests/ismip7_calibration/ais/__init__.py b/compass/landice/tests/ismip7_calibration/ais/__init__.py new file mode 100644 index 0000000000..ec6ea040a4 --- /dev/null +++ b/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. + + ``_`` + 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') diff --git a/compass/landice/tests/ismip7_calibration/ais/aggregate.py b/compass/landice/tests/ismip7_calibration/ais/aggregate.py new file mode 100644 index 0000000000..2cbbaf71b4 --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/ais/aggregate.py @@ -0,0 +1,294 @@ +""" +Aggregate the MALI melt ensemble into the four objective-function terms. +""" + +import numpy as np +import xarray as xr +from mpas_tools.io import write_netcdf + +from compass.landice.tests.ismip7_calibration.ais import melt_model +from compass.landice.tests.ismip7_calibration.terms import ( + average_by_group, + integrate_by_group, +) +from compass.step import Step + +#: cell dimension of an MPAS mesh +CELL_DIMS = ('nCells',) + + +class Aggregate(Step): + """ + A step that turns the MALI melt fields into unit-parameter aggregates of + the four objective-function terms. + + All aggregation is **area-weighted**. On a 4-20 km variable-resolution + mesh a plain mean is wrong: integrals use ``melt * areaCell`` and means + are weighted by ``areaCell``. + + Because melt is exactly proportional to the melt parameter, each + aggregate is divided by the reference parameter the ensemble was run at, + giving a unit aggregate that the calibration scales onto the whole + parameter grid. That scaling is exact, not an approximation. + + This step also reports MALI's per-basin ice-shelf area against the + observed ISMIP7 extent. J1, J2 and J4 are integrals, so any shelf-area + mismatch enters the calibrated parameter directly; the protocol intends + that, but it should be reported alongside the result. + + Attributes + ---------- + melt_forms : list of str + The melt forms to aggregate + + states : list of compass.landice.tests.ismip7_calibration.datasets.OceanState + The ocean states in the ensemble + """ # noqa: E501 + + def __init__(self, test_case, melt_forms, states): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_calibration.ais.Ais + The test case this step belongs to + + melt_forms : list of str + The melt forms to aggregate + + states : list of datasets.OceanState + The ocean states in the ensemble + """ + super().__init__(test_case=test_case, name='aggregate') + self.melt_forms = melt_forms + self.states = states + + self.add_input_file(filename='masks.nc', + target='../remap_masks/ismip7_masks_on_mali.nc') + for melt_form in melt_forms: + for state in states: + self.add_input_file( + filename=f'melt_{melt_form}_{state.name}.nc', + target=f'../{melt_form}_{state.name}/output_melt.nc') + self.add_output_file(filename=f'aggregates_{melt_form}.nc') + self.add_output_file(filename='shelf_area.nc') + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config['ismip7_calibration'] + base_path_mali = section.get('base_path_mali') + mali_mesh_file = section.get('mali_mesh_file') + self.add_input_file( + filename='mesh.nc', + target=f'{base_path_mali}/{mali_mesh_file}') + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + section = config['ismip7_calibration_melt'] + reference = {'ismip7': section.getfloat('reference_k'), + 'ismip6': section.getfloat('reference_gamma0')} + + static = _load_static(logger) + + for melt_form in self.melt_forms: + logger.info(f'Aggregating the {melt_form} ensemble') + ds = _aggregate_form(self.states, melt_form, + reference[melt_form], static, logger) + write_netcdf(ds, f'aggregates_{melt_form}.nc') + + _report_shelf_area(static, logger) + + +def _load_static(logger): + """Load the mesh, masks and MALI's own floating-cell mask.""" + ds_mesh = xr.open_dataset('mesh.nc') + ds_masks = xr.open_dataset('masks.nc') + + area = ds_mesh['areaCell'] + + basins = ds_masks['ismip7BasinNumber'] + bfrn = ds_masks['ismip7BFRNBin'] + region_code = ds_masks['ismip7ShelfRegion'] + observed_floating = ds_masks['ismip7FloatingMask'] == 1 + + # map the region codes onto the labels the targets use + region_label = xr.where(region_code == 1, 'pig', + xr.where(region_code == 2, 'dotson', '')) + + logger.info(f' {int((basins >= 0).sum())} cells carry a basin number') + return dict(area=area, basins=basins.where(basins >= 0), + bfrn=bfrn.where(bfrn >= 0), region_label=region_label, + observed_floating=observed_floating) + + +def _melt_from_run(filename, reference): + """ + Read one melt field, in kg m-2 yr-1 per unit parameter. + + Dividing by the reference parameter the run used gives melt at a unit + parameter, which is exact because melt is proportional to the parameter. + + The contributing cells are the ones MALI actually computed melt for -- + floating *and* connected to the open ocean. Including the rest would + not change an integral, since their melt is zero, but it would dilute + the area-weighted basin means that J3 is built from. + """ + fields = melt_model.read_run(filename, 'mesh.nc') + return fields['melt'] / reference, fields['floating'] + + +def _aggregate_form(states, melt_form, reference, static, logger): + """Build the four unit aggregates for one melt form.""" + area = static['area'] + basins, bfrn = static['basins'], static['bfrn'] + + by_name = {state.name: state for state in states} + melt = {} + floating = {} + for state in states: + filename = f'melt_{melt_form}_{state.name}.nc' + melt[state.name], floating[state.name] = \ + _melt_from_run(filename, reference) + + # J1 and J2 come from the present-day climatology + pd_melt = melt['climatology'] + pd_floating = floating['climatology'] + t1 = integrate_by_group(pd_melt, area, pd_floating, basins, CELL_DIMS, + group_dim='basins') + t2 = integrate_by_group(pd_melt, area, pd_floating, bfrn, CELL_DIMS, + group_dim='BFRN_bins') + + # J3: the warm-minus-cold basin-mean difference, per ocean model + means = {} + for which in ('cold', 'warm'): + per_model = [] + for state in states: + if state.kind != 'model' or state.state != which: + continue + values = melt[state.name] + # regional models only constrain the basins they cover + if state.basins is not None: + values = values.where(basins.isin(list(state.basins))) + agg = average_by_group(values, area, floating[state.name], + basins, CELL_DIMS, group_dim='basins') + per_model.append(agg.expand_dims({'model': [state.label]})) + means[which] = xr.concat(per_model, dim='model', coords='minimal') + t3 = means['warm'] - means['cold'] + + # J4: PIG and Dotson integrated melt, per observation year + per_year = [] + for state in states: + if state.kind != 'obs': + continue + agg = integrate_by_group(melt[state.name], area, + floating[state.name], + static['region_label'], CELL_DIMS, + group_dim='region') + per_year.append(agg.expand_dims({'year': [state.year]})) + t4 = xr.concat(per_year, dim='year', coords='minimal') + # drop the '' label used for cells outside PIG and Dotson + t4 = t4.sel(region=[region for region in t4.region.values if region]) + + logger.info(f' total present-day melt at the reference parameter: ' + f'{float(t1.sum()) * reference:.1f} Gt/yr') + + ds = xr.Dataset({'t1': t1, 't2': t2, 't3': t3, 't4': t4}) + ds.attrs['melt_form'] = melt_form + ds.attrs['reference_parameter'] = reference + ds.attrs['note'] = ( + 'Aggregates at a unit melt parameter. Melt is exactly proportional ' + 'to the parameter, so the whole parameter ensemble is these times ' + 'each parameter value.') + ds.attrs['ocean_states'] = ', '.join(sorted(by_name)) + return ds + + +def _report_shelf_area(static, logger): + """Report MALI's per-basin shelf area against the observed extent.""" + area = static['area'] + basins = static['basins'] + observed = static['observed_floating'] + + modelled_file = None + for name in ('melt_ismip7_climatology.nc', 'melt_ismip6_climatology.nc'): + try: + modelled = melt_model.read_run(name, 'mesh.nc')['floating'] + modelled_file = name + break + except FileNotFoundError: + continue + if modelled_file is None: + logger.warning('No melt output found; skipping the shelf-area report') + return + + grouper = basins.rename('basins') + modelled_area = (area.where(modelled).groupby(grouper).sum() / 1.0e9) + observed_area = (area.where(observed).groupby(grouper).sum() / 1.0e9) + + logger.info('') + logger.info('Ice-shelf area by ISMIP7 basin, 10^3 km^2:') + logger.info(f' {"basin":>5s} {"MALI":>10s} {"observed":>10s} ' + f'{"ratio":>8s}') + for basin in modelled_area.basins.values: + mali = float(modelled_area.sel(basins=basin)) + obs = float(observed_area.sel(basins=basin)) + ratio = mali / obs if obs > 0 else float('nan') + logger.info(f' {int(basin):5d} {mali:10.1f} {obs:10.1f} ' + f'{ratio:8.2f}') + total_mali = float(modelled_area.sum()) + total_obs = float(observed_area.sum()) + logger.info(f' {"total":>5s} {total_mali:10.1f} {total_obs:10.1f} ' + f'{total_mali / total_obs:8.2f}') + logger.info(' (J1, J2 and J4 are integrals, so a shelf-area mismatch ' + 'enters the calibrated parameter directly; J3 is a mean and ' + 'is much less sensitive)') + logger.info('') + + ds = xr.Dataset({'modelled_shelf_area': modelled_area, + 'observed_shelf_area': observed_area}) + ds['modelled_shelf_area'].attrs['units'] = '10^3 km^2' + ds['observed_shelf_area'].attrs['units'] = '10^3 km^2' + write_netcdf(ds, 'shelf_area.nc') + + +def unit_aggregates(filename): + """ + Read the unit aggregates written by this step. + + Parameters + ---------- + filename : str + An ``aggregates_.nc`` file + + Returns + ------- + units : dict + ``t1``, ``t2``, ``t3`` and ``t4`` + """ + ds = xr.load_dataset(filename) + return {name: ds[name] for name in ('t1', 't2', 't3', 't4')} + + +def basin_coordinate(units): + """ + The ISMIP7 basin numbers present in an aggregate, as integers. + + Parameters + ---------- + units : dict + From :py:func:`unit_aggregates` + + Returns + ------- + basins : numpy.ndarray + The basin numbers + """ + return np.asarray(units['t1'].basins.values, dtype=int) diff --git a/compass/landice/tests/ismip7_calibration/ais/calibrate.py b/compass/landice/tests/ismip7_calibration/ais/calibrate.py new file mode 100644 index 0000000000..be799c5d05 --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/ais/calibrate.py @@ -0,0 +1,163 @@ +""" +Select the melt parameter from the MALI ensemble, per the ISMIP7 protocol. +""" + +import numpy as np +import xarray as xr +from mpas_tools.io import write_netcdf + +from compass.landice.tests.ismip7_calibration import datasets +from compass.landice.tests.ismip7_calibration.ais.aggregate import ( + unit_aggregates, +) +from compass.landice.tests.ismip7_calibration.configure import ( + objective_options, + parameter_name, + parameter_values, + weighting, +) +from compass.landice.tests.ismip7_calibration.objective import ( + build_toolbox_terms, + run_optimisation, +) +from compass.step import Step + +#: above this fraction of draws sitting on the largest parameter in the grid, +#: the grid has clipped the distribution badly enough to matter +MAX_AT_UPPER_BOUND = 0.005 + + +class Calibrate(Step): + """ + A step that runs the ISMIP7 parameter selection on the MALI ensemble. + + The protocol draws random term weights and random targets within their + uncertainties, and for each draw picks the parameter that minimises + ``I = sum_i a_i J_i / median(J_i)``. The 5th, 50th and 95th percentiles + of the resulting distribution are what the ISMIP7 projections need. + + Because the objective **normalises each term by its own median** over the + parameter ensemble, a minimised ``I`` is not comparable between melt + forms -- only between parameter values within one form. This step + therefore reports a distribution per melt form and never a cross-form + comparison of the objective. + + Attributes + ---------- + melt_forms : list of str + The melt forms to calibrate + """ + + def __init__(self, test_case, melt_forms): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_calibration.ais.Ais + The test case this step belongs to + + melt_forms : list of str + The melt forms to calibrate + """ + super().__init__(test_case=test_case, name='calibrate') + self.melt_forms = melt_forms + for melt_form in melt_forms: + self.add_input_file( + filename=f'aggregates_{melt_form}.nc', + target=f'../aggregate/aggregates_{melt_form}.nc') + self.add_output_file(filename=f'calibration_{melt_form}.nc') + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + + base_path = config.get('ismip7_calibration', 'base_path_ismip7') + targets = datasets.load_targets(base_path) + t3_models, t4_regions, t4_years = weighting(config) + sample_size, seed = objective_options(config) + + logger.info(f'J3 weighting: ' + f'{"all models" if t3_models is None else t3_models}') + logger.info(f'J4 weighting: ' + f'{"PIG and Dotson" if t4_regions is None else t4_regions}' + f', ' + f'{"all years" if t4_years is None else t4_years}') + + for melt_form in self.melt_forms: + name = parameter_name(melt_form) + values = parameter_values(config, melt_form) + units = unit_aggregates(f'aggregates_{melt_form}.nc') + + terms = build_toolbox_terms(units, targets, values, + t3_models=t3_models, + t4_regions=t4_regions, + t4_years=t4_years) + + logger.info(f'Sampling the objective function {sample_size} ' + f'times for the {melt_form} form...') + result = run_optimisation(terms, values, sample_size=sample_size, + seed=seed) + + _write(result, values, melt_form, name, + f'calibration_{melt_form}.nc') + _report(result, melt_form, name, logger) + _check_grid_brackets_the_distribution(result, values, melt_form, + name, logger) + + +def _check_grid_brackets_the_distribution(result, values, melt_form, name, + logger): + """ + Check that the parameter grid is wide enough to hold the distribution. + + A draw that chose the largest parameter on the grid wanted a larger one + and could not have it. Enough of those and the distribution piles up + against the end of the grid, and its upper percentiles say more about + where the grid stops than about the objective. + + Raises + ------ + ValueError + If too many draws sit on the largest parameter in the grid + """ + fraction = result['at_upper_bound'] + logger.info(f' draws at the top of the {name} grid ' + f'({values[-1]:.4g}): {100.0 * fraction:.3f}%') + if fraction > MAX_AT_UPPER_BOUND: + raise ValueError( + f'{100.0 * fraction:.2f}% of the draws chose the largest ' + f'{name} on the grid, {values[-1]:.4g}, so the {melt_form} ' + f'distribution is clipped and its upper percentiles are not ' + f'meaningful. Raise the maximum of the {name} grid in the ' + f'[ismip7_calibration_melt] config section.') + + +def _write(result, values, melt_form, name, filename): + """Write one parameter distribution to a file.""" + ds = xr.Dataset() + ds['min_p1'] = ('sample', result['min_p1']) + ds['parameter_values'] = ('parameter', np.asarray(values)) + for key in ('p5', 'median', 'p95', 'mode'): + ds[key] = float(result[key]) + ds.attrs['melt_form'] = melt_form + ds.attrs['parameter'] = name + ds.attrs['note'] = ( + 'The objective normalises each term by its own median over the ' + 'parameter ensemble, so the minimised objective is not comparable ' + 'between melt forms -- only between parameter values within one ' + 'form.') + write_netcdf(ds, filename) + + +def _report(result, melt_form, name, logger): + """Log the selected percentiles.""" + logger.info('') + logger.info(f'{melt_form}: {name} percentiles on the MALI mesh') + for key in ('p5', 'median', 'p95'): + logger.info(f' {key:8s} {result[key]:12.5e}') + logger.info(f' {"mode":8s} {result["mode"]:12.5e}') + logger.info('') diff --git a/compass/landice/tests/ismip7_calibration/ais/fit_delta_t.py b/compass/landice/tests/ismip7_calibration/ais/fit_delta_t.py new file mode 100644 index 0000000000..6d642eb8a9 --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/ais/fit_delta_t.py @@ -0,0 +1,223 @@ +""" +Fit the per-basin thermal-forcing correction dT_b, after parameter selection. +""" + +import os + +import numpy as np +import pandas as pd +import xarray as xr +from mpas_tools.io import write_netcdf + +from compass.landice.tests.ismip7_calibration.ais import melt_model +from compass.landice.tests.ismip7_calibration.configure import parameter_name +from compass.step import Step + + +class FitDeltaT(Step): + """ + A step that fits the basin thermal-forcing correction ``dT_b``. + + Protocol Sect. 4.2.1 permits the corrections to be calculated at either + stage: **(1)** before the parameter optimisation, which "can lead to the + problem that the parameter bounds are not constrained when present-day + melt rates are compared to observations (namely, terms J1 and J2)"; or + **(2)** after it, "as in ISMIP6". This step follows **(2)**, which is + what the published quadratic worked example does, and it is why the + calibration itself runs with ``dT_b = 0`` throughout. + + For each basin, ``dT`` is searched on a grid within the configured + bounds, minimising the absolute difference between the basin-integrated + melt and the IMBIE observed value. The protocol recommends bounding + ``dT_b`` by about 2 K, since the warmest-to-coldest spread across + Antarctic ice shelves is only 3-4 K, and recommends avoiding the + correction where possible since it is "only an ad-hoc correction". + + The search is done on the MALI mesh with **area-weighted** integrals. + The toolbox's own ``optimise_deltaT`` cannot be used, because it sums + over structured-grid ``x`` and ``y`` dimensions with a single scalar cell + area. + + Melt is recomputed in Python from MALI's ``TFdraft``, so no further MALI + runs are needed: ``dT`` enters the melt expression only through the + thermal forcing. + + Attributes + ---------- + melt_forms : list of str + The melt forms to fit dT_b for + """ + + def __init__(self, test_case, melt_forms, state_name='climatology'): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_calibration.ais.Ais + The test case this step belongs to + + melt_forms : list of str + The melt forms to fit dT_b for + + state_name : str, optional + The ocean state to fit against; the present-day climatology, to + match the observed present-day melt + """ + super().__init__(test_case=test_case, name='fit_delta_t') + self.melt_forms = melt_forms + self.state_name = state_name + + self.add_input_file(filename='masks.nc', + target='../remap_masks/ismip7_masks_on_mali.nc') + for melt_form in melt_forms: + self.add_input_file( + filename=f'melt_{melt_form}.nc', + target=f'../{melt_form}_{state_name}/output_melt.nc') + self.add_input_file( + filename=f'calibration_{melt_form}.nc', + target=f'../calibrate/calibration_{melt_form}.nc') + self.add_output_file(filename=f'melt_params_{melt_form}.nc') + + def setup(self): + """ + Set up this step of the test case + """ + section = self.config['ismip7_calibration'] + self.add_input_file( + filename='mesh.nc', + target=(f'{section.get("base_path_mali")}/' + f'{section.get("mali_mesh_file")}')) + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + + if not config.getboolean('ismip7_calibration_delta_t', + 'fit_delta_t'): + logger.info('fit_delta_t is False; not fitting dT_b') + return + + base_path = config.get('ismip7_calibration', 'base_path_ismip7') + observed = _load_observed(base_path) + + section = config['ismip7_calibration_delta_t'] + delta_t_grid = np.linspace(section.getfloat('delta_t_min'), + section.getfloat('delta_t_max'), + section.getint('delta_t_count')) + + ds_masks = xr.open_dataset('masks.nc') + + for melt_form in self.melt_forms: + with xr.open_dataset(f'calibration_{melt_form}.nc') as ds_calib: + parameter = float(ds_calib['median']) + name = parameter_name(melt_form) + logger.info(f'Fitting dT_b for the {melt_form} form at the ' + f'median {name} = {parameter:.5e}') + + fields = melt_model.read_run(f'melt_{melt_form}.nc', + 'mesh.nc') + result = _fit(melt_form, parameter, fields, ds_masks, observed, + delta_t_grid, config, logger) + _write_params(result, ds_masks, melt_form, parameter, name, + f'melt_params_{melt_form}.nc') + + +def _load_observed(base_path): + """The IMBIE basin-integrated melt targets, in Gt yr-1.""" + path = os.path.join(base_path, 'parameterisations', 'ocean', 'meltobs', + 'Melt_Paolo_Davison_Adusumilli_imbie2.csv') + return pd.read_csv(path, index_col=0)['BMR (Gt/yr)'] + + +def _fit(melt_form, parameter, fields, ds_masks, observed, delta_t_grid, + config, logger): + """Grid-search dT per basin against the observed integrated melt.""" + area = fields['area'] + floating = fields['floating'] + # the ISMIP7 0-based numbering, which is what the observations use + basin = ds_masks['ismip7BasinNumber'] + + # the basin mean of the ISMIP6 form shifts with dT exactly as the local + # forcing does, because dT is constant within a basin + mean_tf = None + if melt_form == 'ismip6': + mean_tf = melt_model.basin_mean_tf(fields['tf_draft'], area, + floating, fields['basin']) + + basins = sorted(int(value) for value in np.unique(basin.values) + if value >= 0) + + optimal = {} + residual = {} + modelled_zero = {} + for basin_number in basins: + if basin_number not in observed.index: + continue + target = float(observed.loc[basin_number]) + in_basin = (basin == basin_number) & floating + if not bool(in_basin.any()): + continue + + totals = [] + for delta_t in delta_t_grid: + melt = melt_model.melt_from_tf(melt_form, parameter, fields, + config, delta_t=delta_t, + mean_tf=mean_tf) + totals.append(float((melt * area).where(in_basin).sum()) / 1.0e12) + totals = np.asarray(totals) + + best = int(np.argmin(np.abs(totals - target))) + optimal[basin_number] = float(delta_t_grid[best]) + residual[basin_number] = float(abs(totals[best] - target)) + zero = int(np.argmin(np.abs(delta_t_grid))) + modelled_zero[basin_number] = float(totals[zero]) + + logger.info('') + logger.info(f'dT_b per ISMIP7 basin, {melt_form}:') + logger.info(f' {"basin":>5s} {"dT_b (K)":>9s} {"melt(0)":>10s} ' + f'{"observed":>10s} {"residual":>10s}') + for basin_number in sorted(optimal): + logger.info(f' {basin_number:5d} {optimal[basin_number]:9.3f} ' + f'{modelled_zero[basin_number]:10.1f} ' + f'{float(observed.loc[basin_number]):10.1f} ' + f'{residual[basin_number]:10.2f}') + at_bound = [basin_number for basin_number, value in optimal.items() + if abs(value) >= 0.999 * max(abs(delta_t_grid[0]), + abs(delta_t_grid[-1]))] + if at_bound: + logger.warning(f' dT_b hit the bounds in basins {at_bound}; the ' + f'protocol bounds the correction deliberately, so a ' + f'basin at the bound means the melt form cannot match ' + f'the observed integral there.') + logger.info('') + return optimal + + +def _write_params(optimal, ds_masks, melt_form, parameter, name, filename): + """Write a MALI melt-parameter file carrying the fitted dT_b.""" + basin0 = ds_masks['ismip7BasinNumber'] + delta_t = xr.zeros_like(basin0, dtype=float) + for basin_number, value in optimal.items(): + delta_t = xr.where(basin0 == basin_number, value, delta_t) + + ds = xr.Dataset() + ds['ismip6shelfMelt_basin'] = ds_masks['ismip6shelfMelt_basin'] + ds['ismip6shelfMelt_deltaT'] = delta_t + ds['ismip6shelfMelt_deltaT'].attrs = { + 'long_name': 'basin-wide thermal forcing correction', + 'units': 'degC'} + if melt_form == 'ismip6': + ds['ismip6shelfMelt_gamma0'] = parameter + ds['ismip6shelfMelt_gamma0'].attrs = {'units': 'm yr^-1'} + + ds.attrs['melt_form'] = melt_form + ds.attrs['parameter_name'] = name + ds.attrs['parameter_value'] = parameter + ds.attrs['note'] = ( + 'dT_b fitted after parameter selection, per protocol Sect. 4.2.1 ' + 'option 2, against the IMBIE basin-integrated melt observations.') + write_netcdf(ds, filename) diff --git a/compass/landice/tests/ismip7_calibration/ais/make_graph.py b/compass/landice/tests/ismip7_calibration/ais/make_graph.py new file mode 100644 index 0000000000..2680534026 --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/ais/make_graph.py @@ -0,0 +1,58 @@ +""" +Build the graph partition the MALI melt diagnostics run on. +""" + +import os + +from compass.model import make_graph_file, partition +from compass.step import Step + + +class MakeGraph(Step): + """ + A step that builds a graph partition file for the MALI mesh. + + The melt-diagnostic runs all use the same mesh and the same number of + tasks, so the partition is built once here rather than once per run. + Doing it in the test group means the group works on any MALI mesh, + rather than only on meshes that already have a partition file + distributed alongside them. + """ + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_calibration.ais.Ais + The test case this step belongs to + """ + super().__init__(test_case=test_case, name='make_graph') + + def setup(self): + """ + Set up this step of the test case + """ + section = self.config['ismip7_calibration'] + ntasks = section.getint('ntasks') + self.add_input_file( + filename='mesh.nc', + target=(f'{section.get("base_path_mali")}/' + f'{section.get("mali_mesh_file")}')) + self.add_output_file(filename=f'graph.info.part.{ntasks}') + + def run(self): + """ + Run this step of the test case + """ + config = self.config + ntasks = config.getint('ismip7_calibration', 'ntasks') + + if not os.path.exists('graph.info'): + self.logger.info('Building graph.info from the MALI mesh') + make_graph_file(mesh_filename='mesh.nc', + graph_filename='graph.info') + + self.logger.info(f'Partitioning the graph for {ntasks} tasks') + partition(ntasks, config, self.logger, graph_file='graph.info') diff --git a/compass/landice/tests/ismip7_calibration/ais/melt_model.py b/compass/landice/tests/ismip7_calibration/ais/melt_model.py new file mode 100644 index 0000000000..ff8d09f883 --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/ais/melt_model.py @@ -0,0 +1,323 @@ +""" +Recompute MALI's melt from its diagnostic output. + +Both melt forms are functions of the thermal forcing at the ice draft, which +MALI writes as ``ismip6shelfMelt_TFdraft``. That means melt can be +recomputed in Python for any melt parameter and any basin correction +``dT_b`` without re-running MALI, which is what the ``dT_b`` fit needs and +what the verification step compares against. +""" + +import dataclasses + +import numpy as np +import xarray as xr + +from compass.landice.tests.ismip7_calibration.quadratic import ( + MALI, + angle_from_sin_slope, + local_quadratic_melt, + nonlocal_quadratic_melt, +) + +#: bit of ``cellMask`` marking floating ice +FLOATING_MASK_BIT = 4 + +#: seconds in a year, matching MALI's ``scyr`` and its noleap calendar +SECONDS_PER_YEAR = 31536000.0 + +#: rate of change of the freezing temperature of seawater with depth, +#: K m-1, matching MALI's ``oceanFreezingTempDepthDependence`` +FREEZING_TEMP_DEPTH_DEPENDENCE = -7.53e-4 + + +def read_run(filename, mesh_file): + """ + Read the fields a melt diagnostic run writes. + + Parameters + ---------- + filename : str + An ``output_melt.nc`` written by a + :py:class:`~compass.landice.tests.ismip7_calibration.ais.run_state.RunState` + step + + mesh_file : str + The MALI mesh file the run started from, used for the geometry the + melt was computed from + + Returns + ------- + fields : dict + ``melt`` (kg m-2 yr-1, positive for melting), ``tf_draft``, + ``floating``, ``basin`` (MALI's 1-based numbering), ``delta_t`` and + ``area`` + + Notes + ----- + **The contributing cells come from the initial geometry, not from the + output.** MALI computes melt where the ice is floating *and* the cell + is connected to the open ocean, evaluated on the geometry at the start + of the timestep. The ``cellMask`` in the output is *post*-step, and + over a single step some cells with no ice at the start end up with a + trace of it, so the output mask marks thousands of cells as floating + that MALI never computed melt for. Their melt is zero, which does not + change an integral, but it would dilute the area-weighted basin means + that J3 is built from -- and it would make any comparison against an + independent implementation disagree on cells MALI never calculated. + """ # noqa: E501 + thickness, bed = _initial_geometry(mesh_file) + ice = thickness > 0.0 + # MALI counts ice exactly at flotation as floating + floating = ice & (MALI.rho_ice / MALI.rho_ocean * thickness <= -bed) + + with xr.open_dataset(filename) as ds: + bmb = ds['floatingBasalMassBal'].isel(Time=0) + connected = ds['connectedOceanMask'].isel(Time=0) == 1 + fields = dict( + melt=(-bmb * SECONDS_PER_YEAR).compute(), + tf_draft=ds['ismip6shelfMelt_TFdraft'].isel(Time=0).compute(), + floating=(floating & connected).compute(), + basin=ds['ismip6shelfMelt_basin'].compute(), + delta_t=ds['ismip6shelfMelt_deltaT'].compute(), + area=ds['areaCell'].compute()) + return fields + + +def _initial_geometry(mesh_file): + """Ice thickness and bed topography at the start of the timestep.""" + with xr.open_dataset(mesh_file) as ds: + thickness = ds['thickness'] + bed = ds['bedTopography'] + if 'Time' in thickness.dims: + thickness = thickness.isel(Time=0) + if 'Time' in bed.dims: + bed = bed.isel(Time=0) + return thickness.compute(), bed.compute() + + +def initial_draft(mesh_file, constants=None): + """ + Reconstruct the ice draft from the mesh file, before the timestep. + + **The diagnostic run evolves the geometry.** Melt is applied over the + single timestep and thins the ice by up to a metre or so, so the + ``lowerSurface`` and ``thickness`` written to the output are *post*-step + while ``TFdraft`` was computed *pre*-step. Any check that pairs output + geometry with output melt is therefore inconsistent by up to a metre of + draft; the initial draft has to come from the mesh file instead. + + Parameters + ---------- + mesh_file : str + The MALI mesh file the run started from + + constants : compass.landice.tests.ismip7_calibration.quadratic.Constants, optional + The densities to use; defaults to MALI's + + Returns + ------- + draft : xarray.DataArray + The ice draft, negative below sea level + """ # noqa: E501 + if constants is None: + constants = MALI + thickness, bed = _initial_geometry(mesh_file) + + floating_draft = -constants.rho_ice / constants.rho_ocean * thickness + # where the ice is grounded the draft is the bed + return xr.where(floating_draft > bed, floating_draft, bed) + + +def basin_mean_tf(tf_draft, area, floating, basin): + """ + Area-weighted mean thermal forcing per basin, mapped back onto cells. + + This is what MALI's ISMIP6 non-local form computes internally, and it + must be area-weighted: on a 4-20 km mesh a plain mean is wrong. + + Parameters + ---------- + tf_draft : xarray.DataArray + Thermal forcing at the ice draft + + area : xarray.DataArray + Cell area + + floating : xarray.DataArray + True on floating cells + + basin : xarray.DataArray + Basin number per cell + + Returns + ------- + mean_tf : xarray.DataArray + The basin mean, broadcast back onto every cell + """ + weights = area.where(floating) + numer = (tf_draft * weights).groupby(basin.rename('basin')).sum() + denom = weights.groupby(basin.rename('basin')).sum() + means = numer / denom.where(denom != 0.0) + return means.sel(basin=basin.rename('basin')).drop_vars('basin') + + +def melt_from_tf(melt_form, parameter, fields, config, delta_t=None, + mean_tf=None): + """ + Recompute the melt field from the thermal forcing at the draft. + + Parameters + ---------- + melt_form : {'ismip7', 'ismip6'} + Which melt form to evaluate + + parameter : float or xarray.DataArray + ``K`` for ``'ismip7'`` or ``gamma0`` for ``'ismip6'`` + + fields : dict + From :py:func:`read_run` + + config : compass.config.CompassConfigParser + Configuration options, for the constant salinity, slope and Coriolis + parameter + + delta_t : float or xarray.DataArray, optional + The basin correction to apply; defaults to the field the run used + + mean_tf : xarray.DataArray, optional + Precomputed basin-mean thermal forcing, for ``'ismip6'`` + + Returns + ------- + melt : xarray.DataArray + Melt rate in kg m-2 yr-1, positive for melting, zero off the shelves + """ + section = config['ismip7_calibration_melt'] + if delta_t is None: + delta_t = fields['delta_t'] + + tf_draft = fields['tf_draft'] + floating = fields['floating'] + + if melt_form == 'ismip7': + # config gives the *sine* of the slope while local_quadratic_melt + # takes the angle + slope = angle_from_sin_slope(section.getfloat('sin_slope')) + constants = dataclasses.replace( + MALI, coriolis=section.getfloat('coriolis')) + melt = local_quadratic_melt( + parameter, tf_draft, section.getfloat('salinity'), slope, + constants=constants, delta_t=delta_t) + elif melt_form == 'ismip6': + if mean_tf is None: + mean_tf = basin_mean_tf(tf_draft, fields['area'], floating, + fields['basin']) + melt = nonlocal_quadratic_melt(parameter, tf_draft, mean_tf, + constants=MALI, delta_t=delta_t) + else: + raise ValueError(f"melt_form must be 'ismip7' or 'ismip6', but is " + f"'{melt_form}'") + + return xr.where(floating, melt, 0.0) + + +def integrate_by_basin(melt, area, floating, basin): + """ + Integrate a melt rate over each basin, in Gt yr-1. + + Parameters + ---------- + melt : xarray.DataArray + Melt rate in kg m-2 yr-1 + + area : xarray.DataArray + Cell area in m2 + + floating : xarray.DataArray + True on floating cells + + basin : xarray.DataArray + Basin number per cell + + Returns + ------- + total : xarray.DataArray + Integrated melt per basin, in Gt yr-1 + """ + weighted = (melt * area).where(floating) + return weighted.groupby(basin.rename('basin')).sum() / 1.0e12 + + +def interpolate_to_draft(field_3d, z_ocean, draft, bed, + freezing_correction=False): + """ + Interpolate a 3-D ocean field to the ice draft. + + Written from the protocol rather than transliterated from MALI's + Fortran, so that comparing the two is a real check of the + implementation. The four cases match MALI's: above the shallowest layer + centre, below the deepest, where the layer below the draft is beneath + the bed, and linear interpolation between layer centres. + + Parameters + ---------- + field_3d : numpy.ndarray + The field, shaped ``(nCells, nLayers)`` + + z_ocean : numpy.ndarray + Layer centre depths, negative downward and decreasing + + draft : numpy.ndarray + Ice draft per cell, negative below sea level + + bed : numpy.ndarray + Bed topography per cell + + freezing_correction : bool, optional + Whether to apply the depth dependence of the freezing temperature + when the draft is below the layer centre being used. MALI applies + it to the thermal forcing in the two branches that extrapolate + downward -- below the deepest layer centre, and where the layer + below the draft is beneath the bed -- so pass True when comparing + against ``ismip6shelfMelt_TFdraft``. + + Returns + ------- + at_draft : numpy.ndarray + The field interpolated to the draft, per cell + """ + n_cells = field_3d.shape[0] + at_draft = np.full(n_cells, np.nan) + n_layers = len(z_ocean) + rate = FREEZING_TEMP_DEPTH_DEPENDENCE if freezing_correction else 0.0 + + for index in range(n_cells): + # ksup is the deepest layer centre still at or above the draft + above = np.nonzero(z_ocean >= draft[index])[0] + ksup = above[-1] if above.size > 0 else -1 + if ksup < 0: + # above the shallowest centre: take the shallowest layer, with + # no correction, since the draft is above it + at_draft[index] = field_3d[index, 0] + elif ksup == n_layers - 1: + # below the deepest centre: take the deepest layer, corrected + # for the freezing point at the draft + at_draft[index] = ( + field_3d[index, n_layers - 1] - + (z_ocean[n_layers - 1] - draft[index]) * rate) + elif z_ocean[ksup + 1] < bed[index]: + # the layer below the draft is beneath the bed, so there is no + # water there to interpolate into; take the layer above, + # corrected for the freezing point at the draft + at_draft[index] = ( + field_3d[index, ksup] - + (z_ocean[ksup] - draft[index]) * rate) + else: + # between layer centres, interpolate linearly in depth + span = z_ocean[ksup] - z_ocean[ksup + 1] + w_deep = (z_ocean[ksup] - draft[index]) / span + w_shallow = (draft[index] - z_ocean[ksup + 1]) / span + at_draft[index] = (w_deep * field_3d[index, ksup + 1] + + w_shallow * field_3d[index, ksup]) + return at_draft diff --git a/compass/landice/tests/ismip7_calibration/ais/namelist.landice b/compass/landice/tests/ismip7_calibration/ais/namelist.landice new file mode 100644 index 0000000000..8a904f8e2d --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/ais/namelist.landice @@ -0,0 +1,11 @@ +config_velocity_solver = 'none' +config_thermal_solver = 'none' +config_thermal_calculate_bmb = .false. +config_calving = 'none' +config_front_mass_bal_grounded = 'none' +config_ocean_data_extrapolation = .false. +config_do_restart = .false. +config_start_time = '0000-01-01_00:00:00' +config_stop_time = 'none' +config_adaptive_timestep = .false. +config_write_output_on_startup = .false. diff --git a/compass/landice/tests/ismip7_calibration/ais/remap_forcing.py b/compass/landice/tests/ismip7_calibration/ais/remap_forcing.py new file mode 100644 index 0000000000..d26d6d7e89 --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/ais/remap_forcing.py @@ -0,0 +1,251 @@ +""" +Remap the ISMIP7 calibration thermal forcing onto a MALI mesh. +""" + +import os + +import numpy as np +import xarray as xr +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call + +from compass.landice.ismip7.mapping import build_mapping_file +from compass.landice.ismip7.remap import extrapolate_source +from compass.landice.tests.ismip7_calibration import datasets +from compass.step import Step + + +class RemapForcing(Step): + """ + A step that remaps the ISMIP7 calibration thermal forcing onto the MALI + mesh, one file per ocean state. + + This mirrors ``ismip7_forcing/ocean_thermal``, driven by the calibration + ocean states rather than by CMIP scenario files. + + **Only thermal forcing is remapped.** The ISMIP7 ocean forcing processed + for the MALI projections carries thermal forcing and no salinity, so a + projection has no salinity field to read. The melt parameterization is + run with a constant salinity instead + (``config_ismip7_melt_salinity_source = 'constant'``); calibrating + against a spatially varying salinity would tune the melt parameter for + physics the projections cannot run. + + Attributes + ---------- + states : list of compass.landice.tests.ismip7_calibration.datasets.OceanState + The ocean states to remap + """ # noqa: E501 + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_calibration.ais.Ais + The test case this step belongs to + """ + super().__init__(test_case=test_case, name='remap_forcing') + self.states = [] + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config['ismip7_calibration'] + base_path = section.get('base_path_ismip7') + base_path_mali = section.get('base_path_mali') + mali_mesh_file = section.get('mali_mesh_file') + subset = section.get('ocean_state_subset') + + self.add_input_file( + filename=mali_mesh_file, + target=os.path.join(base_path_mali, mali_mesh_file)) + + self.states = datasets.ocean_states(base_path, subset=subset) + missing = datasets.missing_files(self.states) + # only the thermal forcing is needed; a missing salinity file is not + # a problem, since a constant salinity is used + missing = [(name, path) for name, path in missing + if not path.endswith('_S.nc') and '/so/' not in path] + if missing: + listing = '\n '.join(f'{name}: {path}' for name, path in missing) + raise FileNotFoundError( + f'{len(missing)} ISMIP7 thermal-forcing files are ' + f'missing:\n {listing}') + + for state in self.states: + self.add_output_file(filename=f'forcing_{state.name}.nc') + + self.ntasks = section.getint('esmf_ntasks') + self.min_tasks = self.ntasks + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + + section = config['ismip7_calibration'] + base_path = section.get('base_path_ismip7') + mali_mesh_file = section.get('mali_mesh_file') + mali_mesh_name = section.get('mali_mesh_name') + ntasks = section.getint('esmf_ntasks') + + method = config.get('ismip7_calibration_forcing', 'method_remap') + + res = datasets.ISMIP_RESOLUTION_KM + mapping_file = (f'map_ismip{res}km_to_{mali_mesh_name}_' + f'{method}.nc') + # any calibration forcing file defines the source grid; they share it + build_mapping_file(config, logger, self.states[0].tf_file, + mapping_file, mali_mesh_file=mali_mesh_file, + method_remap=method, projection='ais-bedmap2', + ntasks=ntasks) + + clim_tf_file, _ = datasets.climatology_files(base_path) + clim_tf = xr.open_dataset(clim_tf_file)['tf'] + + for state in self.states: + output_file = f'forcing_{state.name}.nc' + if os.path.exists(output_file): + logger.info(f' {output_file} exists, skipping') + continue + logger.info(f'Remapping {state.name}') + _remap_state(state, clim_tf, mapping_file, output_file, logger) + + +def _remap_state(state, clim_tf, mapping_file, output_file, logger): + """Fill, extrapolate, remap and rename one ocean state.""" + filled_file = f'filled_{state.name}.nc' + extrap_file = f'extrap_{state.name}.nc' + remapped_file = f'remapped_{state.name}.nc' + + filled = _fill_from_climatology(state, clim_tf, filled_file, logger) + + # a safety net: after the climatology fill there should be nothing left + # to extrapolate, but a stray gap would otherwise be blended into + # neighbouring cells by the bilinear remap + extrapolate_source(filled_file, extrap_file, 'tf', logger) + + check_call(['ncremap', '-i', extrap_file, '-o', remapped_file, + '-m', mapping_file, '-v', 'tf'], logger=logger) + + _to_mali_form(remapped_file, state, output_file, filled) + + for path in (filled_file, extrap_file, remapped_file): + if os.path.exists(path): + os.remove(path) + + +def _fill_from_climatology(state, clim_tf, filled_file, logger): + """ + Fill gaps in a partially covering ocean state from the climatology. + + The ISMIP7 near-ice-shelf observational datasets apply a single profile + to the Amundsen basin only (protocol Sect. A10) and are undefined + elsewhere; about 6% of the ISMIP grid is valid. The regional *model* + datasets, by contrast, are distributed already filled with the ISMIP7 + climatology outside their domains (protocol Sect. A9), so they arrive + complete. + + MALI needs valid forcing everywhere it has ice, so the observational + states are filled the same way the model states already were. This does + not affect the calibration: J4 aggregates only over Pine Island and + Dotson, both inside the covered basin. + + Filling happens *before* remapping, so that interpolation never blends a + real value with a missing one at the edge of the covered region. + + Returns + ------- + fraction : float + Fraction of the source grid that was filled, recorded in the output + """ + ds = xr.open_dataset(state.tf_file, decode_times=False) + valid = ds['tf'].notnull() + fraction = float((~valid).mean()) + + if fraction > 0.0: + logger.info(f' filling {100.0 * fraction:.1f}% of the source grid ' + f'from the climatology') + attrs = ds['tf'].attrs + ds['tf'] = ds['tf'].where(valid, clim_tf) + ds['tf'].attrs = attrs + + if '_FillValue' in ds['tf'].encoding: + del ds['tf'].encoding['_FillValue'] + write_netcdf(ds, filled_file) + ds.close() + return fraction + + +def _to_mali_form(remapped_file, state, output_file, filled): + """ + Put the remapped thermal forcing into the names, dimensions and order + MALI expects. + + Two traps are avoided here. MALI takes dimension sizes from its input + stream, so ``nISMIP6OceanLayers`` has to be a real dimension of this + file. And a variable whose name matches a dimension breaks MALI's + reader, so the vertical coordinate is written as + ``ismip6shelfMelt_zOcean`` rather than as a coordinate variable named + after the dimension. + """ + ds = xr.open_dataset(remapped_file, decode_times=False) + + z_ocean = ds['z'].values + + rename_dims = {} + if 'ncol' in ds.dims: + rename_dims['ncol'] = 'nCells' + if 'z' in ds.dims: + rename_dims['z'] = 'nISMIP6OceanLayers' + ds = ds.rename(rename_dims) + ds = ds.rename({'tf': 'ismip6shelfMelt_3dThermalForcing'}) + + tf = ds['ismip6shelfMelt_3dThermalForcing'] + if 'Time' not in tf.dims: + tf = tf.expand_dims('Time', axis=0) + # Registry order is ``nISMIP6OceanLayers nCells Time``, which in C order + # is Time, nCells, nISMIP6OceanLayers + tf = tf.transpose('Time', 'nCells', 'nISMIP6OceanLayers') + + ds_out = xr.Dataset() + ds_out['ismip6shelfMelt_3dThermalForcing'] = tf.astype(float) + ds_out['ismip6shelfMelt_3dThermalForcing'].attrs = { + 'long_name': 'thermal forcing for the ISMIP6/ISMIP7 ice-shelf ' + 'melting methods', + 'units': 'degC'} + ds_out['ismip6shelfMelt_3dThermalForcing'].encoding.clear() + + ds_out['ismip6shelfMelt_zOcean'] = ('nISMIP6OceanLayers', z_ocean) + ds_out['ismip6shelfMelt_zOcean'].attrs = { + 'long_name': 'depth coordinate for the ocean thermal forcing', + 'units': 'm'} + + ds_out['xtime'] = ('Time', ['0000-01-01_00:00:00'.ljust(64)]) + ds_out['xtime'] = ds_out.xtime.astype('S') + + ds_out.attrs['ocean_state'] = state.name + ds_out.attrs['objective_term'] = state.term + ds_out.attrs['source_file'] = state.tf_file + ds_out.attrs['climatology_filled_fraction'] = filled + if state.basins is not None: + ds_out.attrs['constrains_ismip7_basins'] = ', '.join( + str(basin) for basin in state.basins) + + finite = np.isfinite(ds_out['ismip6shelfMelt_3dThermalForcing'].values) + if not finite.all(): + raise ValueError( + f'The remapped thermal forcing for {state.name} is not finite ' + f'everywhere: {int((~finite).sum())} of {finite.size} values are ' + f'missing. MALI would produce invalid melt in those cells.') + + # Time must be UNLIMITED or MALI's reader mishandles the file, so this + # is written with xarray directly rather than through write_netcdf + ds_out.to_netcdf(output_file, unlimited_dims=['Time']) + ds.close() diff --git a/compass/landice/tests/ismip7_calibration/ais/remap_masks.py b/compass/landice/tests/ismip7_calibration/ais/remap_masks.py new file mode 100644 index 0000000000..1adde1d845 --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/ais/remap_masks.py @@ -0,0 +1,405 @@ +""" +Remap the ISMIP7 Antarctic masks onto a MALI mesh. +""" + +import os + +import numpy as np +import xarray as xr +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call + +from compass.landice.ismip7.mapping import build_mapping_file +from compass.landice.tests.ismip7_calibration import datasets +from compass.step import Step + +#: codes written to ``ismip7ShelfRegion`` +REGION_CODES = {'none': 0, 'pig': 1, 'dotson': 2} + +#: below this level of agreement, the basin numbering is probably mismatched +MIN_BASIN_AGREEMENT = 80.0 + +#: the ISMIP7 (0-based) basin Pine Island and Dotson both drain into +EASTERN_AMUNDSEN_BASIN = 9 + +#: the fraction of PIG and Dotson cells that must land in that basin +MIN_SHELF_BASIN_AGREEMENT = 90.0 + + +class RemapMasks(Step): + """ + A step that remaps the ISMIP7 masks onto the MALI mesh. + + The calibration aggregates modelled melt over IMBIE2 drainage basins + (J1, J3), over bins of equal buttressing importance (J2), and over the + Pine Island and Dotson ice shelves (J4). Those masks are distributed on + the ISMIP polar stereographic grid; this step puts them on the MALI mesh + so that the aggregation happens in MALI's own discretization. + + Remapping is nearest neighbour throughout, since every field is + categorical. + + **Two basin-numbering conventions are in play and they differ by one.** + ISMIP7's ``basin_numbers_ismip8km_v2.nc`` is 0-based, 0-15, with basin 9 + the Eastern Amundsen and basin 14 Ronne-Filchner, which is how the + protocol refers to them. MALI's ``ismip6shelfMelt_basin`` is 1-based, + 1-16. Both are written here, under distinct names, so that neither can + be silently reinterpreted as the other -- a mistake that produces + plausible-looking but wrong basin aggregates rather than an error. + """ + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_calibration.ais.Ais + The test case this step belongs to + """ + super().__init__(test_case=test_case, name='remap_masks') + self.add_output_file(filename='ismip7_masks_on_mali.nc') + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config['ismip7_calibration'] + base_path_mali = section.get('base_path_mali') + mali_mesh_file = section.get('mali_mesh_file') + + self.add_input_file( + filename=mali_mesh_file, + target=os.path.join(base_path_mali, mali_mesh_file)) + + region_mask_file = section.get('region_mask_file') + if region_mask_file != 'None': + self.add_input_file( + filename=region_mask_file, + target=os.path.join(base_path_mali, region_mask_file)) + + self.ntasks = section.getint('esmf_ntasks') + self.min_tasks = self.ntasks + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + + section = config['ismip7_calibration'] + base_path = section.get('base_path_ismip7') + mali_mesh_file = section.get('mali_mesh_file') + mali_mesh_name = section.get('mali_mesh_name') + region_mask_file = section.get('region_mask_file') + ntasks = section.getint('esmf_ntasks') + + method = config.get('ismip7_calibration_masks', 'method_remap') + + logger.info(f'Loading the ISMIP7 masks from {base_path}') + ds_masks = _load_ismip7_masks(base_path) + + # any of the mask files defines the source grid; they share it + src_grid_file = datasets.mask_files(base_path)['basins'] + + ds_remapped = _remap_to_mali( + config, ds_masks, src_grid_file, mali_mesh_file, mali_mesh_name, + method, ntasks, logger) + + ds_out = _to_integer_masks(ds_remapped) + + # The ISMIP6 non-local form reads gamma0 from this same input + # stream, and its Registry default is zero, which would silently + # zero the melt field. The ensemble runs at a reference value that + # the aggregation divides out again. + reference_gamma0 = config.getfloat('ismip7_calibration_melt', + 'reference_gamma0') + ds_out['ismip6shelfMelt_gamma0'] = reference_gamma0 + ds_out['ismip6shelfMelt_gamma0'].attrs = { + 'long_name': 'gamma0 for the ISMIP6 ice-shelf melting method', + 'units': 'm yr^-1', + 'note': 'a reference value only; melt is proportional to gamma0 ' + 'and the calibration scales this away'} + + _check_shelves_are_in_their_basin(ds_out, logger) + + if region_mask_file != 'None': + _cross_check_basins(ds_out, region_mask_file, logger) + else: + logger.info('No region mask supplied, so the cross-check against ' + 'the ISMIP6-era regions is skipped.') + + ds_out.attrs['source'] = ( + f'ISMIP7 masks from {base_path} remapped onto {mali_mesh_file}') + ds_out.attrs['remap_method'] = method + + write_netcdf(ds_out, 'ismip7_masks_on_mali.nc') + logger.info('Wrote ismip7_masks_on_mali.nc') + + +def _load_ismip7_masks(base_path): + """Load the ISMIP7 masks on the ISMIP grid, as floats for remapping.""" + files = datasets.mask_files(base_path) + + basins = xr.open_dataset(files['basins'])['basinNumber'] + bfrn = xr.open_dataset(files['bfrn'])['BFRN_bins'] + floating = xr.open_dataset(files['floating'])['mask'] + + shelves = xr.open_dataset(files['shelves'])['shelf_mask'] + if 'time' in shelves.dims: + shelves = shelves.isel(time=0) + + # restrict Pine Island to its main trunk, as the worked example does + pig = (shelves == datasets.PIG_ID) & (shelves['x'] > datasets.PIG_X_MAX) + dotson = shelves == datasets.DOTSON_ID + region = xr.where(pig, REGION_CODES['pig'], + xr.where(dotson, REGION_CODES['dotson'], + REGION_CODES['none'])) + + ds = xr.Dataset() + ds['ismip7BasinNumber'] = basins.astype(float) + ds['ismip7BFRNBin'] = bfrn.astype(float) + ds['ismip7FloatingMask'] = floating.astype(float) + ds['ismip7ShelfRegion'] = region.astype(float) + return ds + + +def _remap_to_mali(config, ds_masks, src_grid_file, mali_mesh_file, + mali_mesh_name, method, ntasks, logger): + """ + Remap the ISMIP7 masks onto the MALI mesh. + + This goes through the shared + :py:func:`compass.landice.ismip7.mapping.build_mapping_file` and + ``ncremap``, the same route ``ismip7_forcing`` uses, rather than through + pyremap's ``Remapper``. pyremap invokes ``mpirun`` directly, which + conflicts with the Slurm allocation on machines where compass launches + with ``srun``; ``build_mapping_file`` uses the configured + ``parallel_executable`` instead. + """ + res = datasets.ISMIP_RESOLUTION_KM + mapping_file = f'map_ismip{res}km_to_{mali_mesh_name}_{method}.nc' + source_file = 'ismip7_masks_source.nc' + remapped_file = 'ismip7_masks_remapped.nc' + + # the masks are assembled in memory, so write them back onto the ISMIP + # grid for ncremap; the source grid file supplies the x/y coordinates + # that the SCRIP description needs + with xr.open_dataset(src_grid_file) as ds_grid: + ds_source = ds_masks.assign_coords(x=ds_grid['x'], y=ds_grid['y']) + write_netcdf(ds_source, source_file) + + logger.info(f'Building mapping file {mapping_file}') + build_mapping_file(config, logger, source_file, mapping_file, + mali_mesh_file=mali_mesh_file, method_remap=method, + projection='ais-bedmap2', ntasks=ntasks) + + logger.info('Remapping the masks onto the MALI mesh') + variables = ','.join(sorted(ds_masks.data_vars)) + check_call(['ncremap', '-i', source_file, '-o', remapped_file, + '-m', mapping_file, '-v', variables], logger=logger) + + ds_remapped = xr.load_dataset(remapped_file) + if 'ncol' in ds_remapped.dims: + ds_remapped = ds_remapped.rename({'ncol': 'nCells'}) + + for path in (source_file, remapped_file): + if os.path.exists(path): + os.remove(path) + + return ds_remapped + + +def _to_integer_masks(ds_remapped): + """ + Round the remapped fields to integers and derive the MALI basin field. + + Nearest-neighbour remapping should already return exact source values, + but they come back as floats; rounding makes the intent explicit and + guards against a different method being used. + """ + ds = xr.Dataset() + + basin0 = ds_remapped['ismip7BasinNumber'] + valid = basin0.notnull() + basin0 = basin0.round().fillna(-1).astype(np.int32) + + ds['ismip7BasinNumber'] = basin0 + ds['ismip7BasinNumber'].attrs = { + 'long_name': 'IMBIE2 drainage basin number, ISMIP7 convention', + 'convention': '0-based, 0-15; basin 9 is Eastern Amundsen, ' + 'basin 14 is Ronne-Filchner'} + + # MALI's melt parameterization expects the 1-based convention + ds['ismip6shelfMelt_basin'] = \ + xr.where(valid, basin0 + 1, 0).astype(np.int32) + ds['ismip6shelfMelt_basin'].attrs = { + 'long_name': 'basin number for the MALI melt parameterization', + 'convention': '1-based, 1-16, equal to ismip7BasinNumber + 1; ' + '0 marks cells with no basin'} + + # MALI reads the per-basin thermal-forcing correction from the same input + # stream as the basin numbers, so it has to be present or the run fails. + # The calibration fits dT_b *after* the melt parameter, and uses zero + # throughout; the fit_delta_t step writes calibrated values over this. + ds['ismip6shelfMelt_deltaT'] = \ + xr.zeros_like(ds['ismip6shelfMelt_basin'], dtype=float) + ds['ismip6shelfMelt_deltaT'].attrs = { + 'long_name': 'basin-wide thermal forcing correction', + 'units': 'degC', + 'note': 'zero as written here; the fit_delta_t step fits one value ' + 'per basin and writes them back over this field'} + + bfrn = ds_remapped['ismip7BFRNBin'] + ds['ismip7BFRNBin'] = bfrn.round().fillna(-1).astype(np.int32) + ds['ismip7BFRNBin'].attrs = { + 'long_name': 'buttressing flux response number bin', + 'convention': '0-9; bin 0 is passive ice, bin 9 the most ' + 'buttressing-relevant; -1 marks cells with no bin'} + + ds['ismip7FloatingMask'] = \ + ds_remapped['ismip7FloatingMask'].round().fillna(0).astype(np.int32) + ds['ismip7FloatingMask'].attrs = { + 'long_name': 'ISMIP7 floating-ice mask', + 'convention': '1 where floating, 0 otherwise', + 'note': 'This is the observed ISMIP7 shelf extent, for diagnostics ' + 'such as comparing modelled with observed shelf area. Melt ' + 'aggregation uses MALI own floating cells, since the ' + 'calibration holds the model accountable for its own shelf ' + 'extent.'} + + region = ds_remapped['ismip7ShelfRegion'].round().fillna(0) + ds['ismip7ShelfRegion'] = region.astype(np.int32) + ds['ismip7ShelfRegion'].attrs = { + 'long_name': 'ice-shelf region for calibration term J4', + 'convention': ', '.join(f'{value} = {name}' + for name, value in REGION_CODES.items())} + return ds + + +def _check_shelves_are_in_their_basin(ds_masks, logger): + """ + Check that Pine Island and Dotson land in the basin the protocol says. + + This is the cheap half of the basin-numbering check, and unlike + :py:func:`_cross_check_basins` it needs no extra input file, so it runs + on any mesh. The protocol refers to the Eastern Amundsen as basin 9 and + Ronne-Filchner as basin 14, both 0-based. If the two numbering + conventions were confused anywhere between the ISMIP7 mask and this + file, the ice shelves would land in the wrong basin. + + Raises + ------ + ValueError + If too few PIG and Dotson cells fall in the Eastern Amundsen + """ + basin = ds_masks['ismip7BasinNumber'].values + region = ds_masks['ismip7ShelfRegion'].values + shelves = region > REGION_CODES['none'] + if not shelves.any(): + raise ValueError( + 'No cells were assigned to Pine Island or Dotson, so term J4 ' + 'would have nothing to aggregate over. Check that the ISMIP7 ' + 'shelf mask covers this mesh.') + + in_basin = basin[shelves] == EASTERN_AMUNDSEN_BASIN + percent = 100.0 * in_basin.mean() + logger.info(f'PIG and Dotson cells in ISMIP7 basin ' + f'{EASTERN_AMUNDSEN_BASIN} (Eastern Amundsen): ' + f'{percent:.1f}% of {int(shelves.sum())}') + + if percent < MIN_SHELF_BASIN_AGREEMENT: + raise ValueError( + f'Only {percent:.1f}% of the Pine Island and Dotson cells fall ' + f'in ISMIP7 basin {EASTERN_AMUNDSEN_BASIN}, the Eastern ' + f'Amundsen, which the protocol says they drain into. The basin ' + f'numbering is probably off: ISMIP7 is 0-based and MALI is ' + f'1-based, so MALI basin 10 is ISMIP7 basin 9.') + + +def _cross_check_basins(ds_masks, region_mask_file, logger): + """ + Compare the remapped basin field with MALI's existing region mask. + + The two are built from different sources -- ISMIP7 IMBIE2 v3 here, versus + the ISMIP6 regions rasterized with ``geometric_features`` in 2022 -- so + exact agreement is not expected. What matters is that the *offset* is + the expected one: an off-by-one would show near-zero agreement + everywhere, whereas genuine boundary differences show up in individual + basins. + + Agreement is reported in **both directions**, because the two masks do + not cover the same cells and a one-directional figure is misleading. A + basin that ISMIP7 draws smaller than ISMIP6 scores high conditioned on + ours and low conditioned on theirs; that is a real difference in the + basin outlines, not an error. + + Raises + ------ + ValueError + If overall agreement is below :py:data:`MIN_BASIN_AGREEMENT`, which + means the numbering conventions are almost certainly mismatched + """ + ds_region = xr.open_dataset(region_mask_file) + if 'regionCellMasks' not in ds_region: + logger.warning('No regionCellMasks in the region mask file; skipping ' + 'the basin cross-check') + return + + masks = ds_region['regionCellMasks'].values + # column index + 1 is MALI's 1-based basin number + mali_basin = np.zeros(masks.shape[0], dtype=np.int32) + for col in range(masks.shape[1]): + mali_basin[masks[:, col] == 1] = col + 1 + + ours = ds_masks['ismip6shelfMelt_basin'].values + floating = ds_masks['ismip7FloatingMask'].values == 1 + both = (ours > 0) & (mali_basin > 0) + if not both.any(): + logger.warning('No cells carry both basin numbers; skipping the ' + 'basin cross-check') + return + + logger.info('') + logger.info('Cross-check against the existing MALI region mask:') + for label, sel in (('all cells', both), + ('floating only', both & floating)): + if not sel.any(): + continue + frac = 100.0 * (ours[sel] == mali_basin[sel]).mean() + logger.info(f' {label:14s} agreement {frac:5.1f}% ' + f'(n={int(sel.sum())})') + + overall = 100.0 * (ours[both] == mali_basin[both]).mean() + + logger.info(' per basin, conditioned on each mask in turn:') + logger.info(f' {"basin":>5s} {"ours->theirs":>13s} {"n":>7s} ' + f'{"theirs->ours":>13s} {"n":>7s}') + for basin in range(1, masks.shape[1] + 1): + sel_ours = both & (ours == basin) + sel_theirs = both & (mali_basin == basin) + if not (sel_ours.any() or sel_theirs.any()): + continue + forward = (100.0 * (mali_basin[sel_ours] == basin).mean() + if sel_ours.any() else float('nan')) + backward = (100.0 * (ours[sel_theirs] == basin).mean() + if sel_theirs.any() else float('nan')) + logger.info(f' {basin:5d} {forward:12.1f}% ' + f'{int(sel_ours.sum()):7d} {backward:12.1f}% ' + f'{int(sel_theirs.sum()):7d}') + logger.info(' (a basin drawn smaller by ISMIP7 than by ISMIP6 scores ' + 'high in the first column and low in the second; that is a ' + 'real difference in the outlines, not an error)') + logger.info('') + + if overall < MIN_BASIN_AGREEMENT: + raise ValueError( + f'The remapped ISMIP7 basins agree with the existing MALI region ' + f'mask for only {overall:.1f}% of cells, below the ' + f'{MIN_BASIN_AGREEMENT:.0f}% threshold. The basin numbering ' + f'conventions are probably mismatched: ISMIP7 is 0-based and ' + f'MALI is 1-based, so MALI basin 10 is ISMIP7 basin 9. Using ' + f'one where the other is expected mis-assigns every basin and ' + f'still produces plausible-looking numbers.') diff --git a/compass/landice/tests/ismip7_calibration/ais/report.py b/compass/landice/tests/ismip7_calibration/ais/report.py new file mode 100644 index 0000000000..b55de94141 --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/ais/report.py @@ -0,0 +1,139 @@ +""" +Report the calibrated parameter distributions and per-term diagnostics. +""" + +import matplotlib +import numpy as np +import xarray as xr + +from compass.landice.tests.ismip7_calibration.configure import parameter_name +from compass.step import Step + +matplotlib.use('Agg') +import matplotlib.pyplot as plt # noqa: E402 + + +class Report(Step): + """ + A step that plots and tabulates the calibration result. + + Attributes + ---------- + melt_forms : list of str + The melt forms that were calibrated + """ + + def __init__(self, test_case, melt_forms): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_calibration.ais.Ais + The test case this step belongs to + + melt_forms : list of str + The melt forms that were calibrated + """ + super().__init__(test_case=test_case, name='report') + self.melt_forms = melt_forms + + for melt_form in melt_forms: + self.add_input_file( + filename=f'calibration_{melt_form}.nc', + target=f'../calibrate/calibration_{melt_form}.nc') + self.add_output_file( + filename=f'parameter_distribution_{melt_form}.png') + self.add_input_file(filename='shelf_area.nc', + target='../aggregate/shelf_area.nc') + self.add_output_file(filename='calibration_summary.txt') + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + lines = [] + + lines.append('ISMIP7 melt-module calibration on the MALI mesh') + lines.append('=' * 62) + lines.append('') + + for melt_form in self.melt_forms: + name = parameter_name(melt_form) + with xr.open_dataset(f'calibration_{melt_form}.nc') as ds: + percentiles = {key: float(ds[key]) + for key in ('p5', 'median', 'p95', 'mode')} + samples = ds['min_p1'].values + values = ds['parameter_values'].values + + lines.append(f'{melt_form}: {name}') + for key in ('p5', 'median', 'p95', 'mode'): + lines.append(f' {key:8s} {percentiles[key]:12.5e}') + lines.append('') + + _plot_distribution(samples, values, percentiles, melt_form, name) + + lines.append('The objective normalises each term by its own median ' + 'over the parameter') + lines.append('ensemble, so the minimised objective is NOT comparable ' + 'between melt forms.') + lines.append('Only parameter values within one form can be compared ' + 'by it.') + lines.append('') + + lines.extend(_shelf_area_lines()) + + text = '\n'.join(lines) + '\n' + with open('calibration_summary.txt', 'w') as handle: + handle.write(text) + for line in lines: + logger.info(line) + + +def _plot_distribution(samples, values, percentiles, melt_form, name): + """Plot the distribution of optimal parameter values.""" + step = np.diff(values).min() + edges = np.append(values - 0.5 * step, values[-1] + 0.5 * step) + + fig, axis = plt.subplots(figsize=(7.0, 4.0), constrained_layout=True) + axis.hist(samples, bins=edges, color='0.7', edgecolor='none') + for key, color, style in (('p5', 'C0', '--'), ('median', 'C3', '-'), + ('p95', 'C0', '--')): + axis.axvline(percentiles[key], color=color, linestyle=style, + label=f'{key} = {percentiles[key]:.4g}') + axis.set_xlabel(f'{name} ({melt_form})') + axis.set_ylabel('draws') + axis.set_title(f'Optimal {name} over 100,000 draws of the term weights\n' + f'and targets, {melt_form} form on the MALI mesh') + axis.legend() + fig.savefig(f'parameter_distribution_{melt_form}.png', dpi=150) + plt.close(fig) + + +def _shelf_area_lines(): + """Tabulate modelled against observed ice-shelf area.""" + lines = ['Ice-shelf area by ISMIP7 basin, 10^3 km^2', + '-' * 62, + f' {"basin":>5s} {"MALI":>10s} {"observed":>10s} ' + f'{"ratio":>8s}'] + with xr.open_dataset('shelf_area.nc') as ds: + modelled = ds['modelled_shelf_area'] + observed = ds['observed_shelf_area'] + for basin in modelled.basins.values: + mali = float(modelled.sel(basins=basin)) + obs = float(observed.sel(basins=basin)) + ratio = mali / obs if obs > 0 else float('nan') + lines.append(f' {int(basin):5d} {mali:10.1f} {obs:10.1f} ' + f'{ratio:8.2f}') + total_mali = float(modelled.sum()) + total_obs = float(observed.sum()) + lines.append(f' {"total":>5s} {total_mali:10.1f} {total_obs:10.1f} ' + f'{total_mali / total_obs:8.2f}') + lines.append('') + lines.append('J1, J2 and J4 are integrated melt, so a shelf-area ' + 'mismatch enters the') + lines.append('calibrated parameter directly. J3 is a basin mean and is ' + 'much less sensitive.') + lines.append('') + return lines diff --git a/compass/landice/tests/ismip7_calibration/ais/run_state.py b/compass/landice/tests/ismip7_calibration/ais/run_state.py new file mode 100644 index 0000000000..aa8cd39bbb --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/ais/run_state.py @@ -0,0 +1,215 @@ +""" +A single-timestep MALI melt diagnostic for one ocean state and melt form. +""" + +import os + +import xarray as xr + +from compass.model import run_model +from compass.step import Step + +#: namelist options that must exist in MALI's defaults for each melt form. +#: Compass only *warns* when an option is missing from the defaults, so +#: without an explicit check a MALI build without the ISMIP7 melt method +#: would silently drop these and produce a plausible but wrong calibration. +REQUIRED_OPTIONS = { + 'ismip7': ['config_ismip7_melt_K', 'config_ismip7_melt_sin_slope', + 'config_ismip7_melt_coriolis', 'config_ismip7_melt_salinity', + 'config_ismip7_melt_salinity_source'], + 'ismip6': ['config_basal_mass_bal_float'], +} + + +class RunState(Step): + """ + A step that runs MALI for a single timestep to diagnose the melt field + for one ocean state, using one melt form. + + Melt is computed from the geometry and the forcing alone, so no ice + dynamics, no thermal evolution and no transient are needed. The velocity + solver is off, which also means Albany is not required. + + The run takes **one short timestep rather than none**, because MALI + computes melt inside the timestep rather than in the initial diagnostic + solve; a zero-length run would produce no melt at all. + + One run per ocean state suffices, not one per parameter value: melt is + exactly proportional to the melt parameter, so the whole parameter + ensemble follows by scaling a single run. That is what turns a ~1300-run + campaign into 28. + + Attributes + ---------- + state_name : str + Name of the ocean state this step runs + + melt_form : str + ``'ismip7'`` for the Burgard local quadratic or ``'ismip6'`` for the + non-local quadratic + """ + + def __init__(self, test_case, state_name, melt_form, subdir, + parameter_scale=1.0): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_calibration.ais.Ais + The test case this step belongs to + + state_name : str + Name of the ocean state to run + + melt_form : {'ismip7', 'ismip6'} + Which melt form to run + + subdir : str + Subdirectory for this step + + parameter_scale : float, optional + Multiple of the reference melt parameter to run at. Only the + linearity check uses anything other than 1. + """ + self.state_name = state_name + self.melt_form = melt_form + self.parameter_scale = parameter_scale + name = os.path.basename(subdir) + super().__init__(test_case=test_case, name=name, subdir=subdir) + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config['ismip7_calibration'] + self.ntasks = section.getint('ntasks') + self.min_tasks = self.ntasks + base_path_mali = section.get('base_path_mali') + mali_mesh_file = section.get('mali_mesh_file') + timestep = section.get('timestep') + + _check_namelist_options(config, self.melt_form) + + self.add_input_file(filename='mesh.nc', + target=os.path.join(base_path_mali, + mali_mesh_file)) + self.add_input_file(filename='masks.nc', + target='../remap_masks/' + 'ismip7_masks_on_mali.nc') + self.add_input_file(filename='forcing.nc', + target=f'../remap_forcing/' + f'forcing_{self.state_name}.nc') + + # MALI builds the partition filename from + # config_block_decomp_file_prefix plus the task count; a symlink + # called graph.info alone is silently not found + self.add_input_file( + filename=f'graph.info.part.{self.ntasks}', + target=f'../make_graph/graph.info.part.{self.ntasks}') + + resource_location = 'compass.landice.tests.ismip7_calibration.ais' + + self.add_namelist_file(resource_location, 'namelist.landice', + out_name='namelist.landice') + + options = {'config_basal_mass_bal_float': f"'{self.melt_form}'", + 'config_dt': f"'{timestep}'", + 'config_run_duration': f"'{timestep}'"} + options.update(_melt_namelist_options(config, self.melt_form, + self.parameter_scale)) + self.add_namelist_options(options=options, + out_name='namelist.landice') + + self.add_streams_file( + resource_location, 'streams.landice.template', + out_name='streams.landice', + template_replacements={'mesh_file': 'mesh.nc', + 'masks_file': 'masks.nc', + 'forcing_file': 'forcing.nc', + 'output_interval': timestep}) + + self.add_model_as_input() + self.add_output_file(filename='output_melt.nc') + + def runtime_setup(self): + """ + Set the number of ocean layers from the forcing file + + MPAS takes dimension sizes from the input stream only, and the mesh + file has no ``nISMIP6OceanLayers``. Without this the 3-D forcing + fields are allocated against a zero-length dimension and the run dies + trying to allocate hundreds of GB. It is read from the forcing + rather than hard-coded, so it stays right if ISMIP7 changes the + vertical grid. + """ + super().runtime_setup() + with xr.open_dataset('forcing.nc') as ds: + n_layers = ds.sizes['nISMIP6OceanLayers'] + self.update_namelist_at_runtime( + options={'config_nISMIP6OceanLayers': f'{n_layers}'}, + out_name='namelist.landice') + + def run(self): + """ + Run this step of the test case + """ + # the partition file is supplied ready-made alongside the mesh, so + # there is no plain graph.info for gpmetis to partition + run_model(self, partition_graph=False) + + +def _melt_namelist_options(config, melt_form, parameter_scale=1.0): + """The namelist options specific to one melt form.""" + section = config['ismip7_calibration_melt'] + if melt_form == 'ismip7': + # the run is done at a reference parameter value; melt is exactly + # proportional to it, so the ensemble is formed by scaling afterwards + melt_k = section.getfloat('reference_k') * parameter_scale + return { + 'config_ismip7_melt_K': repr(melt_k), + 'config_ismip7_melt_sin_slope': + repr(section.getfloat('sin_slope')), + 'config_ismip7_melt_coriolis': + repr(section.getfloat('coriolis')), + 'config_ismip7_melt_salinity_source': "'constant'", + 'config_ismip7_melt_salinity': + repr(section.getfloat('salinity'))} + # the ISMIP6 non-local form reads gamma0 from its input stream, so the + # reference value is written into the masks file rather than set here + return {} + + +def _check_namelist_options(config, melt_form): + """ + Check that MALI's default namelist has the options this melt form needs. + + Raises + ------ + ValueError + If any required option is missing, which means the MALI build does + not support this melt form + """ + defaults = config.get('namelists', 'forward') + if not os.path.exists(defaults): + raise FileNotFoundError( + f'MALI default namelist not found at {defaults}. Build MALI, or ' + f'set [paths] mpas_model (or [namelists] forward) in your config ' + f'file.') + + with open(defaults) as handle: + text = handle.read() + + missing = [option for option in REQUIRED_OPTIONS[melt_form] + if option not in text] + if missing: + raise ValueError( + f"MALI's default namelist at\n {defaults}\ndoes not contain " + f"{', '.join(missing)}, so this build does not support " + f"config_basal_mass_bal_float = '{melt_form}'.\n" + f"Compass only warns about namelist options it cannot find, so " + f"without this check the run would silently use MALI's defaults " + f"and produce a plausible but wrong calibration.\n" + f"Build MALI from a branch that has the ISMIP7 melt " + f"parameterization and point [paths] mpas_model at it.") diff --git a/compass/landice/tests/ismip7_calibration/ais/streams.landice.template b/compass/landice/tests/ismip7_calibration/ais/streams.landice.template new file mode 100644 index 0000000000..9f961c3ac8 --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/ais/streams.landice.template @@ -0,0 +1,57 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/compass/landice/tests/ismip7_calibration/ais/verify_melt.py b/compass/landice/tests/ismip7_calibration/ais/verify_melt.py new file mode 100644 index 0000000000..a45a1c860b --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/ais/verify_melt.py @@ -0,0 +1,256 @@ +""" +Verify MALI's melt against the Python reference implementation. +""" + +import numpy as np +import xarray as xr +from mpas_tools.io import write_netcdf + +from compass.landice.tests.ismip7_calibration.ais import melt_model +from compass.step import Step + +#: the melt expression should agree to round-off +MELT_TOLERANCE = 1.0e-10 + +#: so should the vertical interpolation of thermal forcing, in K +INTERPOLATION_TOLERANCE = 1.0e-10 + +#: melt should be exactly proportional to the melt parameter; allow only +#: floating-point noise +LINEARITY_TOLERANCE = 1.0e-12 + + +class VerifyMelt(Step): + """ + A step that checks MALI's melt field against an independent Python + implementation of the same equations. + + Three things are checked, and they are deliberately independent: + + * the **melt expression**, by evaluating the Python reference on MALI's + *own* ``TFdraft``. That isolates the formula from the vertical + interpolation that produced ``TFdraft``. + * the **vertical interpolation**, by interpolating the 3-D forcing to the + ice draft with a plain ``numpy`` implementation written from the + protocol rather than transliterated from the Fortran, and comparing + against MALI's ``TFdraft``. + * the **linearity in the melt parameter**, which is what licenses one + MALI run per ocean state instead of one per (state, parameter) pair. + + The draft used for the interpolation check is reconstructed from the + mesh file, never taken from the run output: melt thins the ice over the + single timestep, so the output geometry is post-step while ``TFdraft`` + was computed pre-step. + + Attributes + ---------- + melt_form : str + The melt form being verified + + state_name : str + The ocean state used for the check + """ + + def __init__(self, test_case, melt_form, state_name, + linearity_scales=()): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_calibration.ais.Ais + The test case this step belongs to + + melt_form : {'ismip7', 'ismip6'} + The melt form to verify + + state_name : str + The ocean state to verify against + + linearity_scales : sequence of float, optional + Multiples of the reference melt parameter that extra MALI runs + were done at, for the linearity check + """ + super().__init__(test_case=test_case, name='verify_melt') + self.melt_form = melt_form + self.state_name = state_name + self.linearity_scales = tuple(linearity_scales) + + for scale in self.linearity_scales: + self.add_input_file( + filename=f'linearity_{scale:g}.nc', + target=f'../{melt_form}_{state_name}_x{scale:g}/' + f'output_melt.nc') + + self.add_input_file( + filename='output_melt.nc', + target=f'../{melt_form}_{state_name}/output_melt.nc') + self.add_input_file( + filename='forcing.nc', + target=f'../remap_forcing/forcing_{state_name}.nc') + self.add_output_file(filename='verification.nc') + + def setup(self): + """ + Set up this step of the test case + """ + section = self.config['ismip7_calibration'] + self.add_input_file( + filename='mesh.nc', + target=(f'{section.get("base_path_mali")}/' + f'{section.get("mali_mesh_file")}')) + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + section = config['ismip7_calibration_melt'] + reference = {'ismip7': section.getfloat('reference_k'), + 'ismip6': section.getfloat('reference_gamma0')} + + fields = melt_model.read_run('output_melt.nc', + 'mesh.nc') + results = {} + + results['melt'] = _check_melt_expression( + self.melt_form, reference[self.melt_form], fields, config, + logger) + results['interpolation'] = _check_interpolation(fields, logger) + results['linearity'] = _check_linearity( + fields, self.linearity_scales, logger) + + ds = xr.Dataset({name: float(value) + for name, value in results.items()}) + ds.attrs['melt_form'] = self.melt_form + ds.attrs['ocean_state'] = self.state_name + write_netcdf(ds, 'verification.nc') + + _raise_on_failure(results, logger) + + +def _check_melt_expression(melt_form, parameter, fields, config, logger): + """Compare MALI's melt with the Python reference on MALI's TFdraft.""" + expected = melt_model.melt_from_tf(melt_form, parameter, fields, config) + actual = fields['melt'] + + melting = fields['floating'] & (np.abs(actual) > 0.0) + count = int(melting.sum()) + if count == 0: + raise ValueError('MALI produced no melt anywhere, so the melt ' + 'expression cannot be checked. Check that the run ' + 'took a real timestep and that the forcing was ' + 'read.') + + difference = np.abs(actual - expected).where(melting) + scale = np.abs(actual).where(melting) + relative = float((difference / scale).max()) + + logger.info('') + logger.info(f'Melt expression, {melt_form}:') + logger.info(f' melting cells {count}') + logger.info(f' MALI melt range ' + f'{float(actual.where(melting).min()):.4g} .. ' + f'{float(actual.where(melting).max()):.4g} kg/m2/yr') + logger.info(f' max relative difference {relative:.3e}') + return relative + + +def _check_interpolation(fields, logger): + """Compare MALI's TFdraft with an independent interpolation.""" + with xr.open_dataset('forcing.nc') as ds: + tf_3d = ds['ismip6shelfMelt_3dThermalForcing'].isel(Time=0).values + z_ocean = ds['ismip6shelfMelt_zOcean'].values + + draft = melt_model.initial_draft('mesh.nc').values + with xr.open_dataset('mesh.nc') as ds_mesh: + bed = ds_mesh['bedTopography'] + if 'Time' in bed.dims: + bed = bed.isel(Time=0) + bed = bed.values + + floating = fields['floating'].values + expected = melt_model.interpolate_to_draft(tf_3d, z_ocean, draft, bed, + freezing_correction=True) + actual = fields['tf_draft'].values + + # all four of MALI's code paths are covered, so nothing is excluded + difference = np.abs(actual[floating] - expected[floating]) + largest = float(np.nanmax(difference)) if difference.size else 0.0 + + logger.info('') + logger.info('Vertical interpolation of thermal forcing:') + logger.info(f' cells compared {int(floating.sum())}') + logger.info(f' max |difference| {largest:.3e} K') + return largest + + +def _check_linearity(fields, scales, logger): + """ + Check that MALI's melt is exactly proportional to the melt parameter. + + This is measured with real MALI runs rather than argued, because it is + what licenses one run per ocean state instead of one per (state, + parameter) pair -- 28 runs rather than about 1300. + + Each extra run is the same ocean state at a different multiple of the + reference melt parameter. Dividing each total by its multiple must give + the same number every time. + """ + if not scales: + logger.info('') + logger.info('Linearity in the melt parameter: no scaled runs were ' + 'set up, so this is not measured. The melt parameter ' + 'enters the melt expression only as a multiplicative ' + 'coefficient, so linearity follows from the code, but ' + 'measuring it is better.') + return 0.0 + + area = fields['area'] + floating = fields['floating'] + + def total(melt): + return float((melt * area).where(floating).sum()) / 1.0e12 + + totals = {1.0: total(fields['melt'])} + for scale in scales: + other = melt_model.read_run(f'linearity_{scale:g}.nc', + 'mesh.nc') + totals[scale] = total(other['melt']) + + per_unit = {scale: value / scale for scale, value in totals.items()} + values = list(per_unit.values()) + spread = (max(values) - min(values)) / abs(per_unit[1.0]) + + logger.info('') + logger.info('Linearity in the melt parameter, measured in MALI:') + logger.info(f' {"scale":>8s} {"total melt":>14s} {"total / scale":>16s}') + for scale in sorted(totals): + logger.info(f' {scale:8.2f} {totals[scale]:14.4f} ' + f'{per_unit[scale]:16.6f}') + logger.info(f' max relative deviation {spread:.3e}') + logger.info('') + return spread + + +def _raise_on_failure(results, logger): + """Raise if any check exceeded its tolerance.""" + checks = (('melt', MELT_TOLERANCE, + "MALI's melt does not match the Python reference evaluated " + "on MALI's own TFdraft"), + ('interpolation', INTERPOLATION_TOLERANCE, + "MALI's vertical interpolation of thermal forcing to the ice " + "draft does not match an independent implementation"), + ('linearity', LINEARITY_TOLERANCE, + 'melt is not exactly proportional to the melt parameter, so ' + 'the ensemble cannot be built by scaling a single run per ' + 'ocean state')) + + failures = [f' {name}: {results[name]:.3e} exceeds {tol:.0e} -- {why}' + for name, tol, why in checks if results[name] > tol] + if failures: + listing = '\n'.join(failures) + raise ValueError(f'MALI melt verification failed:\n{listing}') + + logger.info('All melt verification checks passed.') diff --git a/compass/landice/tests/ismip7_calibration/configure.py b/compass/landice/tests/ismip7_calibration/configure.py new file mode 100644 index 0000000000..872113071a --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/configure.py @@ -0,0 +1,180 @@ +""" +Shared configuration checks for the ISMIP7 calibration test cases. +""" + +import numpy as np + +from compass.landice.tests.ismip7_calibration.datasets import ( + PUBLISHED_T3_MODELS, + PUBLISHED_T4_REGIONS, + PUBLISHED_T4_YEARS, +) + + +def check_options(config, options): + """ + Check that user-supplied config options have been set. + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for the test case + + options : list of str + Options in the ``ismip7_calibration`` section that the user must + supply + + Raises + ------ + ValueError + If any of ``options`` is still ``NotAvailable`` + """ + section = 'ismip7_calibration' + for option in options: + value = config.get(section=section, option=option) + if value == 'NotAvailable': + raise ValueError( + f'You need to supply a user config file containing the ' + f'[{section}] section with the {option} option set.') + + +def parameter_values(config, melt_form): + """ + The grid of parameter values to select from, for one melt form. + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for the test case + + melt_form : {'ismip7', 'ismip6'} + Which melt form the parameter belongs to. ``'ismip7'`` calibrates + ``K`` and ``'ismip6'`` calibrates ``gamma0``. + + Returns + ------- + values : numpy.ndarray + The parameter grid + """ + section = config['ismip7_calibration_melt'] + if melt_form == 'ismip7': + low = section.getfloat('k_min') + high = section.getfloat('k_max') + step = section.getfloat('k_step') + elif melt_form == 'ismip6': + low = section.getfloat('gamma0_min') + high = section.getfloat('gamma0_max') + step = section.getfloat('gamma0_step') + else: + raise ValueError(f"melt_form must be 'ismip7' or 'ismip6', but is " + f"'{melt_form}'") + + # add half a step so that ``high`` itself is included + return np.arange(low, high + 0.5 * step, step) + + +def parameter_name(melt_form): + """ + The name of the parameter a melt form calibrates. + + Parameters + ---------- + melt_form : {'ismip7', 'ismip6'} + The melt form + + Returns + ------- + name : str + ``'K'`` for ``'ismip7'`` and ``'gamma0'`` for ``'ismip6'`` + """ + if melt_form == 'ismip7': + return 'K' + if melt_form == 'ismip6': + return 'gamma0' + raise ValueError(f"melt_form must be 'ismip7' or 'ismip6', but is " + f"'{melt_form}'") + + +def melt_forms(config): + """ + The melt forms to calibrate. + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for the test case + + Returns + ------- + forms : list of str + Each of ``'ismip7'`` and ``'ismip6'`` that was requested + """ + value = config.get('ismip7_calibration', 'melt_forms') + forms = [form.strip() for form in value.split(',') if form.strip()] + for form in forms: + if form not in ('ismip7', 'ismip6'): + raise ValueError(f"melt_forms must contain only 'ismip7' and " + f"'ismip6', but contains '{form}'") + if not forms: + raise ValueError('melt_forms must name at least one melt form') + return forms + + +def weighting(config): + """ + Which summands of J3 and J4 carry non-zero weight. + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for the test case + + Returns + ------- + t3_models : tuple of str or None + Ocean models to weight in J3, or None to weight all of them + + t4_regions : tuple of str or None + Ice shelves to weight in J4, or None to weight all of them + + t4_years : tuple of int or None + Observation years to weight in J4, or None to weight all of them + """ + section = config['ismip7_calibration_objective'] + + def _select(option, published): + value = section.get(option).strip() + if value == 'published': + return published + if value == 'all': + return None + raise ValueError(f"{option} must be 'published' or 'all', but is " + f"'{value}'") + + return (_select('t3_models', PUBLISHED_T3_MODELS), + _select('t4_regions', PUBLISHED_T4_REGIONS), + _select('t4_years', PUBLISHED_T4_YEARS)) + + +def objective_options(config): + """ + The sample size and random seed for the parameter selection. + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for the test case + + Returns + ------- + sample_size : int + Number of random draws + + seed : int or None + Seed for the random draws, or None to leave the random state alone + """ + section = config['ismip7_calibration_objective'] + sample_size = section.getint('sample_size') + seed = section.get('seed').strip() + seed = None if seed.lower() == 'none' else int(seed) + return sample_size, seed diff --git a/compass/landice/tests/ismip7_calibration/datasets.py b/compass/landice/tests/ismip7_calibration/datasets.py new file mode 100644 index 0000000000..ed7a38ec79 --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/datasets.py @@ -0,0 +1,286 @@ +""" +Registry of the ISMIP7 datasets used by the melt-module calibration. + +One source of truth for *which* ocean states feed *which* objective-function +term, and for where the calibration targets live, shared by the 8 km +``replication`` test case and the MALI-mesh ``ais`` test case so the two +cannot drift apart. + +The ocean-state sets follow protocol Table 2 and the focus group's 31 July +2026 update: + +* ``minimal`` -- the states marked X in Table 2: present-day climatology, two + circum-Antarctic cold/warm pairs, and PIG in 2009 and 2012. 7 states. +* ``recommended`` -- adds the two regional cold/warm pairs marked +, which the + focus group strongly suggests. 11 states. +* ``all`` -- every state distributed for ISMIP7: all seven ocean models of + protocol Table 2 (AIS1-AIS7) and all thirteen observation years. 28 states. + This is what the ``model`` coordinate of the J3 target file expects, so only + this subset exercises every target the protocol ships. +""" + +import os +from collections import namedtuple + +import numpy as np +import pandas as pd +import xarray as xr + +#: ISMIP grid resolution of the calibration datasets, in km +ISMIP_RESOLUTION_KM = 8 + +#: ice-shelf ids in ``shelf_mask_ismip8km.nc`` +PIG_ID = 110 +DOTSON_ID = 97 + +#: Pine Island is cut at this x to keep only its main trunk, following the +#: protocol's worked example +PIG_X_MAX = -1.625e6 + +#: ocean-model datasets: file prefix, toolbox label, basins constrained. +#: ``basins`` is None where the simulation is circum-Antarctic; otherwise it +#: lists the ISMIP7 (0-based) basins the regional domain actually covers, and +#: melt outside them must be discarded. +OCEAN_MODELS = [ + ('Mathiot_NEMO_{state}_v3_', 'mathiot', None), + ('Timmermann_FESOM_{state}_v3_', 'timmermann', (14,)), + ('Naughten_FESOM_ACCESS_{state}_v2_', 'naughten_ais_1', None), + ('Naughten_FESOM_MMM_{state}_v2_', 'naughten_ais_2', None), + ('Jourdain-Naughten_NEMO-MITgcm_{state}_', 'jourdain_naughten', (9, 14)), + ('Naughten_MITamu-MITwed_{state}_', 'naughten_naughten', (9, 14)), + ('Haid_FESOM_{state}_', 'haid', (14,)), +] + +#: all years of Amundsen near-ice-shelf ocean observations +OBS_YEARS = [1994, 2000, 2006, 2007, 2009, 2010, 2011, 2012, 2014, 2016, + 2018, 2019, 2020] + +#: minimal J3 set (Table 2, marked X) +MINIMAL_MODELS = ('mathiot', 'naughten_ais_1') + +#: J3 set the focus group recommends (adds those marked +) +RECOMMENDED_MODELS = ('mathiot', 'naughten_ais_1', 'jourdain_naughten', + 'naughten_naughten') + +#: J4 years the protocol suggests as a minimum, a cold and a warm PIG state +RECOMMENDED_OBS_YEARS = (2009, 2012) + +#: the weighting the published quadratic example used: four ocean models, and +#: PIG alone in 2009 and 2012. Narrower than protocol Table 2 makes available +PUBLISHED_T3_MODELS = RECOMMENDED_MODELS +PUBLISHED_T4_REGIONS = ('pig',) +PUBLISHED_T4_YEARS = RECOMMENDED_OBS_YEARS + +#: One ocean state to force the melt module with. +#: +#: ``name`` is a unique, filesystem-safe identifier such as ``mathiot_cold``; +#: ``kind`` is ``'climatology'``, ``'model'`` or ``'obs'``; ``term`` names the +#: objective-function term the state feeds; ``label`` and ``state`` apply to +#: ocean-model states, ``year`` to observational ones; and ``basins`` lists +#: the ISMIP7 (0-based) basins a regional state constrains, or None for a +#: circum-Antarctic one. +OceanState = namedtuple( + 'OceanState', + ['name', 'kind', 'tf_file', 'so_file', 'term', 'label', 'state', 'year', + 'basins']) + + +def climatology_files(base_path): + """ + Paths to the Zhou et al. present-day climatology TF and salinity. + + Parameters + ---------- + base_path : str + Root of the ISMIP7 AIS datasets + + Returns + ------- + tf_file : str + The thermal-forcing file + + so_file : str + The salinity file + """ + clim = os.path.join(base_path, 'obs', 'ocean', 'climatology', + 'zhou_annual_06_nov') + stem = 'AIS_obs_ocean_climatology_zhou_annual_06_nov' + return (os.path.join(clim, 'tf', 'v3', f'tf_{stem}_v3_1972-2024.nc'), + os.path.join(clim, 'so', 'v4', f'so_{stem}_v4_1972-2024.nc')) + + +def ocean_states(base_path, subset='all'): + """ + The ocean states to run the melt module for. + + Parameters + ---------- + base_path : str + Root of the ISMIP7 AIS datasets + + subset : {'minimal', 'recommended', 'all'}, optional + Which set of states to include; see the module docstring + + Returns + ------- + states : list of OceanState + The ocean states in the requested subset + """ + if subset not in ('minimal', 'recommended', 'all'): + raise ValueError(f"subset must be 'minimal', 'recommended' or 'all', " + f"but is '{subset}'") + + param = os.path.join(base_path, 'parameterisations', 'ocean') + model_dir = os.path.join(param, 'ocean_modelling_data') + obs_dir = os.path.join(param, 'ocean_observations_data') + + if subset == 'minimal': + keep_models = MINIMAL_MODELS + keep_years = RECOMMENDED_OBS_YEARS + elif subset == 'recommended': + keep_models = RECOMMENDED_MODELS + keep_years = RECOMMENDED_OBS_YEARS + else: + keep_models = tuple(label for _, label, _ in OCEAN_MODELS) + keep_years = tuple(OBS_YEARS) + + tf_file, so_file = climatology_files(base_path) + states = [OceanState(name='climatology', kind='climatology', + tf_file=tf_file, so_file=so_file, term='J1,J2', + label=None, state=None, year=None, basins=None)] + + for prefix, label, basins in OCEAN_MODELS: + if label not in keep_models: + continue + for state in ('cold', 'warm'): + stem = prefix.format(state=state) + states.append(OceanState( + name=f'{label}_{state}', kind='model', + tf_file=os.path.join(model_dir, f'{stem}TF.nc'), + so_file=os.path.join(model_dir, f'{stem}S.nc'), + term='J3', label=label, state=state, year=None, + basins=basins)) + + for year in OBS_YEARS: + if year not in keep_years: + continue + states.append(OceanState( + name=f'obs_{year}', kind='obs', + tf_file=os.path.join(obs_dir, f'Obs_{year}_TF.nc'), + so_file=os.path.join(obs_dir, f'Obs_{year}_S.nc'), + term='J4', label=None, state=None, year=year, basins=(9,))) + + return states + + +def missing_files(states): + """ + The input files of ``states`` that are not present on disk. + + Parameters + ---------- + states : list of OceanState + The ocean states to check + + Returns + ------- + missing : list of tuple + ``(state name, path)`` pairs for each file that does not exist + """ + missing = [] + for state in states: + for path in (state.tf_file, state.so_file): + if not os.path.exists(path): + missing.append((state.name, path)) + return missing + + +def mask_files(base_path, resolution_km=ISMIP_RESOLUTION_KM): + """ + Paths to the ISMIP7 mask datasets on the ISMIP polar stereographic grid. + + Parameters + ---------- + base_path : str + Root of the ISMIP7 AIS datasets + + resolution_km : int, optional + ISMIP grid resolution in km + + Returns + ------- + files : dict + Paths keyed by ``basins``, ``bfrn``, ``floating`` and ``shelves`` + """ + param = os.path.join(base_path, 'parameterisations', 'ocean') + res = resolution_km + return dict( + basins=os.path.join(param, 'imbie2', + f'basin_numbers_ismip{res}km_v2.nc'), + bfrn=os.path.join(param, 'bfrns', f'BFRN_ismip{res}km_v2.nc'), + floating=os.path.join(param, 'floatingmasks', + f'floatingmask_ismip{res}km.nc'), + shelves=os.path.join(param, 'shelfmask', + f'shelf_mask_ismip{res}km.nc')) + + +def load_targets(base_path): + """ + Load the observational targets for the four objective-function terms. + + These are on the ISMIP grid or are already aggregated, so they are the + same whether the modelled melt came from the 8 km reference implementation + or from MALI. + + Parameters + ---------- + base_path : str + Root of the ISMIP7 AIS datasets + + Returns + ------- + targets : dict + ``t1_mean``, ``t1_sigma``, ``t2_mean``, ``t2_sigma``, ``t2_weights``, + ``t3_mean``, ``t3_sigma``, ``t4_mean`` and ``t4_sigma`` + """ + param = os.path.join(base_path, 'parameterisations', 'ocean') + + melt_imbie = pd.read_csv( + os.path.join(param, 'meltobs', + 'Melt_Paolo_Davison_Adusumilli_imbie2.csv'), + index_col=0) + basin_coord = np.arange(len(melt_imbie)) + t1_mean = xr.DataArray( + melt_imbie['BMR (Gt/yr)'].values.astype(float), + dims=['basin'], coords={'basin': basin_coord}) + t1_sigma = xr.DataArray( + melt_imbie['BMR uncert (Gt/yr)'].values.astype(float), + dims=['basin'], coords={'basin': basin_coord}) + + buttressing = xr.load_dataset( + os.path.join(param, 'meltobs', 'melt_target_term2_v3.nc')) + bfrn = xr.load_dataset(mask_files(base_path)['bfrn']) + t2_weights = xr.DataArray( + (bfrn['BFRN_medians'] / bfrn['BFRN_median']).values, + dims=['BFRN_bins'], + coords={'BFRN_bins': buttressing.BFRN_bins.values}) + + cold = xr.load_dataset( + os.path.join(param, 'ocean_modelling_data', + 'melt_cold_target_term3_v2.nc')) + warm = xr.load_dataset( + os.path.join(param, 'ocean_modelling_data', + 'melt_warm_target_term3_v2.nc')) + t3_mean = warm.melt_rate - cold.melt_rate + t3_sigma = np.sqrt(warm.melt_rate_uncert**2 + cold.melt_rate_uncert**2) + + t4 = xr.load_dataset( + os.path.join(param, 'ocean_observations_data', + 'melt_observations_target_term4.nc')) + + return dict(t1_mean=t1_mean, t1_sigma=t1_sigma, + t2_mean=buttressing['melt_mean'], + t2_sigma=buttressing['melt_mean_err'], + t2_weights=t2_weights, + t3_mean=t3_mean, t3_sigma=t3_sigma, + t4_mean=t4.melt_rate, t4_sigma=t4.melt_rate_uncert) diff --git a/compass/landice/tests/ismip7_calibration/ismip7_calibration.cfg b/compass/landice/tests/ismip7_calibration/ismip7_calibration.cfg new file mode 100644 index 0000000000..1bdffea94b --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/ismip7_calibration.cfg @@ -0,0 +1,157 @@ +# config options for the ISMIP7 melt-module calibration +[ismip7_calibration] + +# Base path to the ISMIP7 AIS datasets, the directory containing +# parameterisations/, obs/ and grid/. User has to supply. +base_path_ismip7 = NotAvailable + +# Base path to the MALI mesh and its ancillary files. User has to supply. +base_path_mali = NotAvailable + +# MALI mesh file, holding the present-day geometry to calibrate on. The +# protocol asks for present-day geometry, and the E3SM inputdata meshes are +# pre-relaxation and already essentially observed, so no geometry replacement +# is needed. User has to supply. +mali_mesh_file = NotAvailable + +# Name of the MALI mesh. Used to name mapping and output files. +mali_mesh_name = NotAvailable + +# MALI region mask file, used to cross-check the remapped ISMIP7 basins +# against the ISMIP6-era regions. Set to None to skip the cross-check, but +# note that this check is what catches the ISMIP7/MALI basin numbering being +# off by one. +region_mask_file = NotAvailable + +# Which ocean states to use. Options: +# minimal - 7 states, those marked X in protocol Table 2 +# recommended - 11 states, adding the regional pairs marked + +# all - 28 states, everything ISMIP7 distributes +# The reduced sets leave J4 structurally under-powered: with PIG alone it +# constrains a single amplitude that any form can match by rescaling its +# parameter. +ocean_state_subset = all + +# Which melt forms to calibrate, comma separated. Options: +# ismip7 - the Burgard et al. (2022) local quadratic, calibrating K +# ismip6 - the ISMIP6 non-local quadratic, calibrating gamma0. With a +# constant salinity this is algebraically identical to the +# Burgard semi-local form, so it is how the semi-local form is +# calibrated. +melt_forms = ismip7, ismip6 + +# Number of MPI tasks for ESMF_RegridWeightGen +esmf_ntasks = 128 + +# Number of MPI tasks for each MALI run +ntasks = 128 + +# Length of the single timestep the melt diagnostic takes. MALI computes +# melt inside the timestep rather than in the initial diagnostic solve, so a +# zero-length run would produce no melt at all. +timestep = 0000-00-01_00:00:00 + + +# config options for remapping the ISMIP7 masks onto the MALI mesh +[ismip7_calibration_masks] + +# Remapping method. Every mask field is categorical, so nearest neighbour is +# the appropriate choice. Options: bilinear, neareststod, conserve +method_remap = neareststod + + +# config options for remapping the ISMIP7 thermal forcing onto the MALI mesh +[ismip7_calibration_forcing] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = bilinear + + +# config options for the melt parameterizations +[ismip7_calibration_melt] + +# Practical salinity at the ice draft, PSU. The ISMIP7 ocean forcing +# processed for the MALI projections carries thermal forcing only, so a +# projection has no salinity field to read; calibrating against a spatially +# varying salinity would tune K for physics the projections cannot run. +# Melt is linear in salinity, so this is a coefficient multiplying K. +salinity = 34.5 + +# The sin(theta) factor of the ISMIP7 quadratic, where theta is the ice-draft +# slope angle. Held constant, so it too simply multiplies K. The default is +# the Antarctic-mean value from Bedmap3 on the ISMIP 8 km grid, matching the +# published calibration. +sin_slope = 0.0051117 + +# Magnitude of the Coriolis parameter, s^-1. Constant, as in the protocol's +# reference implementation. +coriolis = 1.4e-4 + +# Reference parameter values the ensemble is run at. Melt is exactly +# proportional to the melt parameter, so one run per ocean state at a +# reference value is enough and the whole parameter grid follows by scaling. +# The reference values are divided out again, so they do not affect results. +reference_k = 1.0e-4 +reference_gamma0 = 10000.0 + +# The K grid for the ISMIP7 local form: minimum, maximum and step. The +# defaults are the 120 values the protocol's worked example samples. +k_min = 0.25e-5 +k_max = 3.0e-4 +k_step = 0.25e-5 + +# The gamma0 grid for the ISMIP6 non-local form, m yr^-1: minimum, maximum +# and step. The maximum has to be high enough that no draw of the objective +# lands on it: a distribution that piles up against the end of the grid has +# been clipped, and its upper percentiles are not meaningful. +gamma0_min = 250.0 +gamma0_max = 50000.0 +gamma0_step = 250.0 + + +# config options for the objective function and the parameter selection +[ismip7_calibration_objective] + +# Number of random draws of the term weights and the targets +sample_size = 100000 + +# Seed for the random draws, so that the reported percentiles are +# reproducible. Set to None to leave the random state alone. +seed = 0 + +# Which ocean models carry non-zero weight in J3. Options: +# published - the four models of the published quadratic example +# all - every model in the ensemble +t3_models = all + +# Which ice shelves carry non-zero weight in J4. Options: +# published - PIG alone, as the published quadratic example weights it +# all - PIG and Dotson +# The published weighting uses 2 of the 18 available observations and cannot +# discriminate between melt forms, because with one shelf J4 constrains a +# single amplitude that either form matches by rescaling its parameter. +t4_regions = all + +# Which observation years carry non-zero weight in J4. Options: +# published - 2009 and 2012, as the published quadratic example weights it +# all - every year in the ensemble +t4_years = all + + +# config options for the basin thermal-forcing correction, dT_b +[ismip7_calibration_delta_t] + +# Whether to fit dT_b after the parameter selection. Protocol Sect. 4.2.1 +# permits the correction to be fitted before or after; this follows option 2, +# after, as the published quadratic example does. Fitting it first would +# leave the parameter bounds unconstrained where present-day melt is compared +# to observations, which is the drawback Sect. 4.2.1 names. +fit_delta_t = True + +# Bounds on dT_b, in K. The protocol recommends bounding it, since the +# warmest-to-coldest spread across Antarctic ice shelves is only 3-4 K. +delta_t_min = -2.0 +delta_t_max = 2.0 + +# Number of points in the dT_b grid search +delta_t_count = 81 diff --git a/compass/landice/tests/ismip7_calibration/objective.py b/compass/landice/tests/ismip7_calibration/objective.py new file mode 100644 index 0000000000..1a2ecfe60d --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/objective.py @@ -0,0 +1,221 @@ +""" +Driving the ISMIP7 objective function. + +The vendored toolbox +(:py:mod:`compass.landice.tests.ismip7_calibration.toolbox`) implements +protocol Eq. (3): it draws random term weights and random targets within +their uncertainties, and for each draw picks the parameter value that +minimises ``I = sum_i a_i J_i / median(J_i)``. The distribution of those +minimisers gives the 5th / 50th / 95th percentiles the ISMIP7 projections +need. + +This module assembles the arguments the toolbox expects and reduces its +output to those percentiles. It is shared by the ``replication`` test case, +which drives it with melt from the Python reference implementation on the +ISMIP 8 km grid, and by the ``ais`` test case, which drives it with melt +from MALI on a MALI mesh. Keeping one implementation means the MALI +calibration is validated by the replication reproducing published numbers. + +**Melt is exactly linear in the calibration parameter**, so each term is +built once at a unit parameter value and scaled onto the parameter grid. +That is exact, not an approximation, and it is why the ensemble needs one +model run per ocean state rather than one per (state, parameter) pair. +""" + +import numpy as np +import xarray as xr + +from compass.landice.tests.ismip7_calibration.toolbox import ( + parameter_selection_toolbox as toolbox, +) + + +def scale_to_ensemble(unit_aggregate, param_values): + """ + Turn a unit-parameter aggregate into the ``(p1, p2, ...)`` ensemble. + + The toolbox wants each modelled term indexed by ``(p1, p2, ...)``, where + ``p1`` is the melt parameter being selected and ``p2`` a second parameter + that the quadratic parameterizations do not use, so it is a singleton. + + Parameters + ---------- + unit_aggregate : xarray.DataArray + An aggregate (J1, J2, J3 or J4) computed at a parameter value of 1 + + param_values : numpy.ndarray + The parameter grid to scale onto + + Returns + ------- + scaled : xarray.DataArray + ``unit_aggregate`` broadcast over ``p1`` and a singleton ``p2`` + """ + p1 = xr.DataArray(param_values, dims=['p1'], + coords={'p1': param_values}) + out = (unit_aggregate * p1).expand_dims({'p2': np.ones(1)}) + return out.transpose('p1', 'p2', ...) + + +def build_toolbox_terms(units, targets, param_values, t3_models=None, + t4_regions=None, t4_years=None): + """ + Assemble the arguments ``calculate_objective_function`` expects. + + Parameters + ---------- + units : dict + Unit-parameter aggregates ``t1``, ``t2``, ``t3`` and ``t4``, as + produced by :py:mod:`compass.landice.tests.ismip7_calibration.terms` + + targets : dict + Observational targets, from + :py:func:`compass.landice.tests.ismip7_calibration.datasets.load_targets` + + param_values : numpy.ndarray + The parameter grid + + t3_models : sequence of str, optional + Ocean models to give non-zero weight in J3. None weights every model + the ensemble provides. + + t4_regions : sequence of str, optional + Ice shelves to give non-zero weight in J4, from ``'pig'`` and + ``'dotson'``. None weights both. + + The published weighting uses PIG alone, which is 2 of the 18 + available observations. As weighted that way, J4 constrains a single + amplitude that either melt form can match by rescaling its parameter, + so it cannot discriminate between forms; including Dotson makes it a + relative constraint between two shelves, which it can. + + t4_years : sequence of int, optional + Observation years to give non-zero weight in J4. None weights every + year the ensemble provides. + + Returns + ------- + terms : dict + Keyword arguments for + ``toolbox.calculate_objective_function`` + """ # noqa: E501 + t1_model = scale_to_ensemble(units['t1'], param_values) + t2_model = scale_to_ensemble(units['t2'], param_values) + t3_model = scale_to_ensemble(units['t3'], param_values) + + t4_model = scale_to_ensemble(units['t4'], param_values) + t4_model = t4_model.where(targets['t4_mean'].notnull()) + t4_model = t4_model.reindex_like(targets['t4_mean']) + + t1_weights = xr.DataArray( + np.ones(t1_model.sizes['basins']), dims=['basins'], + coords={'basins': t1_model.basins.values}) + + t3_weights = xr.DataArray( + np.ones((t3_model.sizes['model'], t3_model.sizes['basins'])), + dims=['model', 'basins'], + coords={'model': t3_model.model.values, + 'basins': t3_model.basins.values}) + if t3_models is not None: + t3_weights = t3_weights.where( + t3_weights.model.isin(list(t3_models)), other=0) + + t4_weights = xr.DataArray( + np.ones((targets['t4_mean'].sizes['region'], + targets['t4_mean'].sizes['year'])), + dims=['region', 'year'], + coords={'region': targets['t4_mean'].region.values, + 'year': targets['t4_mean'].year.values}) + if t4_regions is not None: + t4_weights = t4_weights.where( + t4_weights.region.isin(list(t4_regions)), other=0) + if t4_years is not None: + t4_weights = t4_weights.where( + t4_weights.year.isin(list(t4_years)), other=0) + # a year the ensemble did not run cannot contribute, whatever the + # weighting asks for + t4_weights = t4_weights.where( + t4_weights.year.isin(list(units['t4'].year.values)), other=0) + + return dict( + t1_model=t1_model, + t1_obs_mean=targets['t1_mean'], + t1_obs_sigma=targets['t1_sigma'], + t1_weights=t1_weights, + t2_model=t2_model, + t2_obs_mean=targets['t2_mean'], + t2_obs_sigma=targets['t2_sigma'], + t2_weights=targets['t2_weights'], + t3_model=t3_model, + t3_obs_mean=targets['t3_mean'].sel(model=t3_model.model.values), + t3_obs_sigma=targets['t3_sigma'].sel(model=t3_model.model.values), + t3_weights=t3_weights, + t4_model=t4_model, + t4_obs_mean=targets['t4_mean'], + t4_obs_sigma=targets['t4_sigma'], + t4_weights=t4_weights) + + +def run_optimisation(terms, param_values, resolution=8000.0, + sample_size=100000, seed=None): + """ + Sample the objective function and return the parameter distribution. + + Parameters + ---------- + terms : dict + From :py:func:`build_toolbox_terms` + + param_values : numpy.ndarray + The parameter grid the terms were built on + + resolution : float, optional + Passed through to the toolbox; it does not use it for these terms, + which arrive already aggregated + + sample_size : int, optional + Number of random draws of the term weights and the targets + + seed : int, optional + Seed for the random draws. Set it so that the reported percentiles + are reproducible; the toolbox draws from the global numpy random + state. + + Returns + ------- + result : dict + ``p5``, ``median``, ``p95``, ``mode``, the raw ``min_p1``, and + ``at_upper_bound``, the fraction of draws that chose the largest + parameter on the grid + """ + if seed is not None: + np.random.seed(seed) + + min_p1, _ = toolbox.calculate_objective_function( + sample_size, + resolution, + terms['t1_model'], terms['t1_obs_mean'], terms['t1_obs_sigma'], + terms['t1_weights'], + terms['t2_model'], terms['t2_obs_mean'], terms['t2_obs_sigma'], + terms['t2_weights'], + terms['t3_model'], terms['t3_obs_mean'], terms['t3_obs_sigma'], + terms['t3_weights'], + terms['t4_model'], terms['t4_obs_mean'], terms['t4_obs_sigma'], + terms['t4_weights']) + min_p1 = np.asarray(min_p1, dtype=float) + + values = np.asarray(param_values, dtype=float) + step = np.diff(values).min() + edges = np.append(values[0] - 0.5 * step, values + 1.0e-7 * step) + counts, _ = np.histogram(min_p1, bins=edges) + + # a draw that chose the largest parameter on the grid wanted a larger + # one and could not have it, so the grid has clipped the distribution + at_upper_bound = float(np.mean(min_p1 >= values[-1] * (1.0 - 1.0e-9))) + + return dict(min_p1=min_p1, + p5=float(np.percentile(min_p1, 5)), + median=float(np.median(min_p1)), + p95=float(np.percentile(min_p1, 95)), + mode=float(values[int(np.argmax(counts))]), + at_upper_bound=at_upper_bound) diff --git a/compass/landice/tests/ismip7_calibration/quadratic.py b/compass/landice/tests/ismip7_calibration/quadratic.py new file mode 100644 index 0000000000..78ce40d902 --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/quadratic.py @@ -0,0 +1,366 @@ +""" +The ISMIP7 "quadratic local" melt parameterisation, Burgard et al. (2022). + +Protocol Eq. (1): + +.. code-block:: none + + m = K sin(theta) (rho_o/rho_i) (c_o/L_i)^2 beta_S S_loc + (g/(2|f|)) |TF_loc| TF_loc + +This is the reference implementation used both to replicate the published +8 km calibration (see the ``replication`` test case) and as the specification +for MALI's Fortran implementation, selected with +``config_basal_mass_bal_float = 'ismip7'``. + +It is written to agree bit-for-bit with ``multimelt.melt_functions +.quadratic_mixed_slope``, which is what the protocol's own worked example uses, +so that any difference in the calibrated ``K`` is attributable to the mesh and +the ice-sheet model rather than to the melt formula. + +Note on the Coriolis parameter: multimelt uses a *constant* ``f = 1.4e-4`` +rather than a latitude-varying one. The MALI mesh carries ``latCell``, so a +latitude-varying ``f`` would be easy, but it would shift ``K`` away from the +published value. The constant is deliberate, so that ``K`` stays comparable. +""" + +from dataclasses import dataclass + +import numpy as np + +#: seconds per year, as used by multimelt +SECONDS_PER_YEAR = 31556926.080000002 + +#: (rho_sw * c_pw) / (rho_i * L_i), K^-1 -- multimelt's ``melt_factor`` +MELT_FACTOR = 0.013338444158574889 + +#: specific heat capacity of seawater, J kg^-1 K^-1 +C_PO = 3974.0 + +#: latent heat of fusion of ice, J kg^-1 +L_I = 334000.0 + +#: haline contraction coefficient, from Lazeroms et al. +BETA_S = 0.000786 + +#: gravitational acceleration, m s^-2 +GRAVITY = 9.81 + +#: Coriolis parameter, s^-1 (constant, as in multimelt) +F_CORIOLIS = 0.00014 + +#: ice density used by the protocol's worked example, kg m^-3 +ICE_DENSITY = 918.0 + + +@dataclass(frozen=True) +class Constants: + """ + The physical constants entering protocol Eq. (1). + + Two sets matter here. ``MULTIMELT`` reproduces the reference + implementation the published calibration used, and ``MALI`` uses the values + MALI itself compiles in, so that a MALI melt field can be checked against + this formula without the constants confounding the comparison. + + Attributes + ---------- + c_o : float + Specific heat capacity of seawater, J kg-1 K-1. + latent_heat : float + Latent heat of fusion of ice, J kg-1. + gravity : float + Gravitational acceleration, m s-2. + beta_s : float + Haline contraction coefficient, PSU-1. + coriolis : float + Coriolis parameter magnitude, s-1. + rho_ice : float + Ice density appearing in the (rho_o / rho_i) factor of Eq. (1), + kg m-3. + rho_ice_flux : float + Ice density used to convert a melt rate in m of ice to a mass flux, + kg m-3. Normally identical to ``rho_ice``, in which case the two + cancel and the result is independent of ice density. They differ in + the protocol's worked example, which takes 917 in the formula and 918 + in the conversion, leaving a spurious factor of 918/917. + rho_ocean : float + Seawater density, kg m-3. + seconds_per_year : float + Year length used to report the melt rate per year. The melt rate + itself is per second and unambiguous; this only sets the units it is + reported in. MALI uses a 365-day year, matching its noleap calendar, + while the protocol's reference implementation uses 365.2422 days -- a + 0.066% difference, which is enough to fail an exact comparison. + """ + + c_o: float + latent_heat: float + gravity: float + beta_s: float + coriolis: float + rho_ice: float + rho_ice_flux: float + rho_ocean: float + seconds_per_year: float + + +#: constants of the protocol's reference implementation (multimelt) +MULTIMELT = Constants( + c_o=3974.0, + latent_heat=334000.0, + gravity=9.81, + beta_s=0.000786, + coriolis=0.00014, + rho_ice=917.0, + rho_ice_flux=918.0, + rho_ocean=1028.0, + seconds_per_year=SECONDS_PER_YEAR, +) + +#: constants MALI compiles in, from li_constants and the default namelist +MALI = Constants( + c_o=3.974e3, + latent_heat=335.0e3, + gravity=9.80616, + beta_s=7.86e-4, + coriolis=1.4e-4, + rho_ice=910.0, + rho_ice_flux=910.0, + rho_ocean=1028.0, + # li_constants scyr: seconds in a 365-day year + seconds_per_year=31536000.0, +) + + +def u_factor(salinity): + """ + The velocity-scale factor of Jenkins et al. (2018). + + ``(c_o / L_i) * beta_S * g / (2 |f|) * S_loc`` + + Parameters + ---------- + salinity : xarray.DataArray or float + Practical salinity at the ice draft. + + Returns + ------- + Same type as ``salinity``. + """ + return ( + (C_PO / L_I) * BETA_S * (GRAVITY / (2.0 * abs(F_CORIOLIS))) * salinity + ) + + +def local_quadratic_melt( + k, + thermal_forcing, + salinity, + slope, + thermal_forcing_avg=None, + constants=None, + delta_t=0.0, +): + """ + Melt rate from the quadratic local parameterisation, in kg m-2 yr-1. + + Parameters + ---------- + k : float or xarray.DataArray + The calibration parameter ``K``. + thermal_forcing : xarray.DataArray + Local thermal forcing at the ice draft, in K. + salinity : xarray.DataArray + Practical salinity at the ice draft. + slope : float or xarray.DataArray + Ice-draft slope angle in radians (positive). A scalar gives the + "constant Antarctic-mean slope" variant; a field gives the + slope-dependent one. + delta_t : float or xarray.DataArray, optional + Basin-wide thermal-forcing correction, K. Protocol §4.2.1 applies it + wherever the thermal forcing appears, so it is added to both the local + forcing and the averaged one. Defaults to zero, which is what the + calibration uses; production runs are expected to use non-zero values. + thermal_forcing_avg : xarray.DataArray, optional + Thermal forcing to use in the ``|TF|`` factor. Defaults to + ``thermal_forcing``, giving the *local* form; pass a shelf- or + basin-average for the *semi-local* form of protocol Eq. (2). + + Returns + ------- + xarray.DataArray + Melt rate in kg m-2 yr-1, positive for melting. + + Notes + ----- + Melt is exactly linear in ``k``, which is what allows the calibration to + use one model run per ocean state rather than one per (state, parameter) + pair; see the linearity of melt in the parameter. + """ + if thermal_forcing_avg is None: + thermal_forcing_avg = thermal_forcing + thermal_forcing = thermal_forcing + delta_t + thermal_forcing_avg = thermal_forcing_avg + delta_t + + if constants is None: + melt_factor = MELT_FACTOR + u = u_factor(salinity) + rho_ice = ICE_DENSITY + seconds_per_year = SECONDS_PER_YEAR + else: + melt_factor = ( + constants.rho_ocean * constants.c_o / + (constants.rho_ice * constants.latent_heat)) + u = ( + (constants.c_o / constants.latent_heat) * + constants.beta_s * + (constants.gravity / (2.0 * abs(constants.coriolis))) * + salinity) + rho_ice = constants.rho_ice_flux + seconds_per_year = constants.seconds_per_year + + melt_m_per_s = ( + k * melt_factor * u * thermal_forcing * + abs(thermal_forcing_avg) * np.sin(slope)) + return melt_m_per_s * seconds_per_year * rho_ice + + +def angle_from_sin_slope(sin_slope): + """ + The slope angle corresponding to a sine, in radians. + + MALI's ``config_ismip7_melt_sin_slope`` is a *sine*, while + :py:func:`local_quadratic_melt` takes an *angle*, so a call that mixes + the two silently rescales the melt. This converts between them. + + Parameters + ---------- + sin_slope : float + The sine of the ice-draft slope angle + + Returns + ------- + slope : float + The slope angle, in radians + """ + return float(np.arcsin(sin_slope)) + + +def nonlocal_quadratic_melt(gamma0, thermal_forcing, thermal_forcing_mean, + constants=None, delta_t=0.0): + """ + Melt rate from the ISMIP6 non-local quadratic, in kg m-2 yr-1. + + This is the parameterization MALI selects with + ``config_basal_mass_bal_float = 'ismip6'``. With a constant salinity it + is algebraically identical to the Burgard *semi-local* form of protocol + Eq. (2): only the decomposition of the constant differs, so every + ``gamma0`` has an exactly equivalent ``K``. That is why the semi-local + form is calibrated through this path rather than as a separate melt + module. + + Parameters + ---------- + gamma0 : float or xarray.DataArray + The calibration parameter, in m yr-1 + + thermal_forcing : xarray.DataArray + Local thermal forcing at the ice draft, in K + + thermal_forcing_mean : xarray.DataArray + Area-weighted mean thermal forcing over the basin each cell belongs + to, in K + + constants : Constants, optional + The physical constants to use; defaults to :py:data:`MALI` + + delta_t : float or xarray.DataArray, optional + Basin-wide thermal-forcing correction, K. Protocol Sect. 4.2.1 + applies it wherever the thermal forcing appears, so it is added to + both the local forcing and the basin mean. + + Returns + ------- + melt : xarray.DataArray + Melt rate in kg m-2 yr-1, positive for melting + """ + if constants is None: + constants = MALI + thermal_forcing = thermal_forcing + delta_t + thermal_forcing_mean = thermal_forcing_mean + delta_t + + cste = (constants.rho_ocean * constants.c_o / + (constants.rho_ice * constants.latent_heat))**2 + return (gamma0 * cste * constants.rho_ice_flux * + thermal_forcing * abs(thermal_forcing_mean)) + + +def draft_slope(draft, dx, dy, x_dim='x', y_dim='y'): + """ + Ice-draft slope angle on a structured grid, in radians. + + Reproduces the centred-difference scheme of + ``multimelt.plume_functions.check_slope_one_dimension``, including its + one-sided fallbacks at NaN neighbours and its substitution of zero where + the slope cannot be computed at all. + + Parameters + ---------- + draft : xarray.DataArray + Ice draft (negative below sea level), on a regular grid. + dx, dy : float + Grid spacing in m. + x_dim, y_dim : str, optional + Names of the horizontal dimensions. + + Returns + ------- + xarray.DataArray + Slope angle in radians. + """ + slope_x = _one_dimensional_slope(draft, x_dim, dx) + slope_y = _one_dimensional_slope(draft, y_dim, dy) + return np.arctan(np.sqrt(slope_x**2 + slope_y**2)) + + +def _one_dimensional_slope(draft, dim, spacing): + """One-sided-tolerant centred difference, as in multimelt.""" + shifted_minus = draft.shift({dim: -1}) + shifted_plus = draft.shift({dim: 1}) + + both = (shifted_minus - shifted_plus) / np.sqrt((2.0 * spacing) ** 2) + right = (draft - shifted_plus) / np.sqrt(spacing**2) + left = (shifted_minus - draft) / np.sqrt(spacing**2) + + slope = both.combine_first(right).combine_first(left) + return slope.where(np.isfinite(slope), 0.0) + + +def mean_slope(draft, floating, dx, dy, x_dim='x', y_dim='y'): + """ + Antarctic-mean ice-draft slope angle over floating ice, in radians. + + This is the "constant slope" used by the protocol's worked example. Note + that it is a mean of *angles*, not of ``sin(theta)``, and that its value is + resolution dependent; both points are raised in + ``protocol-and-toolbox-questions.md`` A4. + + Parameters + ---------- + draft : xarray.DataArray + Ice draft on a regular grid. + floating : xarray.DataArray + Boolean mask, True on floating ice. + dx, dy : float + Grid spacing in m. + x_dim, y_dim : str, optional + Names of the horizontal dimensions. + + Returns + ------- + float + Mean slope angle in radians. + """ + slope = draft_slope(draft, dx, dy, x_dim=x_dim, y_dim=y_dim) + return float(slope.where(floating).mean()) diff --git a/compass/landice/tests/ismip7_calibration/replication/__init__.py b/compass/landice/tests/ismip7_calibration/replication/__init__.py new file mode 100644 index 0000000000..7e30bba1fe --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/replication/__init__.py @@ -0,0 +1,60 @@ +from compass.landice.tests.ismip7_calibration.configure import check_options +from compass.landice.tests.ismip7_calibration.replication.replicate import ( + Replicate, + check_published, +) +from compass.testcase import TestCase +from compass.validate import compare_variables + + +class Replication(TestCase): + """ + A test case that reproduces the published ISMIP7 quadratic calibration on + the ISMIP 8 km grid, through compass's own code path. + + The protocol reports, for the quadratic local parameterization on the + ISMIP 8 km grid (Sect. 4.3.1 and Fig. 5), 5th / 50th / 95th percentiles + of K = 4.75e-5 / 8.5e-5 / 1.375e-4. Reproducing those with compass's + area-weighted, mesh-agnostic terms driving the vendored objective + function validates the whole of the parameter-selection stage before any + MALI run exists, and is the regression test that the unstructured + generalization must not break. + + This test case deliberately uses the **published inputs**, including the + spatially varying 8 km salinity fields. It replicates a published + calculation, so it must not be changed to match the constant-salinity + choice the MALI calibration makes; it validates the parameter-selection + machinery, not MALI's physics. + """ + + 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 = 'replication' + super().__init__(test_group=test_group, name=name, subdir=name) + + self.add_step(Replicate(test_case=self)) + + def configure(self): + """ + Check that the ISMIP7 dataset path has been supplied + """ + check_options(self.config, ['base_path_ismip7']) + + def validate(self): + """ + Check the percentiles against the published values, and against a + baseline if one was provided + """ + filename = 'replicate/replication_8km.nc' + check_published(filename, self.logger) + compare_variables(test_case=self, + variables=['p5', 'median', 'p95', 'mode', + 'min_p1'], + filename1=filename) diff --git a/compass/landice/tests/ismip7_calibration/replication/replicate.py b/compass/landice/tests/ismip7_calibration/replication/replicate.py new file mode 100644 index 0000000000..dbd999567e --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/replication/replicate.py @@ -0,0 +1,278 @@ +""" +Reproduce the published 8 km ISMIP7 quadratic calibration. +""" + +import os + +import numpy as np +import xarray as xr +from mpas_tools.io import write_netcdf + +from compass.landice.tests.ismip7_calibration import datasets +from compass.landice.tests.ismip7_calibration.configure import ( + objective_options, +) +from compass.landice.tests.ismip7_calibration.objective import ( + build_toolbox_terms, + run_optimisation, +) +from compass.landice.tests.ismip7_calibration.quadratic import ( + local_quadratic_melt, + mean_slope, +) +from compass.landice.tests.ismip7_calibration.terms import ( + average_by_group, + integrate_by_group, + uniform_area, +) +from compass.step import Step + +#: the ISMIP grid resolution the published calibration used, m +RESOLUTION = 8000.0 + +#: cell dimensions of the ISMIP structured grid +CELL_DIMS = ('y', 'x') + +#: the 120 K values the protocol's worked example samples +K_VALUES = np.arange(0.25e-5, 3.025e-4, 0.25e-5) + +#: percentiles published in protocol Sect. 4.3.1 and Fig. 5 +PUBLISHED = {'p5': 4.75e-5, 'median': 8.5e-5, 'p95': 1.375e-4} + + +class Replicate(Step): + """ + A step that reproduces the published 8 km quadratic calibration. + + Attributes + ---------- + base_path : str + Root of the ISMIP7 AIS datasets + """ + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_calibration.replication.Replication + The test case this step belongs to + """ # noqa: E501 + super().__init__(test_case=test_case, name='replicate') + self.base_path = None + self.add_output_file(filename='replication_8km.nc') + + def setup(self): + """ + Check that the ISMIP7 datasets this step needs are present + """ + config = self.config + self.base_path = config.get('ismip7_calibration', 'base_path_ismip7') + + states = datasets.ocean_states(self.base_path, subset='all') + missing = datasets.missing_files(states) + if missing: + listing = '\n '.join(f'{name}: {path}' for name, path in missing) + raise FileNotFoundError( + f'{len(missing)} ISMIP7 input files are missing:\n ' + f'{listing}') + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + base_path = config.get('ismip7_calibration', 'base_path_ismip7') + + static = _load_static(base_path, logger) + units = _build_unit_terms(base_path, static, logger) + targets = datasets.load_targets(base_path) + + # The published example weights four ocean models, and PIG alone in + # 2009 and 2012. Reproducing the published numbers means reproducing + # that weighting, so it is fixed here rather than taken from config. + terms = build_toolbox_terms( + units, targets, K_VALUES, + t3_models=datasets.PUBLISHED_T3_MODELS, + t4_regions=datasets.PUBLISHED_T4_REGIONS, + t4_years=datasets.PUBLISHED_T4_YEARS) + + sample_size, seed = objective_options(config) + logger.info(f'Sampling the objective function {sample_size} times...') + result = run_optimisation(terms, K_VALUES, resolution=RESOLUTION, + sample_size=sample_size, seed=seed) + + _write_result(result, K_VALUES, 'replication_8km.nc') + _report(result, logger) + + +def check_published(filename, logger): + """ + Check the replicated percentiles against the published values. + + Parameters + ---------- + filename : str + The file written by :py:class:`Replicate` + + logger : logging.Logger + A logger to report the comparison to + + Raises + ------ + ValueError + If any percentile differs from the published value + """ + result = xr.open_dataset(filename) + mismatched = [] + for name, expected in PUBLISHED.items(): + found = float(result[name]) + # the published values are points on the parameter grid, so this + # should be exact to within floating-point representation rather + # than merely close + if not np.isclose(found, expected, rtol=1.0e-9, atol=0.0): + mismatched.append(f' {name}: expected {expected:.5e}, ' + f'found {found:.5e}') + result.close() + + if mismatched: + listing = '\n'.join(mismatched) + raise ValueError( + f'The replication of the published 8 km ISMIP7 calibration does ' + f'not reproduce the published percentiles:\n{listing}\n' + f'This means the vendored parameter-selection toolbox is not ' + f'being driven as it was for the published numbers, so the MALI ' + f'calibration built on the same code path cannot be trusted ' + f'either.') + + logger.info('The published 8 km percentiles are reproduced exactly.') + + +def _load_static(base_path, logger): + """Load the 8 km topography, masks and derived fields.""" + logger.info('Loading the ISMIP 8 km topography and masks...') + masks = datasets.mask_files(base_path) + + topo = xr.load_dataset( + os.path.join(base_path, 'obs', 'ocean', 'topography', 'bedmap3', 'v3', + 'bedmap3_AIS_obs_ocean_topography_v3.nc')) + basins = xr.load_dataset(masks['basins']).basinNumber.rename('basins') + bfrn = xr.load_dataset(masks['bfrn']) + mask = xr.load_dataset(masks['floating']).mask + + floating = topo['floating_frac'] > 0.5 + slope = mean_slope(topo['draft'], floating, RESOLUTION, RESOLUTION) + logger.info(f' mean draft slope: {slope:.6f} rad ' + f'(sin = {np.sin(slope):.7f})') + + shelves = xr.load_dataset(masks['shelves']).shelf_mask.isel(time=0) + x = shelves['x'] if 'x' in shelves.coords else basins['x'] + # restrict Pine Island to its main trunk, as the worked example does + pig = (shelves == datasets.PIG_ID) & (x > datasets.PIG_X_MAX) + dotson = shelves == datasets.DOTSON_ID + region_label = xr.where(pig, 'pig', xr.where(dotson, 'dotson', '')) + + return dict(draft=topo['draft'], floating=floating, basins=basins, + bfrn=bfrn, mask=mask, slope=slope, region_label=region_label, + area=uniform_area(topo['draft'], RESOLUTION, CELL_DIMS)) + + +def _unit_melt(state, static): + """ + Melt at a unit ``K`` for one ocean state, in kg m-2 yr-1. + + Melt is exactly linear in ``K``, so every member of the parameter + ensemble is this field times its ``K``. Working at ``K = 1`` keeps the + whole replication in single fields rather than 120 copies of each. + """ + draft = static['draft'] + floating = static['floating'] + + tf = xr.load_dataset(state.tf_file)['tf'] + so = xr.load_dataset(state.so_file)['so'] + + tf_draft = tf.sel(z=draft, method='nearest').where(floating) + so_draft = so.sel(z=draft, method='nearest').where(floating) + + melt = local_quadratic_melt(1.0, tf_draft, so_draft, static['slope']) + # the observational files carry a scalar Time coordinate that differs + # between years, which would collide when the per-year aggregates are + # concatenated; nothing downstream uses any non-dimension coordinate + return melt.drop_vars([name for name in melt.coords + if name not in melt.dims], errors='ignore') + + +def _build_unit_terms(base_path, static, logger): + """Aggregate unit melt into the four objective-function terms.""" + area, mask = static['area'], static['mask'] + basins, bfrn = static['basins'], static['bfrn'] + states = datasets.ocean_states(base_path, subset='all') + by_name = {state.name: state for state in states} + + logger.info('Building the present-day ensemble (J1, J2)...') + pd_unit = _unit_melt(by_name['climatology'], static) + t1 = integrate_by_group(pd_unit, area, mask, basins, CELL_DIMS, + group_dim='basins') + t2 = integrate_by_group(pd_unit, area, mask, bfrn['BFRN_bins'], + CELL_DIMS, group_dim='BFRN_bins') + + logger.info('Building the ocean-model ensemble (J3)...') + means = {} + for which in ('cold', 'warm'): + per_model = [] + for state in states: + if state.kind != 'model' or state.state != which: + continue + melt = _unit_melt(state, static) + # regional models only constrain the basins they cover + if state.basins is not None: + melt = melt.where(basins.isin(list(state.basins))) + agg = average_by_group(melt, area, mask, basins, CELL_DIMS, + group_dim='basins') + per_model.append(agg.expand_dims({'model': [state.label]})) + logger.info(f' {which:4s} {state.label}') + means[which] = xr.concat(per_model, dim='model', + coords='minimal') + t3 = means['warm'] - means['cold'] + + logger.info('Building the observational ensemble (J4)...') + per_year = [] + for state in states: + if state.kind != 'obs': + continue + melt = _unit_melt(state, static) + agg = integrate_by_group(melt, area, mask, static['region_label'], + CELL_DIMS, group_dim='region') + per_year.append(agg.expand_dims({'year': [state.year]})) + logger.info(f' {state.year}') + t4 = xr.concat(per_year, dim='year', coords='minimal') + # drop the '' label used for cells outside PIG and Dotson + t4 = t4.sel(region=[region for region in t4.region.values if region]) + + return dict(t1=t1, t2=t2, t3=t3, t4=t4) + + +def _write_result(result, param_values, filename): + """Write the parameter distribution to a file.""" + ds = xr.Dataset() + ds['min_p1'] = ('sample', result['min_p1']) + ds['parameter_values'] = ('parameter', np.asarray(param_values)) + for name in ('p5', 'median', 'p95', 'mode'): + ds[name] = float(result[name]) + ds.attrs['description'] = ( + 'Replication of the published ISMIP7 quadratic calibration on the ' + 'ISMIP 8 km grid') + write_netcdf(ds, filename) + + +def _report(result, logger): + """Log the percentiles beside the published values.""" + logger.info('') + logger.info('K percentiles, ISMIP 8 km grid:') + logger.info(f'{"":12s}{"compass":>12s}{"published":>12s}') + for name, expected in PUBLISHED.items(): + logger.info(f'{name:12s}{result[name]:12.5e}{expected:12.5e}') + logger.info(f'{"mode":12s}{result["mode"]:12.5e}') + logger.info('') diff --git a/compass/landice/tests/ismip7_calibration/terms.py b/compass/landice/tests/ismip7_calibration/terms.py new file mode 100644 index 0000000000..3696764259 --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/terms.py @@ -0,0 +1,304 @@ +""" +Objective-function terms J1-J4 for ISMIP7 melt-module calibration, on any mesh. + +The reference implementation in +``parameterisations/parameter_selection_toolbox.py`` assumes a uniform +structured grid: it converts cell sums to Gt/yr with a single scalar +``cvt_m = reso**2 / 1e12`` and takes unweighted means over basins. That is +correct for the ISMIP 8 km grid but wrong for a variable-resolution +unstructured mesh such as MALI's 4-20 km Antarctic mesh, where cell area +varies by a factor of ~25. + +The functions here take an explicit cell-area array instead, and work over an +arbitrary set of cell dimensions. Passing a uniform area reproduces the +structured-grid result exactly, so the same code serves both cases; see +``tests/test_terms.py``. + +Conventions, matching the toolbox: + +* melt rates are in kg m-2 yr-1, positive for melting, negative for refreezing +* integrated terms (J1, J2, J4) are in Gt yr-1 +* averaged terms (J3) are in kg m-2 yr-1 +* an aggregate that comes out exactly zero is set to NaN, so that regions + with no contributing cells drop out of the objective function rather than + being treated as a real zero +""" + +import numpy as np +import xarray as xr + +#: kg per Gt +KG_PER_GT = 1.0e12 + +#: name of the flattened cell dimension used internally +CELL_DIM = '_cell' + + +def stack_cells(da, cell_dims): + """ + Flatten the cell dimensions of ``da`` into a single dimension. + + Parameters + ---------- + da : xarray.DataArray + Array with one or more cell dimensions. + cell_dims : sequence of str + The dimensions to flatten, e.g. ``('y', 'x')`` for a structured grid or + ``('nCells',)`` for an MPAS mesh. + + Returns + ------- + xarray.DataArray + ``da`` with ``cell_dims`` replaced by :data:`CELL_DIM`. + """ + cell_dims = tuple(cell_dims) + missing = [dim for dim in cell_dims if dim not in da.dims] + if missing: + raise ValueError(f'{missing} not found in dimensions {tuple(da.dims)}') + if len(cell_dims) == 1 and cell_dims[0] == CELL_DIM: + return da + stacked = da.stack({CELL_DIM: cell_dims}) + # drop the MultiIndex coordinate; we only need positional alignment and the + # index would otherwise have to be carried through every groupby + return stacked.drop_vars( + [CELL_DIM, *cell_dims], errors='ignore' + ).assign_coords({CELL_DIM: np.arange(stacked.sizes[CELL_DIM])}) + + +def _prepare(melt, area, mask, groups, cell_dims): + """Flatten melt, area, mask and groups onto a common cell dimension.""" + melt = stack_cells(melt, cell_dims) + area = stack_cells(area, cell_dims) + groups = stack_cells(groups, cell_dims) + if mask is None: + valid = xr.ones_like(area, dtype=bool) + else: + valid = stack_cells(mask, cell_dims).astype(bool) + # cells with no group assignment never contribute + if groups.dtype.kind == 'f': + valid = valid & groups.notnull() + return melt, area, valid, groups + + +def integrate_by_group( + melt, area, mask, groups, cell_dims, group_dim='basins' +): + """ + Integrate a melt rate over each group, in Gt yr-1. + + This is the mesh-agnostic form of the toolbox's + ``melt.where(mask).groupby(groups).sum() * reso**2 / 1e12``. + + Parameters + ---------- + melt : xarray.DataArray + Melt rate in kg m-2 yr-1. May carry any number of extra dimensions + (parameter values, ocean model, year, ...); they are preserved. + area : xarray.DataArray + Cell area in m2, over ``cell_dims`` only. + mask : xarray.DataArray or None + Boolean; True where the cell should contribute (e.g. floating ice). + groups : xarray.DataArray + Group id per cell (basin number, buttressing bin, ...). + cell_dims : sequence of str + Dimensions of ``melt`` that index cells. + group_dim : str, optional + Name to give the resulting group dimension. + + Returns + ------- + xarray.DataArray + Integrated melt in Gt yr-1, with ``cell_dims`` replaced by + ``group_dim``. Groups whose total is exactly zero are set to NaN. + """ + melt, area, valid, groups = _prepare(melt, area, mask, groups, cell_dims) + + weighted = (melt * area).where(valid) + total = weighted.groupby(groups.rename(group_dim)).sum(skipna=True) + total = total / KG_PER_GT + return total.where(total != 0.0) + + +def average_by_group(melt, area, mask, groups, cell_dims, group_dim='basins'): + """ + Area-weighted mean melt rate over each group, in kg m-2 yr-1. + + This is the mesh-agnostic form of the toolbox's + ``melt.where(mask).groupby(groups).mean()``. On a uniform grid an + area-weighted mean and a plain mean coincide; on a variable-resolution mesh + they do not, and the area-weighted one is the physically meaningful choice. + + Parameters and returns are as for :func:`integrate_by_group`, except + that the result is a mean rather than an integral. + """ + melt, area, valid, groups = _prepare(melt, area, mask, groups, cell_dims) + + grouper = groups.rename(group_dim) + weighted = (melt * area).where(valid) + # only count area where melt is defined, so that NaN melt does not bias the + # denominator + weights = area.where(valid & melt.notnull()) + + numer = weighted.groupby(grouper).sum(skipna=True) + denom = weights.groupby(grouper).sum(skipna=True) + mean = numer / denom.where(denom != 0.0) + return mean.where(mean != 0.0) + + +def calculate_term1( + ensemble, area, mask, basins, melt_obs, cell_dims, var='melt_rate' +): + """ + J1: basin-integrated present-day melt, in Gt yr-1. + + Parameters + ---------- + ensemble : xarray.Dataset + Present-day ensemble, with ``var`` indexed by at least ``p1`` and + ``p2``. + area, mask, basins : xarray.DataArray + Cell area (m2), contributing-cell mask, and IMBIE2 basin number per + cell. Basin numbers follow the **ISMIP7 0-based convention**, + which differs by one from MALI's; see the basin-numbering note in + the developer guide for why that matters. + melt_obs : pandas.DataFrame + Observational targets, with columns ``'BMR (Gt/yr)'`` and + ``'BMR uncert (Gt/yr)'``, one row per basin in basin order. + cell_dims : sequence of str + Cell dimensions of the ensemble. + var : str, optional + Name of the melt-rate variable. + + Returns + ------- + model, obs_mean, obs_sigma : xarray.DataArray + """ + model = integrate_by_group( + ensemble[var], area, mask, basins, cell_dims, group_dim='basins' + ) + obs_mean = xr.DataArray( + data=np.asarray(melt_obs['BMR (Gt/yr)'].values, dtype=float), + name='melt_Gt_per_y', + dims=['basin'], + coords={'basin': np.arange(len(melt_obs))}, + ) + obs_sigma = xr.DataArray( + data=np.asarray(melt_obs['BMR uncert (Gt/yr)'].values, dtype=float), + name='melt_unc_Gt_per_y', + dims=['basin'], + coords={'basin': np.arange(len(melt_obs))}, + ) + return model, obs_mean, obs_sigma + + +def calculate_term2( + ensemble, area, mask, bfrn_bins, target, cell_dims, var='melt_rate' +): + """ + J2: melt integrated over buttressing (BFRN) bins, in Gt yr-1. + + ``target`` is the dataset holding ``melt_mean`` and ``melt_mean_err`` per + bin. Other parameters are as for :func:`calculate_term1`. + """ + model = integrate_by_group( + ensemble[var], area, mask, bfrn_bins, cell_dims, group_dim='BFRN_bins' + ) + return model, target['melt_mean'], target['melt_mean_err'] + + +def calculate_term3( + cold_ensemble, + warm_ensemble, + cold_target, + warm_target, + area, + mask, + basins, + cell_dims, + var='melt_rate', +): + """ + J3: warm-minus-cold basin-mean melt difference, in kg m-2 yr-1. + + The model term is the difference of area-weighted basin means between the + warm and cold ocean states; the target is the corresponding difference from + the ocean models, with uncertainties combined in quadrature. + """ + cold = average_by_group( + cold_ensemble[var], area, mask, basins, cell_dims, group_dim='basins' + ) + warm = average_by_group( + warm_ensemble[var], area, mask, basins, cell_dims, group_dim='basins' + ) + model = warm - cold + + obs_mean = warm_target.melt_rate - cold_target.melt_rate + obs_sigma = np.sqrt( + warm_target.melt_rate_uncert**2 + cold_target.melt_rate_uncert**2 + ) + return model, obs_mean, obs_sigma + + +def calculate_term4( + obs_ensemble, + area, + mask, + region_label, + target, + cell_dims, + var='melt_rate', +): + """ + J4: PIG/Dotson integrated melt per observation year, in Gt yr-1. + + ``region_label`` labels each cell with an ice-shelf name (``'pig'``, + ``'dotson'``); ``target`` holds ``melt_rate`` and ``melt_rate_uncert`` + indexed by region and year. The result is reindexed onto the target so the + region ordering matches. + """ + model = integrate_by_group( + obs_ensemble[var], + area, + mask, + region_label, + cell_dims, + group_dim='region', + ) + model = model.where(target.melt_rate.notnull()) + model = model.reindex_like(target.melt_rate) + return model, target.melt_rate, target.melt_rate_uncert + + +def uniform_area(template, resolution, cell_dims): + """ + Build a uniform cell-area array, for structured grids. + + Convenience so that a structured-grid case can call the same functions: + ``area = uniform_area(melt, 8000.0, ('y', 'x'))``. + + Parameters + ---------- + template : xarray.DataArray + Any array carrying the cell dimensions and their coordinates. + resolution : float + Grid spacing in m; cell area is ``resolution**2``. + cell_dims : sequence of str + The cell dimensions. + + Returns + ------- + xarray.DataArray + Constant array of ``resolution**2``, over ``cell_dims`` only. + """ + coords = { + dim: template.coords[dim] + for dim in cell_dims + if dim in template.coords + } + shape = tuple(template.sizes[dim] for dim in cell_dims) + return xr.DataArray( + np.full(shape, float(resolution) ** 2), + dims=tuple(cell_dims), + coords=coords, + name='area', + ) diff --git a/compass/landice/tests/ismip7_calibration/toolbox/PROVENANCE.md b/compass/landice/tests/ismip7_calibration/toolbox/PROVENANCE.md new file mode 100644 index 0000000000..7760bdeb83 --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/toolbox/PROVENANCE.md @@ -0,0 +1,47 @@ +# Vendored `parameter_selection_toolbox.py` + +`parameter_selection_toolbox.py` in this directory is an **unmodified, verbatim +copy** of the ISMIP7 parameter-selection toolbox. It implements the objective +function of the ISMIP7 AIS ice-ocean protocol (Reese et al., Sect. 4.2) and the +downstream `deltaT` fit. + +## Where it came from + +| | | +|---|---| +| Repository | | +| Path | `parameterisations/parameter_selection_toolbox.py` | +| Last commit to touch the file | `132beb155e63fd07957e28407a83cb9d2c954447` ("Clean up toolbox code", 2026-06-23) | +| Repository `main` when copied | `84d17692935867d62b9c8e26f076a9869d74be31` (2026-09-08) | +| SHA256 | `3d016c04987d7c66c2603e2c3110d8a6a68cee0b6193c7f47bccf64aaa1ce029` | + +**The upstream tags do not track the toolbox.** The newest toolbox tag, +`param-toolbox-v1`, points at a commit from 2026-04-20, which is *older* than +the last change to this file. Pinning to that tag would pin the wrong code, so +the commit is recorded instead. If upstream starts tagging toolbox updates +reliably, record the tag here as well. + +## Why it is vendored rather than depended upon + +There is no conda-forge or PyPI package that contains this module. It lives at +the top level of `parameterisations/`, outside the `i7aof` package that the +upstream repository does distribute, so even installing that package would not +provide it. The alternatives considered were a git submodule for a single file, +and reimplementing the objective function in Compass. Vendoring keeps exact +provenance without either cost. + +## The integrity check + +`__init__.py` verifies the SHA256 above on import, and +`compass/landice/tests/ismip7_calibration/tests/test_toolbox.py` asserts it in +CI. A local edit or a partial update therefore fails loudly, rather than +silently changing published calibration numbers. + +## How to update + +1. Copy the new upstream file over this one, unmodified. +2. Update the commit, `main` and SHA256 rows above (`sha256sum` the file). +3. Update `_EXPECTED_SHA256` in `__init__.py`. +4. Run the `replication` test case. It must still reproduce the published 8 km + percentiles 4.75e-5 / 8.5e-5 / 1.375e-4 exactly. If it does not, the change + is not a refactor and the calibration results need revisiting. diff --git a/compass/landice/tests/ismip7_calibration/toolbox/__init__.py b/compass/landice/tests/ismip7_calibration/toolbox/__init__.py new file mode 100644 index 0000000000..a8e7e82b8d --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/toolbox/__init__.py @@ -0,0 +1,94 @@ +""" +The ISMIP7 parameter-selection toolbox, vendored from upstream. + +:py:mod:`~compass.landice.tests.ismip7_calibration.toolbox.parameter_selection_toolbox` +is a verbatim copy of ``parameterisations/parameter_selection_toolbox.py`` from +https://github.com/ismip/ismip7-antarctic-ocean-forcing. It implements the +objective function of the ISMIP7 AIS ice-ocean protocol (Reese et al., +Sect. 4.2) and the downstream ``deltaT`` fit. ``PROVENANCE.md`` in this +directory records exactly which upstream revision it came from and how to +update it. + +The copy must stay byte-for-byte identical to upstream, so that a future +refresh is a clean replace and so that the published calibration numbers +cannot change without anyone noticing. :py:func:`check_integrity` enforces +that, and is called on import. +""" # noqa: E501 + +import hashlib +import os + +#: SHA256 of the vendored toolbox; see ``PROVENANCE.md`` +EXPECTED_SHA256 = \ + '3d016c04987d7c66c2603e2c3110d8a6a68cee0b6193c7f47bccf64aaa1ce029' + +#: upstream commit that last modified the vendored file +UPSTREAM_COMMIT = '132beb155e63fd07957e28407a83cb9d2c954447' + +#: upstream repository the file was copied from +UPSTREAM_URL = 'https://github.com/ismip/ismip7-antarctic-ocean-forcing' + + +def toolbox_path(): + """ + The path to the vendored toolbox source file. + + Returns + ------- + path : str + Absolute path to ``parameter_selection_toolbox.py`` + """ + return os.path.join(os.path.dirname(os.path.abspath(__file__)), + 'parameter_selection_toolbox.py') + + +def file_sha256(): + """ + The SHA256 of the vendored toolbox source file, as it is on disk. + + Returns + ------- + digest : str + Hex digest of ``parameter_selection_toolbox.py`` + """ + with open(toolbox_path(), 'rb') as handle: + return hashlib.sha256(handle.read()).hexdigest() + + +def check_integrity(): + """ + Verify that the vendored toolbox still matches the recorded upstream copy. + + Raises + ------ + RuntimeError + If the file on disk does not match :py:data:`EXPECTED_SHA256`. The + vendored copy must be an unmodified upstream file, because the + published ISMIP7 calibration numbers depend on it; a local edit or a + partial update should fail loudly rather than quietly change results. + """ + digest = file_sha256() + if digest != EXPECTED_SHA256: + raise RuntimeError( + f'The vendored ISMIP7 parameter-selection toolbox at\n' + f' {toolbox_path()}\n' + f'does not match the recorded upstream copy.\n' + f' expected SHA256: {EXPECTED_SHA256}\n' + f' found SHA256: {digest}\n' + f'The file must be a verbatim copy of\n' + f' {UPSTREAM_URL}\n' + f' parameterisations/parameter_selection_toolbox.py\n' + f'at commit {UPSTREAM_COMMIT}. If you are deliberately updating ' + f'it, follow the instructions in PROVENANCE.md in this ' + f'directory.') + + +check_integrity() + +from compass.landice.tests.ismip7_calibration.toolbox import ( # noqa: E402 + parameter_selection_toolbox, +) + +__all__ = ['parameter_selection_toolbox', 'check_integrity', 'file_sha256', + 'toolbox_path', 'EXPECTED_SHA256', 'UPSTREAM_COMMIT', + 'UPSTREAM_URL'] diff --git a/compass/landice/tests/ismip7_calibration/toolbox/parameter_selection_toolbox.py b/compass/landice/tests/ismip7_calibration/toolbox/parameter_selection_toolbox.py new file mode 100644 index 0000000000..053d53069e --- /dev/null +++ b/compass/landice/tests/ismip7_calibration/toolbox/parameter_selection_toolbox.py @@ -0,0 +1,701 @@ +import os + +import numpy as np +import xarray as xr + + +def calculate_objective_function( + sample_size, + reso, + t1_model, + t1_obs_mean, + t1_obs_sigma, + t1_weights, + t2_model, + t2_obs_mean, + t2_obs_sigma, + t2_weights, + t3_model, + t3_obs_mean, + t3_obs_sigma, + t3_weights, + t4_model, + t4_obs_mean, + t4_obs_sigma, + t4_weights, +): + """ + Input: + To create input datasets, please see calc_term1, etc + - sample_size: number of samples to take + - reso: native model resolution + - t1_model: xarray dataset containing modelled melt rates, aggregated + and indexed by basins + - t1_obs_mean, sigma: target melt rates in basins, and uncertainties + - t2_model: xarray dataset containing modelled melt rates, aggregated + and indexed by buttressing bins + - t2_obs_mean, sigma: target melt rates in buttressing bins, + and uncertainties + - t2_weights: based on buttressing values + - t3_model: xarray dataset containing modelled melt rates, + averaged and indexed by basins, and ocean models + - t3_obs_mean, sigma: target melt rates in basins and ocean models, + and uncertainties + - t3_weights: weighting of different terms + - t4_model: xarray dataset containing modelled melt rates, + aggregated and indexed by region = pine island, dotson, + and observations year + - t4_obs_mean, sigma: target melt rates in regions, years, + and uncertainties + - t4_weights: weighting of different terms + Output: + -min_p1: list of p1 values that minimise randomly sampled objective + function, length of sample size + -min_p2: list of cooresponding p2 values + + This function finds the pair of p1,p2 for which the parameterised melt + optimises three terms: + -J1: basin-integrated melt for present-day + -J2: buttressing-bin integrated melt for present-day, + weighted by buttressing + -J3: basin-averaged melt for cold and warm cases of the + ocean modelling datasets + -J4: basin-integrated melt from observations of PIG and Dotson + """ + + nBasins = int(t1_model.basins.values.max()) + + ################################################ + # randomly sample the weights of the terms + a1 = xr.DataArray( + np.random.uniform(0, 1, size=sample_size), dims=['sample'] + ) + a2 = xr.DataArray( + np.random.uniform(0, 1, size=sample_size), dims=['sample'] + ) + a3 = xr.DataArray( + np.random.uniform(0, 1, size=sample_size), dims=['sample'] + ) + a4 = xr.DataArray( + np.random.uniform(0, 1, size=sample_size), dims=['sample'] + ) + + ############################################### + # Sample uncertainties in terms + + # Sample uncertainty in term 1 + t1_target_s = [] + for b in range(nBasins + 1): + t1_target_s = t1_target_s + [ + np.random.normal( + loc=t1_obs_mean.values[b], + scale=t1_obs_sigma.values[b], + size=sample_size, + ) + ] + t1_target_s = xr.DataArray( + t1_target_s, + dims=['basins', 'sample'], + coords={ + 'basins': t1_model.basins.values, + 'sample': np.arange(sample_size), + }, + ) + + # Sample uncertainty in term 2 + t2_target_s = [] + nBins = 10 + for b in range(nBins): + t2_target_s = t2_target_s + [ + np.random.normal( + loc=t2_obs_mean[b], + scale=t2_obs_sigma[b], + size=sample_size, + ) + ] + t2_target_s = xr.DataArray( + t2_target_s, + dims=['BFRN_bins', 'sample'], + coords={ + 'BFRN_bins': t2_model.BFRN_bins.values, + 'sample': np.arange(sample_size), + }, + ) + + # Term 3 targets + samples = np.random.normal( + loc=t3_obs_mean.values[..., np.newaxis], + scale=t3_obs_sigma.values[..., np.newaxis], + size=(*t3_obs_mean.values.shape, sample_size), + ) + t3_target_s = xr.DataArray( + samples, + dims=[*t3_obs_mean.dims, 'sample'], + coords={**t3_obs_mean.coords, 'sample': np.arange(sample_size)}, + ) + + t3_weights_samples = xr.DataArray( + np.random.uniform( + 0, + 1, + size=( + len(t3_obs_mean.model.values), + len(t3_obs_mean.basins.values), + sample_size, + ), + ), + dims=['model', 'basins', 'sample'], + coords={ + 'model': t3_obs_mean.model.values, + 'basins': t3_obs_mean.basins.values, + 'sample': np.arange(sample_size), + }, + ) + + # Term 4 + samples = np.random.normal( + loc=t4_obs_mean.values[..., np.newaxis], + scale=t4_obs_sigma.values[..., np.newaxis], + size=(*t4_obs_mean.values.shape, sample_size), + ) + t4_target_s = xr.DataArray( + samples, + dims=[*t4_obs_mean.dims, 'sample'], + coords={**t4_obs_mean.coords, 'sample': np.arange(sample_size)}, + ) + t4_weights_samples = xr.DataArray( + np.random.uniform( + 0, + 1, + size=( + len(t4_obs_mean.region.values), + len(t4_obs_mean.year.values), + sample_size, + ), + ), + dims=['region', 'year', 'sample'], + coords={ + 'region': t4_obs_mean.region.values, + 'year': t4_obs_mean.year.values, + 'sample': np.arange(sample_size), + }, + ) + + ################### + # Calculate objective function + + term1 = mae(t1_model, t1_target_s, t1_weights, 'basins') + term2 = mae(t2_model, t2_target_s, t2_weights, ['BFRN_bins']) + # important to use skipna here + term3 = mae( + t3_model, + t3_target_s, + t3_weights * t3_weights_samples, + ['basins', 'model'], + True, + ) + # important to use skipna here + term4 = mae( + t4_model, + t4_target_s, + t4_weights * t4_weights_samples, + ['region', 'year'], + True, + ) + + # Sample the objective function + eps = 0.000001 # to avoid divison by 0 + objective_function = ( + a1 * term1 / (term1.median(dim=['p1', 'p2']) + eps) + + a2 * term2 / (term2.median(dim=['p1', 'p2']) + eps) + + a3 * term3 / (term3.median(dim=['p1', 'p2']) + eps) + + a4 * term4 / (term4.median(dim=['p1', 'p2']) + eps) + ) + + # Identify minimum values + objective_function_stacked = objective_function.stack(params=('p1', 'p2')) + min_params = objective_function_stacked.idxmin(dim='params').values + + min_p1 = np.array([v[0] for v in min_params]) + min_p2 = np.array([v[1] for v in min_params]) + + return min_p1, min_p2 + + +def mae(predicted=None, observed=None, weights=1, dims='basins', skipna=False): + """ + Calculates mean absolute error + """ + return ( + abs(weights * (predicted - observed)) + .mean(dims, skipna=skipna) + .rename('result') + ) + + +def calculate_term1(pd_ensemble, mask_m, basins_m, nBasins, cvt_m, MeltData): + ######## + # TERM 1 + # parameterisaition melt, aggregate to Gt/a per basin + t1_model = ( + pd_ensemble['melt_rate'] + .where(mask_m, np.nan) + .groupby(basins_m) + .sum(skipna=True) + * cvt_m + ) # convert to Gt/a + # make sure to remove regions that do not have optimal dT for any basin + t1_model = t1_model.where(t1_model != 0, np.nan) + + # Observed melt in Gt/a per basin, observed melt is "sample_size"-times + # randomly sampled assuming normal distribution + t1_obs_mean = MeltData['BMR (Gt/yr)'].values + t1_obs_mean = xr.DataArray( + data=t1_obs_mean, + name='melt_Gt_per_y', + dims=['basin'], + coords={'basin': range(len(MeltData))}, + ) + + t1_obs_sigma = MeltData['BMR uncert (Gt/yr)'].values + t1_obs_sigma = xr.DataArray( + data=t1_obs_sigma, + dims=['basin'], + coords={'basin': range(len(MeltData))}, + name='melt_unc_Gt_per_y', + ) + + return t1_model, t1_obs_mean, t1_obs_sigma + + +def calculate_term2(pd_ensemble, mask_m, bfrn_m, cvt_m, buttressing_target): + ######## + # TERM 2 + + # parameterisation melt aggregated per buttressing bin, in Gt/a + t2_model = ( + pd_ensemble['melt_rate'] + .where(mask_m, np.nan) + .groupby(bfrn_m['BFRN_bins']) + .sum(skipna=True) + * cvt_m + ) + t2_model = t2_model.where(t2_model != 0, np.nan) + + t2_obs_mean = buttressing_target['melt_mean'] + t2_obs_sigma = buttressing_target['melt_mean_err'] + + return t2_model, t2_obs_mean, t2_obs_sigma + + +def calculate_term3( + cold_ensemble, + warm_ensemble, + cold_target, + warm_target, + mask_m, + basins_m, +): + ######## + # TERM 3 + + # parameterisation melt, average in kg/m2/a per basin for cold ocean + t3_model_cold = ( + (cold_ensemble['melt_rate'].where(mask_m, np.nan)) + .groupby(basins_m) + .mean() + ) + t3_model_cold = t3_model_cold.where(t3_model_cold != 0, np.nan) + + # parameterisation melt, average in kg/m2/a per basin for warm ocean + t3_model_warm = ( + (warm_ensemble['melt_rate'].where(mask_m, np.nan)) + .groupby(basins_m) + .mean() + ) + t3_model_warm = t3_model_warm.where(t3_model_warm != 0, np.nan) + + t3_model = t3_model_warm - t3_model_cold + + t3_obs_mean = warm_target.melt_rate - cold_target.melt_rate + + t3_obs_sigma = np.sqrt( + warm_target.melt_rate_uncert**2 + cold_target.melt_rate_uncert**2 + ) + + return t3_model, t3_obs_mean, t3_obs_sigma + + +def calculate_term4(obs_ensemble, region_label_m, t4_obs, mask_m, cvt_m): + """ + Calculate Term 4 which is modelled melt in different years + for PIG and Dotson + """ + region_label_m = region_label_m.reindex_like(obs_ensemble) + + obs_ensemble_stacked = obs_ensemble.stack(grid=('x', 'y')) + mask_stacked = region_label_m.stack(grid=('x', 'y')) + t4_model = ( + obs_ensemble_stacked['melt_rate'].groupby(mask_stacked).sum(dim='grid') + * cvt_m + ) # Gt/a + t4_model = t4_model.rename({'group': 'region'}) + t4_model = t4_model.where(t4_model != 0, np.nan) + t4_model = t4_model.where(t4_obs.melt_rate.notnull()) + + t4_obs_mean = t4_obs.melt_rate + t4_obs_sigma = t4_obs.melt_rate_uncert + # make to order region index to pig, dotson + t4_model = t4_model.reindex_like(t4_obs_mean) + + return t4_model, t4_obs_mean, t4_obs_sigma + + +def optimise_deltaT(dT_ensemble, basins, reso, MeltDataImbie): + """ + Calculate optimal deltaT for basin-wide present-day melt. + """ + + number_of_basins = int(basins.max().values) + cvt = reso**2 / 1e12 # to convert to Gt/a + + optimal_deltaT_per_basin = [] + residual_per_basin = [] + sensitivity_per_basin = [] + + param_melt_rate = dT_ensemble.sel(deltaT=0).copy(deep=True) * np.nan + + for basin_i in range(number_of_basins + 1): + bmr = dT_ensemble.where(basins == basin_i, 0.0).sum(['x', 'y']) * cvt + + # only use deltaT between -2 and 2 + bmr = bmr.where( + np.logical_and(bmr['deltaT'] <= 2.0, bmr['deltaT'] >= -2.0), np.nan + ) + + optimal_deltaT_per_basin.append( + np.round( + (abs(bmr - MeltDataImbie.loc[basin_i, 'BMR (Gt/yr)'])) + .idxmin() + .item(), + 3, + ) + ) + residual_per_basin.append( + (abs(bmr - MeltDataImbie.loc[basin_i, 'BMR (Gt/yr)'])).min().item() + ) + + # if an optimal delatT exists, save melt and calc melt sensitivity + if not np.isnan(optimal_deltaT_per_basin[-1]): + param_melt_rate = param_melt_rate.where( + basins != basin_i, + dT_ensemble.sel( + deltaT=optimal_deltaT_per_basin[-1], method='nearest' + ), + ) + # Calc approx melt sensitivity. + # otherwise approximate with higher values, ideally +1deg C + # Note that this is only approximate + sensitivity_per_basin.append( + np.round( + ( + ( + dT_ensemble.sel( + deltaT=optimal_deltaT_per_basin[-1] + 1, + method='nearest', + ) + .where(basins == basin_i, np.nan) + .mean() + - param_melt_rate.where( + basins == basin_i, np.nan + ).mean() + ).item() + / 1 + ), + 2, + ) + ) + else: + sensitivity_per_basin.append(np.nan) + result_ds = xr.Dataset( + data_vars=dict( + melt_rate=(['y', 'x'], param_melt_rate.values), + optimal_deltaT_per_basin=( + ['basin'], + np.array(optimal_deltaT_per_basin), + ), + sensitivity_per_basin=(['basin'], np.array(sensitivity_per_basin)), + residual_per_basin=(['basin'], np.array(residual_per_basin)), + ), + coords=dict( + x=(['x'], dT_ensemble['x'].values), + y=(['y'], dT_ensemble['y'].values), + basin=(['basin'], np.arange(0, number_of_basins + 1)), + ), + ) + + return result_ds + + +def select_optimal_deltaT( + ds, basins, boxes, obs_data, param_type, outname, reso, ice_density, dT +): + """ + Only used for PICO + Input: + - ds: xarray dataset containing melt rates (m.i.e/a), with a dimension + deltaT that will be optimised over, on a regular grid + - basins: xarray dataset containing basin numbers, starting at 1 + - obs_data: data frame containing basin-aggregated melt rates in Gt/a + - param_type: pico, quadratic, ... + - outname: output file to save to + - reso: resolution of model output + - dT: allowed adjustment range + Output: + - saves a netcdf file to "outname" containing the melt rates for each + basin based on optimal deltaT, and arrays of optimal_deltaT, residuals + and melt sensitivities + """ + + # print('Identifying optimal delta T for each basin...') + + number_of_basins = int(basins.max().values) + cvt = reso**2 * ice_density / 1e12 # to convert to Gt/a + + optimal_deltaT_per_basin = [] + residual_per_basin = [] + sensitivity_per_basin = [] + + param_melt_rate = ds['melt_rate'].sel(deltaT=0).copy(deep=True) * np.nan + + for basin_i in range(1, number_of_basins + 1): + bmr = ( + ds['melt_rate'].where(basins == basin_i, 0.0).sum(['x', 'y']) * cvt + ) + if param_type == 'pico': + # Add physical constraints from Reese et al., 2018 + bmrBox1 = ( + ds['melt_rate'] + .where(np.logical_and(basins == basin_i, boxes == 1), 0.0) + .sum(['x', 'y']) + * cvt + ) + bmrBox2 = ( + ds['melt_rate'] + .where(np.logical_and(basins == basin_i, boxes == 2), 0.0) + .sum(['x', 'y']) + * cvt + ) + bmr = bmr.where( + np.logical_and(bmrBox1 > 0, bmrBox1 > bmrBox2), np.nan + ) + + # only use deltaT between +-dT + bmr = bmr.where( + np.logical_and(bmr['deltaT'] <= dT, bmr['deltaT'] >= -1 * dT), + np.nan, + ) + + optimal_deltaT_per_basin.append( + (abs(bmr - obs_data.loc[basin_i, 'BMR (Gt/yr)'])).idxmin() + ) + residual_per_basin.append( + (abs(bmr - obs_data.loc[basin_i, 'BMR (Gt/yr)'])).min() + ) + + # if an optimal delatT exists, save melt and calc melt sensitivity + if not np.isnan(optimal_deltaT_per_basin[-1]): + param_melt_rate = param_melt_rate.where( + basins != basin_i, + ds['melt_rate'].sel(deltaT=optimal_deltaT_per_basin[-1]), + ) + # Calc approx melt sensitivity. + # If this is the max deltaT, use half a degree colder, + # otherwise approximate with higher values, ideally +1deg C + # Note that this is only approximate + if optimal_deltaT_per_basin[-1] == ds.deltaT.max(): + sensitivity_per_basin.append( + ( + param_melt_rate.where(basins == basin_i, np.nan).mean() + - ds['melt_rate'] + .sel( + deltaT=optimal_deltaT_per_basin[-1] - 0.5, + method='nearest', + ) + .where(basins == basin_i, np.nan) + .mean() + ) + / 0.5 + ) + else: + sensitivity_per_basin.append( + ( + ds['melt_rate'] + .sel( + deltaT=optimal_deltaT_per_basin[-1] + 1, + method='nearest', + ) + .where(basins == basin_i, np.nan) + .mean() + - param_melt_rate.where( + basins == basin_i, np.nan + ).mean() + ) + / 1 + ) + else: + sensitivity_per_basin.append(np.nan) + + result_ds = xr.Dataset( + data_vars=dict( + melt_rate=(['y', 'x'], param_melt_rate.values), + optimal_deltaT_per_basin=( + ['basin'], + np.array(optimal_deltaT_per_basin), + ), + sensitivity_per_basin=(['basin'], np.array(sensitivity_per_basin)), + residual_per_basin=(['basin'], np.array(residual_per_basin)), + ), + coords=dict( + x=(['x'], ds['x'].values), + y=(['y'], ds['y'].values), + basin=(['basin'], np.arange(1, number_of_basins + 1)), + ), + ) + result_ds.to_netcdf(outname) + ds.close() + result_ds.close() + return result_ds.drop_vars('melt_rate') + + +def select_subensemble_using_optimal_deltaT( + ds, basins, opt_ensemble, outname, p1, p2 +): + """ + Input: + - ds: xarray dataset containing melt rates (m.i.e/a), + wit dimension deltaT that will be selected from, on a regular grid + - basins: xarray dataset containing basin numbers, starting at 1 + - opt_ensemble: array containing optimised deltaT's for present-day + - outname: output file to save to + Output: + - saves a netcdf file to "outname" containing the melt rates for each + basin based on optimal deltaT + """ + + # print('Select sub-ensemble...') + + number_of_basins = int(basins.max().values) + + # Create melt rate dataset based on optimal deltaT + melt_rate = ds['melt_rate'].sel(deltaT=0).copy(deep=True) * np.nan + + for basin in range(1, number_of_basins + 1): + optimal_deltaT = opt_ensemble['optimal_deltaT_per_basin'].loc[ + dict(p1=p1, p2=p2, basin=basin) + ] + + if np.isnan(optimal_deltaT.values): + melt_rate = melt_rate.where( + basins != basin, ds['melt_rate'].sel(deltaT=0) * np.nan + ) + else: + melt_rate = melt_rate.where( + basins != basin, ds['melt_rate'].sel(deltaT=optimal_deltaT) + ) + + result_ds = xr.Dataset( + data_vars=dict( + melt_rate=(['y', 'x'], melt_rate.values), + ), + coords=dict(x=(['x'], ds['x'].values), y=(['y'], ds['y'].values)), + ) + + result_ds.to_netcdf(outname) + + +def load_melt_rates_into_dataset( + ensemble_name, ensemble_table, ensemble_path, p1_name, p2_name, identifier +): + print('Loading ' + ensemble_name + ' into one dataset...') + members = [] + p1s = [] + p2s = [] + + for _i, ehash in enumerate(ensemble_table.index): + p1 = ensemble_table.loc[ehash, p1_name] + p2 = ensemble_table.loc[ehash, p2_name] + p1s.append(p1) + p2s.append(p2) + + print( + os.path.join( + ensemble_path, + ensemble_name + + '_' + + str(ehash) + + '/optimised' + + identifier + + '.nc', + ) + ) + + if os.path.isfile( + os.path.join( + ensemble_path, + ensemble_name + + '_' + + str(ehash) + + '/optimised' + + identifier + + '.nc', + ) + ): + ds = xr.load_dataset( + os.path.join( + ensemble_path, + ensemble_name + + '_' + + str(ehash) + + '/optimised' + + identifier + + '.nc', + ) + ) + elif os.path.isfile( + os.path.join( + ensemble_path, + ensemble_name + + '_' + + str(ehash) + + '_optimised' + + identifier + + '.nc', + ) + ): + ds = xr.load_dataset( + os.path.join( + ensemble_path, + ensemble_name + + '_' + + str(ehash) + + '_optimised' + + identifier + + '.nc', + ) + ) + else: + print('Error: Cannot find dataset') + + ds = ds.assign_coords(ehash=ehash) + members.append(ds) + + print('Combining datasets') + ensemble = xr.concat(members, dim='ehash', coords='minimal') + + ensemble = ( + ensemble.assign_coords({'p1': ('ehash', p1s), 'p2': ('ehash', p2s)}) + .set_index(ehash=['p1', 'p2']) + .unstack('ehash') + ) + return ensemble diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py index e03e85ecd6..04a6e5c400 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -2,16 +2,13 @@ import os import shutil -import numpy as np import xarray as xr from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call -from scipy.ndimage import distance_transform_edt -from compass.landice.tests.ismip7_forcing.create_mapfile import ( - build_mapping_file, -) -from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.ismip7.ice_sheet_params import get_params +from compass.landice.ismip7.mapping import build_mapping_file +from compass.landice.ismip7.remap import extrapolate_source from compass.step import Step @@ -123,8 +120,8 @@ def run(self): # so they don't pollute neighboring cells during interpolation extrap_file = f"extrap_{basename}" if not os.path.exists(extrap_file): - self._extrapolate_source(input_file, extrap_file, "mrro", - logger) + extrapolate_source(input_file, extrap_file, "mrro", + logger) logger.info(f" Remapping: {basename}") args = ["ncremap", @@ -225,56 +222,3 @@ def _combine_and_rename(self, remapped_files, output_file): ds = ds.drop_vars("Time") write_netcdf(ds, output_file) - - def _extrapolate_source(self, input_file, output_file, varname, logger): - """ - Extrapolate fill/missing values on the source polar stereographic - grid using nearest-neighbor via distance_transform_edt. This must - be done before remapping so that fill values don't contaminate the - interpolation stencil. - - Parameters - ---------- - input_file : str - Path to the input NetCDF file on the source grid - - output_file : str - Path to write the extrapolated file - - varname : str - Name of the variable to extrapolate (e.g., "mrro") - - logger : logging.Logger - Logger for status messages - """ - logger.info(f" Extrapolating fill values on source grid: " - f"{os.path.basename(input_file)}") - - ds = xr.open_dataset(input_file, engine="netcdf4") - data = ds[varname] - - # Process each time step - # Source files have dims like (time, y, x) - values = data.values.copy() - non_spatial_shape = values.shape[:-2] # (time,) - - for idx in np.ndindex(non_spatial_shape): - slab = values[idx] # shape (ny, nx) - valid_mask = np.isfinite(slab) - if valid_mask.all() or not valid_mask.any(): - continue - nearest_inds = distance_transform_edt( - ~valid_mask, return_distances=False, return_indices=True) - invalid = ~valid_mask - values[idx][invalid] = slab[ - nearest_inds[0, invalid], - nearest_inds[1, invalid]] - - ds[varname] = (data.dims, values) - ds[varname].attrs = data.attrs - - # Remove _FillValue encoding so output has no masked values - if "_FillValue" in ds[varname].encoding: - del ds[varname].encoding["_FillValue"] - - write_netcdf(ds, output_file) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py index 222fb6e92e..44174b44a6 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -6,10 +6,8 @@ from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call -from compass.landice.tests.ismip7_forcing.create_mapfile import ( - build_mapping_file, -) -from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.ismip7.ice_sheet_params import get_params +from compass.landice.ismip7.mapping import build_mapping_file from compass.step import Step diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py index 2f3e4c2cb4..e7ca7a1f3e 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -6,10 +6,8 @@ from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call -from compass.landice.tests.ismip7_forcing.create_mapfile import ( - build_mapping_file, -) -from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.ismip7.ice_sheet_params import get_params +from compass.landice.ismip7.mapping import build_mapping_file from compass.step import Step diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index caa8f6be9f..bbf2c16b48 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -6,10 +6,8 @@ from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call -from compass.landice.tests.ismip7_forcing.create_mapfile import ( - build_mapping_file, -) -from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.ismip7.ice_sheet_params import get_params +from compass.landice.ismip7.mapping import build_mapping_file from compass.step import Step diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py index 5dec7b78bc..f12e212753 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -6,10 +6,8 @@ from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call -from compass.landice.tests.ismip7_forcing.create_mapfile import ( - build_mapping_file, -) -from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.ismip7.ice_sheet_params import get_params +from compass.landice.ismip7.mapping import build_mapping_file from compass.step import Step diff --git a/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py b/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py index fbf3f75b83..bc55cebc99 100644 --- a/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py +++ b/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py @@ -7,10 +7,8 @@ from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call -from compass.landice.tests.ismip7_forcing.create_mapfile import ( - build_mapping_file, -) -from compass.landice.tests.ismip7_forcing.fracture.remap_utils import ( +from compass.landice.ismip7.mapping import build_mapping_file +from compass.landice.ismip7.remap import ( add_xtime_and_write, extrapolate_source, open_rename_and_trim, @@ -122,7 +120,8 @@ def run(self): # Extrapolate fill values on the source grid before remapping so # they don't pollute neighboring cells during interpolation extrap_file = f"extrap_{basename}" - extrapolate_source(gridded_file, extrap_file, "excess_melt", logger) + extrapolate_source(gridded_file, extrap_file, "excess_melt", + logger, decode_times=False) # Remap the excess melt onto the MALI mesh remapped_file = f"remapped_{basename}" diff --git a/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py b/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py index 10df254c69..d9b3838a14 100644 --- a/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py +++ b/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py @@ -4,10 +4,8 @@ from mpas_tools.logging import check_call -from compass.landice.tests.ismip7_forcing.create_mapfile import ( - build_mapping_file, -) -from compass.landice.tests.ismip7_forcing.fracture.remap_utils import ( +from compass.landice.ismip7.mapping import build_mapping_file +from compass.landice.ismip7.remap import ( add_xtime_and_write, extrapolate_source, open_rename_and_trim, @@ -124,7 +122,8 @@ def run(self): # they don't pollute neighboring cells during interpolation extrap_file = f"extrap_{basename}" extrapolate_source(input_file, extrap_file, - list(self._variables.keys()), logger) + list(self._variables.keys()), logger, + decode_times=False) # Remap both lake property variables onto the MALI mesh remapped_file = f"remapped_{basename}" diff --git a/compass/landice/tests/ismip7_forcing/fracture/process_shelf_collapse.py b/compass/landice/tests/ismip7_forcing/fracture/process_shelf_collapse.py index e9a8e1cfb5..6bf71fd40c 100644 --- a/compass/landice/tests/ismip7_forcing/fracture/process_shelf_collapse.py +++ b/compass/landice/tests/ismip7_forcing/fracture/process_shelf_collapse.py @@ -4,10 +4,8 @@ from mpas_tools.logging import check_call -from compass.landice.tests.ismip7_forcing.create_mapfile import ( - build_mapping_file, -) -from compass.landice.tests.ismip7_forcing.fracture.remap_utils import ( +from compass.landice.ismip7.mapping import build_mapping_file +from compass.landice.ismip7.remap import ( add_xtime_and_write, open_rename_and_trim, ) diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py index a1469c8f2d..f2dd357de4 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -2,16 +2,13 @@ import os import shutil -import numpy as np import xarray as xr from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call -from scipy.ndimage import distance_transform_edt -from compass.landice.tests.ismip7_forcing.create_mapfile import ( - build_mapping_file, -) -from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.ismip7.ice_sheet_params import get_params +from compass.landice.ismip7.mapping import build_mapping_file +from compass.landice.ismip7.remap import extrapolate_source from compass.step import Step @@ -150,8 +147,8 @@ def _run_scenario(self): # so they don't pollute neighboring cells during interpolation extrap_file = f"extrap_{basename}" if not os.path.exists(extrap_file): - self._extrapolate_source(input_file, extrap_file, "tf", - logger) + extrapolate_source(input_file, extrap_file, "tf", + logger) logger.info(f" Remapping: {basename}") args = ["ncremap", @@ -246,8 +243,8 @@ def _run_climatology(self): if not os.path.exists(remapped_file): extrap_file = f"extrap_{basename}" if not os.path.exists(extrap_file): - self._extrapolate_source(input_file, extrap_file, "tf", - logger) + extrapolate_source(input_file, extrap_file, "tf", + logger) logger.info(f" Remapping: {basename}") args = ["ncremap", @@ -563,59 +560,3 @@ def _rename_climatology_3d(self, remapped_file, output_file): ds = ds.drop_vars("nISMIP6OceanLayers") write_netcdf(ds, output_file) - - def _extrapolate_source(self, input_file, output_file, varname, logger): - """ - Extrapolate fill/missing values on the source polar stereographic - grid using nearest-neighbor interpolation from valid cells. This - must be done before remapping so that fill values don't contaminate - the interpolation stencil. - - Parameters - ---------- - input_file : str - Path to the input NetCDF file on the source grid - - output_file : str - Path to write the extrapolated file - - varname : str - Name of the variable to extrapolate (e.g., "tf") - - logger : logging.Logger - Logger for status messages - """ - logger.info(f" Extrapolating fill values on source grid: " - f"{os.path.basename(input_file)}") - - ds = xr.open_dataset(input_file, engine="netcdf4") - data = ds[varname] - - # Process each time step (and z level if 3D) - # Source files have dims like (time, z, y, x) or (time, y, x) - values = data.values.copy() - non_spatial_shape = values.shape[:-2] # (time,) or (time, z) - - # Use distance_transform_edt with return_indices to find the - # nearest valid cell index for each invalid cell. This is O(n) - # on the grid and much faster than KD-tree approaches. - for idx in np.ndindex(non_spatial_shape): - slab = values[idx] # shape (ny, nx) - valid_mask = np.isfinite(slab) - if valid_mask.all() or not valid_mask.any(): - continue - nearest_inds = distance_transform_edt( - ~valid_mask, return_distances=False, return_indices=True) - invalid = ~valid_mask - values[idx][invalid] = slab[ - nearest_inds[0, invalid], - nearest_inds[1, invalid]] - - ds[varname] = (data.dims, values) - ds[varname].attrs = data.attrs - - # Remove _FillValue encoding so output has no masked values - if "_FillValue" in ds[varname].encoding: - del ds[varname].encoding["_FillValue"] - - write_netcdf(ds, output_file) diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index 781c261909..6656b4627b 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -18,6 +18,20 @@ Utilities calculate_decomp_core_pair +ISMIP7 framework +^^^^^^^^^^^^^^^^ + +.. currentmodule:: compass.landice.ismip7 + +.. autosummary:: + :toctree: generated/ + + ice_sheet_params.get_params + mapping.build_mapping_file + remap.extrapolate_source + remap.open_rename_and_trim + remap.add_xtime_and_write + Test Groups ^^^^^^^^^^^ @@ -373,6 +387,93 @@ ismip6_run ismip6_ais_proj2300.set_up_experiment.SetUpExperiment.setup ismip6_ais_proj2300.set_up_experiment.SetUpExperiment.run +ismip7_calibration +~~~~~~~~~~~~~~~~~~ + +.. currentmodule:: compass.landice.tests.ismip7_calibration + +.. autosummary:: + :toctree: generated/ + + Ismip7Calibration + configure.check_options + configure.melt_forms + configure.objective_options + configure.parameter_name + configure.parameter_values + configure.weighting + datasets.climatology_files + datasets.load_targets + datasets.mask_files + datasets.missing_files + datasets.ocean_states + objective.build_toolbox_terms + objective.run_optimisation + objective.scale_to_ensemble + quadratic.angle_from_sin_slope + quadratic.draft_slope + quadratic.local_quadratic_melt + quadratic.mean_slope + quadratic.nonlocal_quadratic_melt + quadratic.u_factor + terms.average_by_group + terms.calculate_term1 + terms.calculate_term2 + terms.calculate_term3 + terms.calculate_term4 + terms.integrate_by_group + terms.stack_cells + terms.uniform_area + toolbox.check_integrity + toolbox.file_sha256 + toolbox.toolbox_path + + replication.Replication + replication.Replication.configure + replication.Replication.validate + replication.replicate.Replicate + replication.replicate.Replicate.setup + replication.replicate.Replicate.run + replication.replicate.check_published + + ais.Ais + ais.Ais.configure + ais.Ais.validate + ais.aggregate.Aggregate + ais.aggregate.Aggregate.setup + ais.aggregate.Aggregate.run + ais.aggregate.unit_aggregates + ais.aggregate.basin_coordinate + ais.calibrate.Calibrate + ais.calibrate.Calibrate.run + ais.fit_delta_t.FitDeltaT + ais.fit_delta_t.FitDeltaT.setup + ais.fit_delta_t.FitDeltaT.run + ais.make_graph.MakeGraph + ais.make_graph.MakeGraph.setup + ais.make_graph.MakeGraph.run + ais.melt_model.basin_mean_tf + ais.melt_model.initial_draft + ais.melt_model.integrate_by_basin + ais.melt_model.interpolate_to_draft + ais.melt_model.melt_from_tf + ais.melt_model.read_run + ais.remap_forcing.RemapForcing + ais.remap_forcing.RemapForcing.setup + ais.remap_forcing.RemapForcing.run + ais.remap_masks.RemapMasks + ais.remap_masks.RemapMasks.setup + ais.remap_masks.RemapMasks.run + ais.report.Report + ais.report.Report.run + ais.run_state.RunState + ais.run_state.RunState.setup + ais.run_state.RunState.runtime_setup + ais.run_state.RunState.run + ais.verify_melt.VerifyMelt + ais.verify_melt.VerifyMelt.setup + ais.verify_melt.VerifyMelt.run + ismip7_forcing ~~~~~~~~~~~~~~ @@ -383,8 +484,6 @@ ismip7_forcing Ismip7Forcing configure.configure - ice_sheet_params.get_params - create_mapfile.build_mapping_file atmosphere.Atmosphere atmosphere.Atmosphere.configure @@ -412,9 +511,6 @@ ismip7_forcing fracture.Fracture fracture.Fracture.configure - fracture.remap_utils.extrapolate_source - fracture.remap_utils.open_rename_and_trim - fracture.remap_utils.add_xtime_and_write fracture.process_excess_melt.ProcessExcessMelt fracture.process_excess_melt.ProcessExcessMelt.setup fracture.process_excess_melt.ProcessExcessMelt.run diff --git a/docs/developers_guide/landice/framework.rst b/docs/developers_guide/landice/framework.rst index 712eed33ba..009b2dd42a 100644 --- a/docs/developers_guide/landice/framework.rst +++ b/docs/developers_guide/landice/framework.rst @@ -18,6 +18,35 @@ The landice framework module ``compass/landice/extrapolate.py`` provides a function for extrapolating variables into undefined regions. It is copied from a similar script in MPAS-Tools. +ismip7 +~~~~~~ + +The landice framework package :py:mod:`compass.landice.ismip7` holds code +shared by the ISMIP7 test groups -- ``ismip7_forcing``, ``ismip7_run`` and +``ismip7_calibration``. All of them remap data from an ISMIP7 polar +stereographic grid onto a MALI mesh, so the helpers to do that live here +rather than in any one test group. + +:py:func:`compass.landice.ismip7.ice_sheet_params.get_params()` returns the +parameters that differ between the Antarctic and Greenland ice sheets: the +projection, the filename prefix, dataset versions and resolutions, and whether +the ocean forcing is 3-D. + +:py:func:`compass.landice.ismip7.mapping.build_mapping_file()` builds an ESMF +mapping file from an ISMIP7 polar stereographic grid to a MALI mesh, using +``mpas_tools.scrip.from_mpas`` and ``ESMF_RegridWeightGen``. + +:py:func:`compass.landice.ismip7.remap.extrapolate_source()` fills missing +values on the source grid by nearest neighbour before remapping, so that fill +values do not contaminate the interpolation stencil. Pass +``decode_times=False`` for source files whose time coordinate is not +CF-compliant, such as the fracture forcing. + +:py:func:`compass.landice.ismip7.remap.open_rename_and_trim()` and +:py:func:`compass.landice.ismip7.remap.add_xtime_and_write()` convert a +remapped file to MALI's variable and dimension names, restrict it to a range +of years, and add the ``xtime`` variable MALI needs. + iceshelf_melt ~~~~~~~~~~~~~ The landice framework module ``compass/landice/iceshelf_melt.py`` provides diff --git a/docs/developers_guide/landice/test_groups/index.rst b/docs/developers_guide/landice/test_groups/index.rst index b0871ac2b5..374719877c 100644 --- a/docs/developers_guide/landice/test_groups/index.rst +++ b/docs/developers_guide/landice/test_groups/index.rst @@ -20,6 +20,7 @@ Test groups hydro_radial ismip6_forcing ismip6_run + ismip7_calibration ismip7_forcing ismip7_run isunnguata_sermia diff --git a/docs/developers_guide/landice/test_groups/ismip7_calibration.rst b/docs/developers_guide/landice/test_groups/ismip7_calibration.rst new file mode 100644 index 0000000000..f0ae151a3a --- /dev/null +++ b/docs/developers_guide/landice/test_groups/ismip7_calibration.rst @@ -0,0 +1,284 @@ +.. _dev_landice_ismip7_calibration: + +ismip7_calibration +================== + +The ``ismip7_calibration`` test group +(:py:class:`compass.landice.tests.ismip7_calibration.Ismip7Calibration`) +calibrates 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 fit the free parameter of its melt +module against four objective-function terms and to report the 5th, 50th and +95th percentiles of the resulting parameter distribution: + +===== ========================================================== +term what it constrains +===== ========================================================== +J1 present-day melt integrated over each IMBIE2 drainage basin +J2 present-day melt integrated over bins of equal buttressing importance +J3 the warm-minus-cold basin-mean melt difference of ocean models +J4 observed melt of Pine Island and Dotson, per observation year +===== ========================================================== + +The user's guide describes the config options in +:ref:`landice_ismip7_calibration`. + +The test group has two test cases: ``replication`` and ``ais``. + +.. _dev_landice_ismip7_calibration_framework: + +framework +--------- + +Code shared with the other ISMIP7 test groups lives in the landice framework +package :py:mod:`compass.landice.ismip7`, described in +:ref:`dev_landice_framework`. The modules below are specific to the +calibration. + +datasets +~~~~~~~~ + +:py:mod:`compass.landice.tests.ismip7_calibration.datasets` is the single +registry of *which* ocean states feed *which* objective-function term, and of +where the calibration targets live. Both test cases drive it, so the 8 km +replication and the MALI-mesh calibration cannot drift apart. + +:py:func:`compass.landice.tests.ismip7_calibration.datasets.ocean_states` +returns the ``minimal`` (7), ``recommended`` (11) or ``all`` (28) subset of +protocol Table 2. Regional ocean models record the basins they cover, so +that melt outside their domains is discarded. + +terms +~~~~~ + +:py:mod:`compass.landice.tests.ismip7_calibration.terms` implements J1-J4 in +a form that works on any mesh. The upstream toolbox assumes a uniform +structured grid: it converts cell sums to Gt/yr with a single scalar +``reso**2`` and takes unweighted means over basins. That is correct on the +ISMIP 8 km grid and wrong on MALI's 4-20 km mesh, where cell area varies by a +factor of about 25. These functions take an explicit cell-area array +instead, so passing a uniform area reproduces the structured-grid answer +exactly while a real mesh gets an area-weighted one. + +quadratic +~~~~~~~~~ + +:py:mod:`compass.landice.tests.ismip7_calibration.quadratic` implements both +melt formulas in Python: the Burgard et al. (2022) local quadratic of +protocol Eq. (1), and the ISMIP6 non-local quadratic. It serves as the +specification MALI's Fortran is checked against and as what the replication +drives. + +Note that ``local_quadratic_melt`` takes the slope *angle* while MALI's +``config_ismip7_melt_sin_slope`` is a *sine*; +:py:func:`compass.landice.tests.ismip7_calibration.quadratic.angle_from_sin_slope` +converts between them. + +toolbox +~~~~~~~ + +:py:mod:`compass.landice.tests.ismip7_calibration.toolbox` holds a verbatim +copy of the upstream ISMIP7 parameter-selection toolbox, which implements the +objective function itself. ``PROVENANCE.md`` alongside it records the +upstream commit and a SHA256 that is verified on import and asserted by a +unit test, so an accidental edit or a partial update fails loudly rather than +quietly changing published numbers. + +objective +~~~~~~~~~ + +:py:mod:`compass.landice.tests.ismip7_calibration.objective` assembles the +toolbox's arguments and reduces its output to percentiles. Each term is +built once at a unit parameter value and scaled onto the parameter grid, +which is exact because melt is proportional to the parameter. + +.. _dev_landice_ismip7_calibration_replication: + +replication +----------- + +``landice/ismip7_calibration/replication`` reproduces the published 8 km +calibration through compass's own code path, using the Python reference melt +implementation on the ISMIP grid. It runs in about a minute, needs no MALI +run, and must reproduce the published percentiles + +.. code-block:: none + + K = 4.75e-5 / 8.5e-5 / 1.375e-4 + +exactly. ``validate()`` raises if it does not. This is the check that the +vendored toolbox is being driven correctly; the MALI calibration is built on +the same code path, so if this fails the MALI numbers cannot be trusted +either. + +This test case deliberately uses the **published inputs**, including the +spatially varying 8 km salinity fields. It replicates a published +calculation, so it must not be changed to match the constant-salinity choice +the MALI calibration makes. + +.. _dev_landice_ismip7_calibration_ais: + +ais +--- + +``landice/ismip7_calibration/ais`` is the calibration itself, on an Antarctic +MALI mesh. + +``make_graph`` + Builds the graph partition file the melt diagnostics run on, once, from + the mesh. Doing it here means the test group works on any MALI mesh + rather than only on meshes that already have a partition file + distributed alongside them. + +``remap_masks`` + Remaps the ISMIP7 IMBIE2 basins, buttressing (BFRN) bins, floating mask + and PIG/Dotson regions onto the MALI mesh, nearest-neighbour throughout + since every field is categorical. + + **Two basin-numbering conventions are in play and they differ by one.** + ISMIP7's ``basinNumber`` is 0-based, 0-15, with basin 9 the Eastern + Amundsen and basin 14 Ronne-Filchner, which is how the protocol refers to + them. MALI's ``ismip6shelfMelt_basin`` is 1-based, 1-16, so MALI basin + 10 is ISMIP7 basin 9. Both are written, under distinct names, so that + neither can be silently reinterpreted as the other -- a mistake that + mis-assigns every basin while still producing plausible-looking numbers. + The step cross-tabulates the result against the mesh's existing + ``regionCellMasks`` in both directions and fails below 80% agreement. + +``remap_forcing`` + Remaps the calibration thermal forcing for each ocean state. The + observational states cover only about 6% of the ISMIP grid (protocol + Sect. A10), so they are filled from the present-day climatology *before* + remapping, exactly as the regional model datasets are distributed + (Sect. A9); interpolation then never blends a real value with a missing + one. The remapped field is required to be finite everywhere. + + **Only thermal forcing is remapped.** See the salinity note below. + +``_`` + One single-timestep MALI melt diagnostic per ocean state and melt form, + with the velocity solver off, so Albany is not required. + + The run takes **one short timestep rather than none**: MALI computes melt + inside the timestep, not in the initial diagnostic solve, so a zero-length + run would produce no melt at all. + + There is one run per ocean state, **not** one per parameter value. Melt + is exactly proportional to the melt parameter, so the whole parameter + ensemble follows by scaling a single run. That is what makes this 28 runs + per melt form rather than about 1300. Do not "improve" this away. + +``__x`` + Two extra runs of the reference ocean state at other multiples of the + melt parameter, so that the linearity the ensemble design relies on is + measured in MALI rather than argued from the code. Only the ``ismip7`` + form is scaled this way: its parameter is a namelist option, while the + ISMIP6 ``gamma0`` is read from an input file. + +``verify_melt`` + Checks MALI's melt against the Python reference evaluated on MALI's *own* + ``TFdraft``, which isolates the melt expression from the vertical + interpolation; checks that interpolation against an independent + implementation written from the protocol, over all four of MALI's code + paths; and measures the linearity in the melt parameter from the scaled + runs above. + + The comparison is restricted to the cells MALI itself computes melt for + -- floating *and* connected to the open ocean. MALI leaves ``TFdraft`` + and the melt at zero elsewhere, so including those cells would compare + against values MALI never calculated. + + The draft is reconstructed from the mesh file, never taken from the run + output. Melt thins the ice over the single timestep, so the output + ``lowerSurface`` and ``thickness`` are *post*-step while ``TFdraft`` was + computed *pre*-step; pairing them is inconsistent by up to a metre of + draft. + + MALI applies the depth dependence of the freezing temperature to the + thermal forcing in the two branches that extrapolate downward -- below + the deepest layer centre, and where the layer below the draft is beneath + the bed. The Python reference must do the same, or the two disagree by + several kelvin over thousands of cells. + +``aggregate`` + Aggregates melt to basins, buttressing bins and shelf regions. All + aggregation is **area-weighted**: integrals use ``melt * areaCell`` and + means are weighted by ``areaCell``. Also reports MALI's per-basin + ice-shelf area against the observed ISMIP7 extent, since J1, J2 and J4 + are integrals and any shelf-area mismatch enters the calibrated parameter + directly. + + The contributing cells are the ones MALI computed melt for, as in + ``verify_melt``. Including the rest would not change an integral, since + their melt is zero, but it would dilute the area-weighted basin means + that J3 is built from. + +``calibrate`` + Runs the 100,000-sample parameter selection per melt form. The objective + **normalises each term by its own median** over the parameter ensemble, so + a minimised objective is *not* comparable between melt forms -- only + between parameter values within one form. + +``fit_delta_t`` + Fits the per-basin thermal-forcing correction ``dT_b`` **after** parameter + selection, following protocol Sect. 4.2.1 option 2, as the published + quadratic worked example does. Fitting it first would leave the parameter + bounds unconstrained where present-day melt is compared with observations, + which is the drawback Sect. 4.2.1 names. The toolbox's own + ``optimise_deltaT`` cannot be used because it sums over structured-grid + ``x`` and ``y`` with a single scalar cell area, so the search is + reimplemented area-weighted. ``dT`` enters the melt only through the + thermal forcing, so melt is recomputed in Python from ``TFdraft`` and no + further MALI runs are needed. + +``report`` + Parameter-distribution plots and a summary table. + +melt forms and salinity +----------------------- + +Two melt forms are calibrated, selected with the ``melt_forms`` config +option: + +``ismip7`` + The Burgard et al. (2022) **local** quadratic of protocol Eq. (1), + calibrating ``K``. + +``ismip6`` + The **non-local** quadratic MALI already had, calibrating ``gamma0``. + With a constant salinity this is algebraically identical to the Burgard + *semi-local* form of protocol Eq. (2) -- only the decomposition of the + constant differs -- so this is how the semi-local form is calibrated, + rather than as a third melt module. + +**The calibration uses a constant salinity**, set by the ``salinity`` config +option and passed to MALI as +``config_ismip7_melt_salinity_source = 'constant'``. MALI can read a 3-D +salinity field, but the ISMIP7 ocean forcing already processed for the MALI +projections carries thermal forcing only, so a projection has no salinity +field to read. Calibrating against a spatially varying salinity would tune +``K`` for physics the projections cannot run. + +Melt is linear in salinity and basin-mean salinity spans about 34.1 to 34.7 +against the constant 34.5, so this shifts the calibrated parameter by roughly +1%. Note that the published ``K`` was obtained with a locally varying +salinity, so comparison with it is approximate at about that level. + +using a MALI build with the ISMIP7 melt method +---------------------------------------------- + +The ``ismip7`` melt form needs a MALI build that has +``config_basal_mass_bal_float = 'ismip7'``. If the compass ``MALI-Dev`` +submodule does not yet have it, point compass at another build with + +.. code-block:: cfg + + [paths] + mpas_model = /path/to/E3SM/components/mpas-albany-landice + +Compass only *warns* when a namelist option is missing from the model's +defaults, so a build without the ISMIP7 melt method would silently drop +``config_ismip7_melt_K`` and produce a plausible but wrong calibration. The +``run_state`` step therefore checks the default namelist at setup and fails +with a clear message instead. diff --git a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst index b535bf081a..648eb2c49a 100644 --- a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst @@ -19,15 +19,15 @@ framework The shared config options for the ``ismip7_forcing`` test group are described in :ref:`landice_ismip7_forcing` in the User's Guide. -ice_sheet_params -~~~~~~~~~~~~~~~~ - -The module :py:mod:`compass.landice.tests.ismip7_forcing.ice_sheet_params` -defines a dictionary of ice-sheet-specific parameters (projection, file -naming prefix, grid resolution, data version, ocean dimensionality) and -provides the function -:py:func:`compass.landice.tests.ismip7_forcing.ice_sheet_params.get_params` -to retrieve them based on the ``ice_sheet`` config option. +Code shared with the other ISMIP7 test groups lives in the landice framework +package :py:mod:`compass.landice.ismip7`, described in +:ref:`dev_landice_framework`. This test group uses +:py:func:`compass.landice.ismip7.ice_sheet_params.get_params` for the +ice-sheet-specific parameters (projection, file naming prefix, grid +resolution, data version, ocean dimensionality), +:py:func:`compass.landice.ismip7.mapping.build_mapping_file` to create the +SCRIP and ESMF mapping files, and the remapping helpers in +:py:mod:`compass.landice.ismip7.remap`. configure ~~~~~~~~~ @@ -42,18 +42,6 @@ Repository-local example user configs are available at These are intended for development/testing and include environment-specific paths. -create_mapfile -~~~~~~~~~~~~~~ - -The module :py:mod:`compass.landice.tests.ismip7_forcing.create_mapfile` -defines a unified framework for creating SCRIP and mapping files. The function -:py:func:`compass.landice.tests.ismip7_forcing.create_mapfile.build_mapping_file` -creates a SCRIP file from the input polar stereographic grid using the -``create_scrip_file_from_planar_rectangular_grid`` command from MPAS-Tools, -then generates a mapping file via ``ESMF_RegridWeightGen``. The projection -is automatically determined from the ``ice_sheet`` config option using -``ice_sheet_params``. - Test cases ---------- @@ -163,7 +151,7 @@ The annual source fields use an integer ``year``/``time`` coordinate with ``decode_times=False`` and constructs ``xtime`` at January 1st of each year. Shared remapping helpers used by the fracture steps live in -:py:mod:`compass.landice.tests.ismip7_forcing.fracture.remap_utils`: +:py:mod:`compass.landice.ismip7.remap`: ``extrapolate_source`` (nearest-neighbor fill of NaNs on the source grid), ``open_rename_and_trim`` (open a remapped file, rename dimensions/variables to MALI conventions, and restrict to the requested year range), and diff --git a/docs/users_guide/landice/test_groups/index.rst b/docs/users_guide/landice/test_groups/index.rst index 6a49090c00..ea52774f02 100644 --- a/docs/users_guide/landice/test_groups/index.rst +++ b/docs/users_guide/landice/test_groups/index.rst @@ -25,6 +25,7 @@ physics but that are not run routinely. hydro_radial ismip6_forcing ismip6_run + ismip7_calibration ismip7_forcing ismip7_run isunnguata_sermia diff --git a/docs/users_guide/landice/test_groups/ismip7_calibration.rst b/docs/users_guide/landice/test_groups/ismip7_calibration.rst new file mode 100644 index 0000000000..8d79f42c7f --- /dev/null +++ b/docs/users_guide/landice/test_groups/ismip7_calibration.rst @@ -0,0 +1,176 @@ +.. _landice_ismip7_calibration: + +ismip7_calibration +================== + +The ``landice/ismip7_calibration`` test group calibrates MALI's sub-shelf +melt parameterization against the Ice Sheet Model Intercomparison for CMIP7 +(ISMIP7) Antarctic ice-ocean protocol (Reese et al., Sect. 4.2), and produces +the 5th, 50th and 95th percentiles of the melt parameter that the ISMIP7 MALI +projections need. + +The protocol fits the free parameter of the melt module against four terms: + +===== ========================================================== +term what it constrains +===== ========================================================== +J1 present-day melt integrated over each IMBIE2 drainage basin +J2 present-day melt integrated over bins of equal buttressing importance +J3 the warm-minus-cold basin-mean melt difference of ocean models +J4 observed melt of Pine Island and Dotson, per observation year +===== ========================================================== + +It then draws 100,000 random samples of the term weights and of the targets +within their uncertainties, minimising the objective for each draw. The +distribution of the minimisers gives the percentiles. + +The test group includes two test cases. + +* ``replication`` reproduces the published 8 km calibration through compass's + own code path. It runs in about a minute on one core, needs no MALI run, + and needs only the ISMIP7 datasets. Use it to check that everything is + wired up before running the full calibration. + +* ``ais`` is the calibration itself on an Antarctic MALI mesh: it remaps the + ISMIP7 masks and forcing, runs one single-timestep MALI melt diagnostic per + ocean state, aggregates the melt, selects the parameter and fits the + per-basin thermal-forcing correction ``dT_b``. + +Two melt forms can be calibrated, and by default both are: + +``ismip7`` + the Burgard et al. (2022) **local** quadratic recommended for ISMIP7, + calibrating ``K`` + +``ismip6`` + the **non-local** quadratic MALI already had, calibrating ``gamma0``. + With a constant salinity this is algebraically the same as the Burgard + *semi-local* form, so this is how that form is calibrated. + +.. _landice_ismip7_calibration_salinity: + +A note on salinity +------------------ + +The calibration runs MALI with a **constant** ocean salinity. MALI can read +a 3-D salinity field, but the ISMIP7 ocean forcing already processed for the +MALI projections carries thermal forcing only, so a projection has no +salinity field to read. Calibrating against a spatially varying salinity +would tune the melt parameter for physics the projections cannot run. + +Melt is linear in salinity, and basin-mean salinity spans about 34.1 to 34.7 +against the default constant of 34.5, so this shifts the calibrated parameter +by roughly 1%. The published ``K`` was obtained with a locally varying +salinity, so comparison with it is approximate at about that level. + +.. _landice_ismip7_calibration_usage: + +Usage +----- + +The test group needs the ISMIP7 AIS datasets, a MALI mesh with its graph +partition file, and a MALI build that supports the ISMIP7 melt method. +Supply them in a user config file: + +.. code-block:: cfg + + [paths] + # only needed if the compass MALI-Dev submodule does not yet have the + # ISMIP7 melt method + mpas_model = /path/to/E3SM/components/mpas-albany-landice + + [ismip7_calibration] + base_path_ismip7 = /path/to/ISMIP7/data/AIS + base_path_mali = /path/to/inputdata/glc/mpasli/mpas.ais4to20km + mali_mesh_file = ais_4to20km.20250625.nc + mali_mesh_name = ais_4to20km + # optional; set to None if no ISMIP6-era region mask exists on this + # mesh, which skips the basin cross-check + region_mask_file = ais_4to20km_region_mask.20230105.nc + +Then set up and run, for example: + +.. code-block:: bash + + compass setup -t landice/ismip7_calibration/replication \ + -f my.cfg -w $WORKDIR + cd $WORKDIR/landice/ismip7_calibration/replication + sbatch job_script.sh + +Setting up needs nothing but the config file above; every part of the +calibration is a step of the test case, so ``compass run`` from the test case +directory does the whole thing in order, and the job script ``compass setup`` +writes runs exactly that. + +The ``ais`` test case creates one step per ocean state and melt form -- 56 +with the default settings -- plus two more for the linearity check. Each is a +short MALI run, but compass runs steps one after another, so the ensemble +takes a few hours of wall clock. The runs do not depend on each other, so if +that matters, run ``compass run`` from each step directory in separate jobs +instead, and then run the ``verify_melt``, ``aggregate``, ``calibrate``, +``fit_delta_t`` and ``report`` steps in that order. + +Config options +-------------- + +The default config options are: + +.. code-block:: cfg + + [ismip7_calibration] + + # Which ocean states to use. Options: + # minimal - 7 states, those marked X in protocol Table 2 + # recommended - 11 states, adding the regional pairs marked + + # all - 28 states, everything ISMIP7 distributes + ocean_state_subset = all + + # Which melt forms to calibrate, comma separated + melt_forms = ismip7, ismip6 + + # Number of MPI tasks for ESMF_RegridWeightGen + esmf_ntasks = 128 + + # Number of MPI tasks for each MALI run + ntasks = 128 + + # Length of the single timestep the melt diagnostic takes + timestep = 0000-00-01_00:00:00 + + [ismip7_calibration_melt] + + # Practical salinity at the ice draft, PSU + salinity = 34.5 + + # The sin(theta) factor of the ISMIP7 quadratic + sin_slope = 0.0051117 + + # Magnitude of the Coriolis parameter, s^-1 + coriolis = 1.4e-4 + + [ismip7_calibration_objective] + + # Number of random draws of the term weights and the targets + sample_size = 100000 + + # Seed for the random draws, so the percentiles are reproducible + seed = 0 + + # Which summands carry non-zero weight: 'published' or 'all' + t3_models = all + t4_regions = all + t4_years = all + +See ``compass/landice/tests/ismip7_calibration/ismip7_calibration.cfg`` for +the full set with comments. + +The ``t4_regions`` option is worth understanding. The published weighting +uses Pine Island alone, which is 2 of the 18 available observations. So +weighted, J4 constrains a single amplitude that either melt form can match by +rescaling its parameter, so it cannot discriminate between them. Including +Dotson makes J4 a relative constraint between two shelves, which no rescaling +can satisfy if the ratio is wrong. The default is therefore ``all``; set +``published`` to reproduce the published weighting. + +Similarly, the reduced ocean-state subsets leave J4 structurally +under-powered, which is why ``ocean_state_subset`` defaults to ``all``. diff --git a/pyproject.toml b/pyproject.toml index a1d690588d..db4a509da9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -86,6 +86,12 @@ compass = [ [tool.setuptools.dynamic] version = { attr = "compass.version.__version__" } +[tool.pytest.ini_options] +testpaths = ["tests"] +# unit tests only: they must run in seconds with no input datasets, no MPAS +# build and no network access, because CI has none of those +addopts = "--strict-markers" + [tool.mypy] python_version = "3.14" check_untyped_defs = true @@ -95,6 +101,12 @@ warn_redundant_casts = true warn_unused_configs = true [tool.isort] +# a verbatim copy of the upstream ISMIP7 parameter-selection toolbox; it must +# not be modified, so its imports are not sorted. See PROVENANCE.md alongside +# it. +extend_skip = [ + "compass/landice/tests/ismip7_calibration/toolbox/parameter_selection_toolbox.py", +] multi_line_output = 3 include_trailing_comma = true force_grid_wrap = 0 diff --git a/tests/README.md b/tests/README.md new file mode 100644 index 0000000000..a31653d1f2 --- /dev/null +++ b/tests/README.md @@ -0,0 +1,19 @@ +# Compass unit tests + +Unit tests for compass's own Python code, run with `pytest` from the root of +the repository: + +```bash +pytest +``` + +These are **unit** tests: they must run in seconds, on any machine, with no +input datasets, no MPAS build and no network access. Anything that needs +those belongs in a compass test case instead, where the framework can stage +inputs and compare against a baseline. + +The tests are collected by CI on every pull request, so a test that depends on +data outside the repository will fail there even if it passes locally. + +The layout mirrors the package: `tests/landice/ismip7_calibration/test_terms.py` +tests `compass/landice/tests/ismip7_calibration/terms.py`. diff --git a/tests/landice/ismip7/test_remap.py b/tests/landice/ismip7/test_remap.py new file mode 100644 index 0000000000..8cee3ae923 --- /dev/null +++ b/tests/landice/ismip7/test_remap.py @@ -0,0 +1,131 @@ +""" +Tests for the shared ISMIP7 remapping helpers. +""" + +import logging + +import numpy as np +import xarray as xr + +from compass.landice.ismip7.remap import extrapolate_source + + +def _logger(): + logger = logging.getLogger('test_remap') + logger.addHandler(logging.NullHandler()) + return logger + + +def _write_source(path, values, time_units=None): + """Write a small (time, y, x) source file.""" + ds = xr.Dataset() + ds['field'] = (('time', 'y', 'x'), values) + ds['time'] = ('time', np.arange(values.shape[0])) + if time_units is not None: + ds['time'].attrs['units'] = time_units + ds.to_netcdf(path) + ds.close() + + +def test_extrapolate_fills_all_gaps(tmp_path): + """Every missing value should be replaced by a valid neighbour.""" + values = np.arange(18.0).reshape(2, 3, 3) + values[0, 1, 1] = np.nan + values[1, 0, 0] = np.nan + source = tmp_path / 'source.nc' + output = tmp_path / 'output.nc' + _write_source(source, values) + + extrapolate_source(str(source), str(output), 'field', _logger()) + + with xr.open_dataset(output) as ds: + filled = ds['field'].values + assert np.isfinite(filled).all() + + +def test_extrapolate_takes_the_nearest_value(tmp_path): + """A gap is filled from its nearest valid cell, not by interpolation.""" + values = np.full((1, 1, 5), np.nan) + values[0, 0, 0] = 10.0 + values[0, 0, 4] = 20.0 + source = tmp_path / 'source.nc' + output = tmp_path / 'output.nc' + _write_source(source, values) + + extrapolate_source(str(source), str(output), 'field', _logger()) + + with xr.open_dataset(output) as ds: + filled = ds['field'].values[0, 0] + # nearest-neighbour fill, so every cell takes one of the two valid values + assert set(np.unique(filled)) <= {10.0, 20.0} + assert filled[1] == 10.0 + assert filled[3] == 20.0 + + +def test_extrapolate_leaves_valid_data_untouched(tmp_path): + """Cells that were already valid must not change.""" + values = np.arange(9.0).reshape(1, 3, 3) + values[0, 1, 1] = np.nan + source = tmp_path / 'source.nc' + output = tmp_path / 'output.nc' + _write_source(source, values) + + extrapolate_source(str(source), str(output), 'field', _logger()) + + with xr.open_dataset(output) as ds: + filled = ds['field'].values + valid = np.isfinite(values) + assert np.array_equal(filled[valid], values[valid]) + + +def test_extrapolate_accepts_several_variables(tmp_path): + """A list of variable names extrapolates each of them.""" + ds = xr.Dataset() + first = np.array([[[1.0, np.nan]]]) + second = np.array([[[np.nan, 2.0]]]) + ds['first'] = (('time', 'y', 'x'), first) + ds['second'] = (('time', 'y', 'x'), second) + source = tmp_path / 'source.nc' + output = tmp_path / 'output.nc' + ds.to_netcdf(source) + ds.close() + + extrapolate_source(str(source), str(output), ['first', 'second'], + _logger()) + + with xr.open_dataset(output) as result: + assert np.isfinite(result['first'].values).all() + assert np.isfinite(result['second'].values).all() + + +def test_extrapolate_handles_non_cf_time(tmp_path): + """ + The fracture forcing carries units="year", which xarray cannot decode. + + Those callers pass decode_times=False; the default stays True so that + the ocean and atmosphere callers keep their previous behaviour. + """ + values = np.array([[[1.0, np.nan]]]) + source = tmp_path / 'source.nc' + output = tmp_path / 'output.nc' + _write_source(source, values, time_units='year') + + extrapolate_source(str(source), str(output), 'field', _logger(), + decode_times=False) + + with xr.open_dataset(output, decode_times=False) as ds: + assert np.isfinite(ds['field'].values).all() + assert ds['time'].attrs['units'] == 'year' + + +def test_extrapolate_is_a_no_op_without_gaps(tmp_path): + """Data with no missing values passes through unchanged.""" + values = np.arange(8.0).reshape(2, 2, 2) + source = tmp_path / 'source.nc' + output = tmp_path / 'output.nc' + _write_source(source, values) + + extrapolate_source(str(source), str(output), 'field', _logger()) + + with xr.open_dataset(output) as ds: + assert np.array_equal(ds['field'].values, values) diff --git a/tests/landice/ismip7_calibration/test_configure.py b/tests/landice/ismip7_calibration/test_configure.py new file mode 100644 index 0000000000..e986f7cdc3 --- /dev/null +++ b/tests/landice/ismip7_calibration/test_configure.py @@ -0,0 +1,162 @@ +""" +Tests for the ISMIP7 calibration config helpers. +""" + +import numpy as np +import pytest + +from compass.config import CompassConfigParser +from compass.landice.tests.ismip7_calibration.configure import ( + check_options, + melt_forms, + objective_options, + parameter_name, + parameter_values, + weighting, +) + + +def _config(): + """The test group's default config options.""" + config = CompassConfigParser() + config.add_from_package('compass.landice.tests.ismip7_calibration', + 'ismip7_calibration.cfg') + return config + + +def test_check_options_rejects_unset_paths(): + """A path the user must supply should fail loudly, not silently.""" + config = _config() + + with pytest.raises(ValueError, match='base_path_ismip7'): + check_options(config, ['base_path_ismip7']) + + +def test_check_options_accepts_a_supplied_path(): + """Once set, the same option passes.""" + config = _config() + config.set('ismip7_calibration', 'base_path_ismip7', '/some/path') + + check_options(config, ['base_path_ismip7']) + + +def test_parameter_grid_matches_the_published_k_values(): + """ + The default K grid must be the 120 values the protocol's worked example + samples, or the replicated percentiles land on different grid points. + """ + values = parameter_values(_config(), 'ismip7') + + assert values[0] == pytest.approx(0.25e-5) + assert values[-1] == pytest.approx(3.0e-4) + assert len(values) == 120 + np.testing.assert_allclose(np.diff(values), 0.25e-5, rtol=1.0e-9) + + +def test_parameter_grid_for_the_nonlocal_form(): + """ + gamma0 has its own grid, in m/yr. The maximum has to stay well above + where the distribution ends, or the objective's draws pile up against + it and the upper percentiles become a property of the grid. + """ + values = parameter_values(_config(), 'ismip6') + + assert values[0] == pytest.approx(250.0) + assert values[-1] == pytest.approx(50000.0) + + +def test_parameter_values_rejects_an_unknown_form(): + with pytest.raises(ValueError, match="must be 'ismip7' or 'ismip6'"): + parameter_values(_config(), 'pico') + + +def test_parameter_name_per_form(): + """The two forms calibrate different parameters.""" + assert parameter_name('ismip7') == 'K' + assert parameter_name('ismip6') == 'gamma0' + with pytest.raises(ValueError): + parameter_name('pico') + + +def test_melt_forms_parses_a_comma_separated_list(): + config = _config() + config.set('ismip7_calibration', 'melt_forms', 'ismip7, ismip6') + + assert melt_forms(config) == ['ismip7', 'ismip6'] + + +def test_melt_forms_accepts_a_single_form(): + config = _config() + config.set('ismip7_calibration', 'melt_forms', 'ismip7') + + assert melt_forms(config) == ['ismip7'] + + +def test_melt_forms_rejects_an_unknown_form(): + config = _config() + config.set('ismip7_calibration', 'melt_forms', 'ismip7, pico') + + with pytest.raises(ValueError, match='pico'): + melt_forms(config) + + +def test_melt_forms_rejects_an_empty_list(): + config = _config() + config.set('ismip7_calibration', 'melt_forms', '') + + with pytest.raises(ValueError, match='at least one'): + melt_forms(config) + + +def test_weighting_all_gives_no_restriction(): + """'all' means every summand the ensemble provides carries weight.""" + config = _config() + config.set('ismip7_calibration_objective', 't3_models', 'all') + config.set('ismip7_calibration_objective', 't4_regions', 'all') + config.set('ismip7_calibration_objective', 't4_years', 'all') + + assert weighting(config) == (None, None, None) + + +def test_weighting_published_restricts_j4_to_pig(): + """ + The published weighting uses PIG alone in 2009 and 2012 -- 2 of the 18 + available observations. That is a deliberate option, not the default, + because so weighted J4 cannot discriminate between melt forms. + """ + config = _config() + config.set('ismip7_calibration_objective', 't3_models', 'published') + config.set('ismip7_calibration_objective', 't4_regions', 'published') + config.set('ismip7_calibration_objective', 't4_years', 'published') + + t3_models, t4_regions, t4_years = weighting(config) + + assert t4_regions == ('pig',) + assert t4_years == (2009, 2012) + assert len(t3_models) == 4 + + +def test_weighting_rejects_an_unknown_choice(): + config = _config() + config.set('ismip7_calibration_objective', 't3_models', 'some') + + with pytest.raises(ValueError, match="must be 'published' or 'all'"): + weighting(config) + + +def test_objective_options_defaults_are_reproducible(): + """A fixed seed is the default, so the percentiles do not drift.""" + sample_size, seed = objective_options(_config()) + + assert sample_size == 100000 + assert seed == 0 + + +def test_objective_options_allows_no_seed(): + """'None' leaves the global random state alone.""" + config = _config() + config.set('ismip7_calibration_objective', 'seed', 'None') + + _, seed = objective_options(config) + + assert seed is None diff --git a/tests/landice/ismip7_calibration/test_datasets.py b/tests/landice/ismip7_calibration/test_datasets.py new file mode 100644 index 0000000000..abc3f3b9e3 --- /dev/null +++ b/tests/landice/ismip7_calibration/test_datasets.py @@ -0,0 +1,148 @@ +""" +Tests for the ISMIP7 ocean-state and dataset registry. + +These check the registry's internal consistency, which is what keeps the +8 km replication and the MALI-mesh calibration driving the same states. They +do not touch the ISMIP7 datasets themselves, which are not in the repository. +""" + +import pytest + +from compass.landice.tests.ismip7_calibration import datasets + +BASE_PATH = '/not/a/real/path' + + +def test_subset_sizes_match_the_protocol(): + """ + 7 states for the minimal set, 11 for the recommended one and 28 for + everything ISMIP7 distributes. + """ + assert len(datasets.ocean_states(BASE_PATH, 'minimal')) == 7 + assert len(datasets.ocean_states(BASE_PATH, 'recommended')) == 11 + assert len(datasets.ocean_states(BASE_PATH, 'all')) == 28 + + +def test_subsets_are_nested(): + """A larger subset must contain everything a smaller one does.""" + minimal = {state.name for state in + datasets.ocean_states(BASE_PATH, 'minimal')} + recommended = {state.name for state in + datasets.ocean_states(BASE_PATH, 'recommended')} + every = {state.name for state in datasets.ocean_states(BASE_PATH, 'all')} + + assert minimal <= recommended <= every + + +def test_an_unknown_subset_is_rejected(): + with pytest.raises(ValueError, match="must be 'minimal'"): + datasets.ocean_states(BASE_PATH, 'everything') + + +def test_state_names_are_unique(): + """Names become directory names, so a clash would silently overwrite.""" + names = [state.name for state in datasets.ocean_states(BASE_PATH, 'all')] + + assert len(names) == len(set(names)) + + +def test_every_state_feeds_a_term(): + """A state that constrains nothing would just waste a MALI run.""" + for state in datasets.ocean_states(BASE_PATH, 'all'): + assert state.term in ('J1,J2', 'J3', 'J4') + + +def test_the_climatology_is_the_only_present_day_state(): + """J1 and J2 are present-day terms, driven by the one climatology.""" + states = datasets.ocean_states(BASE_PATH, 'all') + climatology = [state for state in states if state.kind == 'climatology'] + + assert len(climatology) == 1 + assert climatology[0].name == 'climatology' + assert climatology[0].term == 'J1,J2' + + +def test_ocean_models_come_in_cold_and_warm_pairs(): + """J3 is a warm-minus-cold difference, so both must be present.""" + states = [state for state in datasets.ocean_states(BASE_PATH, 'all') + if state.kind == 'model'] + by_label = {} + for state in states: + by_label.setdefault(state.label, set()).add(state.state) + + assert len(by_label) == 7 + for label, which in by_label.items(): + assert which == {'cold', 'warm'}, label + + +def test_regional_models_declare_the_basins_they_cover(): + """ + Melt outside a regional model's domain must be discarded, so the + covered basins have to be recorded. Basin numbers are the ISMIP7 + 0-based convention. + """ + states = {state.label: state + for state in datasets.ocean_states(BASE_PATH, 'all') + if state.kind == 'model'} + + assert states['mathiot'].basins is None + assert states['naughten_ais_1'].basins is None + assert states['jourdain_naughten'].basins == (9, 14) + assert states['naughten_naughten'].basins == (9, 14) + assert states['timmermann'].basins == (14,) + + +def test_observations_constrain_the_eastern_amundsen(): + """ + PIG and Dotson are both in ISMIP7 basin 9, the Eastern Amundsen, in the + 0-based convention the protocol uses. + """ + states = [state for state in datasets.ocean_states(BASE_PATH, 'all') + if state.kind == 'obs'] + + assert len(states) == 13 + for state in states: + assert state.basins == (9,) + assert state.year in datasets.OBS_YEARS + + +def test_the_published_weighting_is_a_strict_subset(): + """ + The published J4 weighting is PIG in 2009 and 2012 only -- 2 of the 18 + available observations. Pinning it here guards the replication. + """ + assert datasets.PUBLISHED_T4_REGIONS == ('pig',) + assert datasets.PUBLISHED_T4_YEARS == (2009, 2012) + assert set(datasets.PUBLISHED_T4_YEARS) < set(datasets.OBS_YEARS) + assert len(datasets.PUBLISHED_T3_MODELS) == 4 + + +def test_mask_file_names_follow_the_resolution(): + """The ISMIP masks are distributed at several resolutions.""" + files = datasets.mask_files(BASE_PATH, resolution_km=8) + + assert files['basins'].endswith('basin_numbers_ismip8km_v2.nc') + assert files['bfrn'].endswith('BFRN_ismip8km_v2.nc') + assert files['floating'].endswith('floatingmask_ismip8km.nc') + # the shelf mask is only distributed at 8 km + assert files['shelves'].endswith('shelf_mask_ismip8km.nc') + + +def test_shelf_ids_and_the_pine_island_cut(): + """ + The worked example keeps only Pine Island's main trunk, cutting it at a + fixed x. Getting the ids or the cut wrong would silently change J4. + """ + assert datasets.PIG_ID == 110 + assert datasets.DOTSON_ID == 97 + assert datasets.PIG_X_MAX == -1.625e6 + + +def test_missing_files_reports_absent_inputs(): + """The registry points at a fake root, so everything is missing.""" + states = datasets.ocean_states(BASE_PATH, 'minimal') + + missing = datasets.missing_files(states) + + # both the thermal-forcing and the salinity file of each state + assert len(missing) == 2 * len(states) diff --git a/tests/landice/ismip7_calibration/test_melt_model.py b/tests/landice/ismip7_calibration/test_melt_model.py new file mode 100644 index 0000000000..3d3c2be815 --- /dev/null +++ b/tests/landice/ismip7_calibration/test_melt_model.py @@ -0,0 +1,331 @@ +""" +Tests for recomputing MALI's melt from its diagnostic output. + +The vertical interpolation here is written from the protocol rather than +transliterated from MALI's Fortran, so that comparing the two in the +``verify_melt`` step is a real check. These tests pin its four code paths +and the area-weighted basin mean. +""" + +import numpy as np +import pytest +import xarray as xr + +from compass.landice.tests.ismip7_calibration.ais.melt_model import ( + FREEZING_TEMP_DEPTH_DEPENDENCE, + basin_mean_tf, + initial_draft, + integrate_by_basin, + interpolate_to_draft, + read_run, +) +from compass.landice.tests.ismip7_calibration.quadratic import MALI + +#: three layer centres, negative downward as MALI expects +Z_OCEAN = np.array([-30.0, -90.0, -150.0]) + + +def _field(values): + """One cell with a value per layer.""" + return np.array([values], dtype=float) + + +def test_interpolation_above_the_shallowest_centre(): + """A draft above the top layer centre takes the top layer's value.""" + result = interpolate_to_draft(_field([1.0, 2.0, 3.0]), Z_OCEAN, + np.array([-10.0]), np.array([-500.0])) + + assert result[0] == pytest.approx(1.0) + + +def test_interpolation_below_the_deepest_centre(): + """A draft below the bottom layer centre takes the bottom value.""" + result = interpolate_to_draft(_field([1.0, 2.0, 3.0]), Z_OCEAN, + np.array([-300.0]), np.array([-500.0])) + + assert result[0] == pytest.approx(3.0) + + +def test_interpolation_when_the_layer_below_is_beneath_the_bed(): + """ + Where the next layer down is below the bed there is no water there, so + the shallower layer is used rather than interpolating into rock. + """ + # the draft is at -60 and the bed at -80, so the -90 layer centre is + # below the bed and there is no water there to interpolate into + result = interpolate_to_draft(_field([1.0, 2.0, 3.0]), Z_OCEAN, + np.array([-60.0]), np.array([-80.0])) + + assert result[0] == pytest.approx(1.0) + + +def test_interpolation_between_layer_centres(): + """Half way between two centres gives the average of their values.""" + result = interpolate_to_draft(_field([1.0, 2.0, 3.0]), Z_OCEAN, + np.array([-60.0]), np.array([-500.0])) + + assert result[0] == pytest.approx(1.5) + + +def test_interpolation_is_linear_in_depth(): + """A quarter of the way down gives a quarter of the difference.""" + result = interpolate_to_draft(_field([0.0, 4.0, 8.0]), Z_OCEAN, + np.array([-45.0]), np.array([-500.0])) + + assert result[0] == pytest.approx(1.0) + + +def test_interpolation_reproduces_layer_centre_values_exactly(): + """At a layer centre the result must be that layer's value.""" + values = [1.0, 2.0, 3.0] + for index, depth in enumerate(Z_OCEAN): + result = interpolate_to_draft(_field(values), Z_OCEAN, + np.array([depth]), + np.array([-500.0])) + assert result[0] == pytest.approx(values[index]) + + +def test_interpolation_handles_many_cells_independently(): + """Each cell uses its own draft and bed.""" + field = np.array([[1.0, 2.0, 3.0], [1.0, 2.0, 3.0]]) + result = interpolate_to_draft(field, Z_OCEAN, + np.array([-10.0, -300.0]), + np.array([-500.0, -500.0])) + + assert result[0] == pytest.approx(1.0) + assert result[1] == pytest.approx(3.0) + + +def test_freezing_correction_applies_below_the_deepest_centre(): + """ + MALI corrects the thermal forcing for the depth dependence of the + freezing point where it extrapolates downward from the deepest layer. + """ + draft = -300.0 + result = interpolate_to_draft(_field([1.0, 2.0, 3.0]), Z_OCEAN, + np.array([draft]), np.array([-500.0]), + freezing_correction=True) + + expected = 3.0 - (Z_OCEAN[-1] - draft) * FREEZING_TEMP_DEPTH_DEPENDENCE + assert result[0] == pytest.approx(expected) + + +def test_freezing_correction_applies_below_the_bed(): + """ + The same correction applies in the other downward-extrapolating branch, + where the layer below the draft is beneath the bed. Omitting it here + was a real bug the verification step caught. + """ + draft = -60.0 + result = interpolate_to_draft(_field([1.0, 2.0, 3.0]), Z_OCEAN, + np.array([draft]), np.array([-80.0]), + freezing_correction=True) + + expected = 1.0 - (Z_OCEAN[0] - draft) * FREEZING_TEMP_DEPTH_DEPENDENCE + assert result[0] == pytest.approx(expected) + + +def test_freezing_correction_does_not_apply_when_interpolating(): + """ + Between two layer centres there is no extrapolation, so the correction + must not be applied. + """ + plain = interpolate_to_draft(_field([1.0, 2.0, 3.0]), Z_OCEAN, + np.array([-60.0]), np.array([-500.0])) + corrected = interpolate_to_draft(_field([1.0, 2.0, 3.0]), Z_OCEAN, + np.array([-60.0]), np.array([-500.0]), + freezing_correction=True) + + assert plain[0] == pytest.approx(corrected[0]) + + +def test_freezing_correction_does_not_apply_above_the_shallowest_centre(): + """The draft is above the layer, so there is nothing to correct.""" + plain = interpolate_to_draft(_field([1.0, 2.0, 3.0]), Z_OCEAN, + np.array([-10.0]), np.array([-500.0])) + corrected = interpolate_to_draft(_field([1.0, 2.0, 3.0]), Z_OCEAN, + np.array([-10.0]), np.array([-500.0]), + freezing_correction=True) + + assert plain[0] == pytest.approx(corrected[0]) + + +def test_basin_mean_is_area_weighted(): + """ + MALI's non-local form uses an area-weighted basin mean; a plain mean + would be wrong on a variable-resolution mesh. + """ + tf = xr.DataArray([0.0, 10.0], dims='nCells') + area = xr.DataArray([9.0, 1.0], dims='nCells') + floating = xr.DataArray([True, True], dims='nCells') + basin = xr.DataArray([1, 1], dims='nCells') + + result = basin_mean_tf(tf, area, floating, basin) + + np.testing.assert_allclose(result.values, [1.0, 1.0]) + + +def test_basin_mean_excludes_grounded_cells(): + """Only floating cells contribute to the mean.""" + tf = xr.DataArray([100.0, 2.0], dims='nCells') + area = xr.DataArray([1.0, 1.0], dims='nCells') + floating = xr.DataArray([False, True], dims='nCells') + basin = xr.DataArray([1, 1], dims='nCells') + + result = basin_mean_tf(tf, area, floating, basin) + + np.testing.assert_allclose(result.values, [2.0, 2.0]) + + +def test_basin_mean_keeps_basins_separate(): + """Each cell gets the mean of its own basin.""" + tf = xr.DataArray([1.0, 3.0, 10.0], dims='nCells') + area = xr.DataArray([1.0, 1.0, 1.0], dims='nCells') + floating = xr.DataArray([True, True, True], dims='nCells') + basin = xr.DataArray([1, 1, 2], dims='nCells') + + result = basin_mean_tf(tf, area, floating, basin) + + np.testing.assert_allclose(result.values, [2.0, 2.0, 10.0]) + + +def test_integrate_by_basin_converts_to_gigatonnes(): + """kg m-2 yr-1 times m2 is kg yr-1; the result must be Gt yr-1.""" + melt = xr.DataArray([1.0], dims='nCells') + area = xr.DataArray([1.0e12], dims='nCells') + floating = xr.DataArray([True], dims='nCells') + basin = xr.DataArray([3], dims='nCells') + + result = integrate_by_basin(melt, area, floating, basin) + + assert float(result.sel(basin=3)) == pytest.approx(1.0) + + +def test_initial_draft_uses_flotation_where_floating(tmp_path): + """ + The draft must be reconstructed from the mesh, not from the run output, + because melt thins the ice over the timestep. + """ + path = tmp_path / 'mesh.nc' + ds = xr.Dataset() + ds['thickness'] = ('nCells', np.array([1000.0])) + ds['bedTopography'] = ('nCells', np.array([-2000.0])) + ds.to_netcdf(path) + ds.close() + + draft = initial_draft(str(path)) + + expected = -MALI.rho_ice / MALI.rho_ocean * 1000.0 + assert float(draft[0]) == pytest.approx(expected) + + +def test_initial_draft_uses_the_bed_where_grounded(tmp_path): + """Grounded ice sits on the bed, above its flotation draft.""" + path = tmp_path / 'mesh.nc' + ds = xr.Dataset() + ds['thickness'] = ('nCells', np.array([100.0])) + ds['bedTopography'] = ('nCells', np.array([-50.0])) + ds.to_netcdf(path) + ds.close() + + draft = initial_draft(str(path)) + + assert float(draft[0]) == pytest.approx(-50.0) + + +def test_initial_draft_handles_a_time_dimension(tmp_path): + """MALI mesh files may carry a singleton Time dimension.""" + path = tmp_path / 'mesh.nc' + ds = xr.Dataset() + ds['thickness'] = (('Time', 'nCells'), np.array([[1000.0]])) + ds['bedTopography'] = (('Time', 'nCells'), np.array([[-2000.0]])) + ds.to_netcdf(path) + ds.close() + + draft = initial_draft(str(path)) + + assert draft.dims == ('nCells',) + + +def _write_run(tmp_path, thickness, bed, connected, name='output_melt.nc'): + """A minimal melt-diagnostic output and its mesh file.""" + mesh = tmp_path / 'mesh.nc' + ds_mesh = xr.Dataset() + ds_mesh['thickness'] = (('Time', 'nCells'), np.array([thickness])) + ds_mesh['bedTopography'] = (('Time', 'nCells'), np.array([bed])) + ds_mesh.to_netcdf(mesh) + ds_mesh.close() + + n_cells = len(thickness) + run = tmp_path / name + ds = xr.Dataset() + ds['floatingBasalMassBal'] = (('Time', 'nCells'), + np.zeros((1, n_cells))) + ds['ismip6shelfMelt_TFdraft'] = (('Time', 'nCells'), + np.zeros((1, n_cells))) + ds['connectedOceanMask'] = (('Time', 'nCells'), + np.array([connected], dtype=np.int32)) + ds['ismip6shelfMelt_basin'] = ('nCells', + np.ones(n_cells, dtype=np.int32)) + ds['ismip6shelfMelt_deltaT'] = ('nCells', np.zeros(n_cells)) + ds['areaCell'] = ('nCells', np.ones(n_cells)) + ds.to_netcdf(run) + ds.close() + return str(run), str(mesh) + + +def test_contributing_cells_come_from_the_initial_geometry(tmp_path): + """ + Melt is computed on the geometry at the *start* of the timestep, so a + cell with no ice then must not contribute, however the post-step output + looks. Over one step some ice-free cells pick up a trace of ice, and + counting them would dilute the area-weighted basin means. + """ + run, mesh = _write_run(tmp_path, + thickness=[0.0, 1000.0], + bed=[-500.0, -2000.0], + connected=[1, 1]) + + fields = read_run(run, mesh) + + assert not bool(fields['floating'][0]) + assert bool(fields['floating'][1]) + + +def test_grounded_ice_does_not_contribute(tmp_path): + """Thick ice on a shallow bed is grounded, so it has no shelf melt.""" + run, mesh = _write_run(tmp_path, + thickness=[1000.0], + bed=[-100.0], + connected=[1]) + + fields = read_run(run, mesh) + + assert not bool(fields['floating'][0]) + + +def test_cells_cut_off_from_the_ocean_do_not_contribute(tmp_path): + """ + MALI computes melt only where the cell is connected to the open ocean, + leaving the melt at zero elsewhere. + """ + run, mesh = _write_run(tmp_path, + thickness=[1000.0], + bed=[-2000.0], + connected=[0]) + + fields = read_run(run, mesh) + + assert not bool(fields['floating'][0]) + + +def test_ice_exactly_at_flotation_counts_as_floating(tmp_path): + """MALI treats ice exactly at flotation as floating.""" + thickness = 1000.0 + bed = -MALI.rho_ice / MALI.rho_ocean * thickness + run, mesh = _write_run(tmp_path, thickness=[thickness], bed=[bed], + connected=[1]) + + fields = read_run(run, mesh) + + assert bool(fields['floating'][0]) diff --git a/tests/landice/ismip7_calibration/test_objective.py b/tests/landice/ismip7_calibration/test_objective.py new file mode 100644 index 0000000000..dd3822768a --- /dev/null +++ b/tests/landice/ismip7_calibration/test_objective.py @@ -0,0 +1,244 @@ +""" +Tests for assembling the ISMIP7 objective function's arguments. + +The weighting decides which observations constrain the calibration at all, +so getting it wrong changes the answer silently. These tests pin the +scaling that the one-run-per-state ensemble design relies on, and the +weighting rules. +""" + +import numpy as np +import pytest +import xarray as xr + +from compass.landice.tests.ismip7_calibration.objective import ( + build_toolbox_terms, + scale_to_ensemble, +) + +PARAMETERS = np.array([1.0, 2.0, 4.0]) +BASINS = np.arange(3) +MODELS = ['mathiot', 'haid'] +REGIONS = ['pig', 'dotson'] +YEARS = [2009, 2012, 2014] + + +def _units(): + """Unit-parameter aggregates with the shapes the toolbox expects.""" + t1 = xr.DataArray(np.ones(3), dims='basins', coords={'basins': BASINS}) + t2 = xr.DataArray(np.ones(4), dims='BFRN_bins', + coords={'BFRN_bins': np.arange(4)}) + t3 = xr.DataArray(np.ones((2, 3)), dims=('model', 'basins'), + coords={'model': MODELS, 'basins': BASINS}) + t4 = xr.DataArray(np.ones((3, 2)), dims=('year', 'region'), + coords={'year': YEARS, 'region': REGIONS}) + return dict(t1=t1, t2=t2, t3=t3, t4=t4) + + +def _targets(): + """Targets matching the aggregate shapes.""" + t4_mean = xr.DataArray(np.ones((2, 3)), dims=('region', 'year'), + coords={'region': REGIONS, 'year': YEARS}) + return dict( + t1_mean=xr.DataArray(np.ones(3), dims='basin'), + t1_sigma=xr.DataArray(np.ones(3), dims='basin'), + t2_mean=xr.DataArray(np.ones(4), dims='BFRN_bins'), + t2_sigma=xr.DataArray(np.ones(4), dims='BFRN_bins'), + t2_weights=xr.DataArray(np.ones(4), dims='BFRN_bins'), + t3_mean=xr.DataArray(np.ones((2, 3)), dims=('model', 'basins'), + coords={'model': MODELS, 'basins': BASINS}), + t3_sigma=xr.DataArray(np.ones((2, 3)), dims=('model', 'basins'), + coords={'model': MODELS, 'basins': BASINS}), + t4_mean=t4_mean, + t4_sigma=t4_mean.copy()) + + +def test_scaling_is_exact_and_proportional(): + """ + The ensemble is one run per ocean state scaled onto the parameter grid. + That is exact because melt is proportional to the parameter, so the + scaling must be too. + """ + unit = xr.DataArray([3.0, 5.0], dims='basins') + + scaled = scale_to_ensemble(unit, PARAMETERS) + + for index, value in enumerate(PARAMETERS): + np.testing.assert_allclose( + scaled.isel(p1=index, p2=0).values, value * unit.values) + + +def test_scaling_adds_the_singleton_second_parameter(): + """ + The toolbox indexes terms by (p1, p2); the quadratic forms use only one + parameter, so p2 is a singleton. + """ + unit = xr.DataArray([1.0, 2.0], dims='basins') + + scaled = scale_to_ensemble(unit, PARAMETERS) + + assert scaled.dims[:2] == ('p1', 'p2') + assert scaled.sizes['p2'] == 1 + assert scaled.sizes['p1'] == len(PARAMETERS) + + +def test_all_summands_are_weighted_by_default(): + """Passing None weights everything the ensemble provides.""" + terms = build_toolbox_terms(_units(), _targets(), PARAMETERS) + + assert (terms['t3_weights'] == 1).all() + assert (terms['t4_weights'] == 1).all() + + +def test_restricting_j3_to_named_models(): + """A model left out of the weighting must not constrain the fit.""" + terms = build_toolbox_terms(_units(), _targets(), PARAMETERS, + t3_models=('mathiot',)) + + weights = terms['t3_weights'] + assert (weights.sel(model='mathiot') == 1).all() + assert (weights.sel(model='haid') == 0).all() + + +def test_restricting_j4_to_pig_zeroes_dotson(): + """ + The published weighting excludes Dotson entirely. With one shelf J4 + constrains a single amplitude that either melt form can match by + rescaling its parameter, so it stops discriminating between them -- + which is exactly why this is an option and not the default. + """ + terms = build_toolbox_terms(_units(), _targets(), PARAMETERS, + t4_regions=('pig',)) + + weights = terms['t4_weights'] + assert (weights.sel(region='pig') == 1).all() + assert (weights.sel(region='dotson') == 0).all() + + +def test_restricting_j4_to_named_years(): + terms = build_toolbox_terms(_units(), _targets(), PARAMETERS, + t4_years=(2009, 2012)) + + weights = terms['t4_weights'] + assert float(weights.sel(region='pig', year=2009)) == 1.0 + assert float(weights.sel(region='pig', year=2014)) == 0.0 + + +def test_a_year_the_ensemble_did_not_run_carries_no_weight(): + """ + Asking for a year with no MALI run must not silently contribute; the + weighting is intersected with what the ensemble actually provides. + """ + units = _units() + units['t4'] = units['t4'].sel(year=[2009, 2012]) + targets = _targets() + + terms = build_toolbox_terms(units, targets, PARAMETERS, t4_years=None) + + assert float(terms['t4_weights'].sel(region='pig', year=2014)) == 0.0 + + +def test_published_and_full_weightings_differ(): + """ + The two weightings must not be accidentally identical, or the config + option would be meaningless. + """ + full = build_toolbox_terms(_units(), _targets(), PARAMETERS) + published = build_toolbox_terms(_units(), _targets(), PARAMETERS, + t3_models=('mathiot',), + t4_regions=('pig',), + t4_years=(2009, 2012)) + + assert float(full['t4_weights'].sum()) > \ + float(published['t4_weights'].sum()) + assert float(full['t3_weights'].sum()) > \ + float(published['t3_weights'].sum()) + + +def test_the_toolbox_receives_every_argument_it_needs(): + """A missing key would surface as an opaque TypeError at call time.""" + terms = build_toolbox_terms(_units(), _targets(), PARAMETERS) + + for term in ('t1', 't2', 't3', 't4'): + for suffix in ('model', 'obs_mean', 'obs_sigma', 'weights'): + assert f'{term}_{suffix}' in terms + + +def test_j3_targets_are_aligned_with_the_ensemble_models(): + """ + The targets carry every ocean model; selecting the ensemble's models + keeps the two aligned rather than broadcasting against each other. + """ + terms = build_toolbox_terms(_units(), _targets(), PARAMETERS) + + assert list(terms['t3_obs_mean'].model.values) == \ + list(terms['t3_model'].model.values) + + +def test_scaling_preserves_nan_for_empty_groups(): + """ + An aggregate that came out NaN because no cells contributed must stay + NaN through the scaling, so it drops out of the objective. + """ + unit = xr.DataArray([1.0, np.nan], dims='basins') + + scaled = scale_to_ensemble(unit, PARAMETERS) + + assert np.isnan(scaled.isel(p1=0, p2=0).values[1]) + + +def test_scaling_by_zero_gives_zero_melt(): + """A zero parameter must give zero aggregate, not NaN.""" + unit = xr.DataArray([2.0], dims='basins') + + scaled = scale_to_ensemble(unit, np.array([0.0, 1.0])) + + assert float(scaled.isel(p1=0, p2=0)[0]) == pytest.approx(0.0) + + +def _fake_toolbox_result(monkeypatch, minimisers): + """Make calculate_objective_function return chosen minimisers.""" + from compass.landice.tests.ismip7_calibration import objective + + def fake(*args, **kwargs): + return np.asarray(minimisers, dtype=float), None + + monkeypatch.setattr(objective.toolbox, 'calculate_objective_function', + fake) + + +def _placeholder_terms(): + """The keys run_optimisation passes through to the toolbox.""" + return {f'{term}_{suffix}': None + for term in ('t1', 't2', 't3', 't4') + for suffix in ('model', 'obs_mean', 'obs_sigma', 'weights')} + + +def test_no_draws_at_the_top_of_a_wide_enough_grid(monkeypatch): + """A distribution well inside the grid is not clipped.""" + values = np.array([1.0, 2.0, 3.0, 4.0]) + _fake_toolbox_result(monkeypatch, [1.0, 2.0, 2.0, 3.0]) + + from compass.landice.tests.ismip7_calibration.objective import ( + run_optimisation, + ) + result = run_optimisation(_placeholder_terms(), values, sample_size=4) + + assert result['at_upper_bound'] == pytest.approx(0.0) + + +def test_draws_on_the_largest_parameter_are_counted(monkeypatch): + """ + A draw that picked the largest parameter wanted a larger one, so the + grid has clipped it. Counting them is what makes the clipping visible + instead of silent. + """ + values = np.array([1.0, 2.0, 3.0, 4.0]) + _fake_toolbox_result(monkeypatch, [1.0, 4.0, 4.0, 2.0]) + + from compass.landice.tests.ismip7_calibration.objective import ( + run_optimisation, + ) + result = run_optimisation(_placeholder_terms(), values, sample_size=4) + + assert result['at_upper_bound'] == pytest.approx(0.5) diff --git a/tests/landice/ismip7_calibration/test_quadratic.py b/tests/landice/ismip7_calibration/test_quadratic.py new file mode 100644 index 0000000000..b2562f8f18 --- /dev/null +++ b/tests/landice/ismip7_calibration/test_quadratic.py @@ -0,0 +1,184 @@ +""" +Tests for the Python reference implementations of the melt formulas. + +These are the specification MALI's Fortran is checked against, so the +properties the calibration depends on are asserted here directly. +""" + +import numpy as np +import pytest +import xarray as xr + +from compass.landice.tests.ismip7_calibration.quadratic import ( + MALI, + MULTIMELT, + angle_from_sin_slope, + local_quadratic_melt, + nonlocal_quadratic_melt, +) + +SALINITY = 34.5 +SLOPE = angle_from_sin_slope(0.0051117) + + +def _thermal_forcing(values): + return xr.DataArray(np.asarray(values, dtype=float), dims='nCells') + + +def test_angle_from_sin_slope_round_trips(): + """ + MALI's config option is a sine while local_quadratic_melt takes an + angle; mixing them would silently rescale the melt. + """ + for sin_slope in (0.0051117, 0.005, 0.1): + assert np.sin(angle_from_sin_slope(sin_slope)) == \ + pytest.approx(sin_slope) + + +def test_local_melt_is_exactly_linear_in_k(): + """ + The whole ensemble design rests on this: one MALI run per ocean state, + scaled onto the parameter grid, instead of one run per (state, K) pair. + """ + tf = _thermal_forcing([0.5, 1.0, 2.0, -0.5]) + base = local_quadratic_melt(1.0, tf, SALINITY, SLOPE, constants=MALI) + + for factor in (0.25, 3.0, 17.0): + scaled = local_quadratic_melt(factor, tf, SALINITY, SLOPE, + constants=MALI) + np.testing.assert_allclose(scaled.values, factor * base.values, + rtol=1.0e-14, atol=0.0) + + +def test_nonlocal_melt_is_exactly_linear_in_gamma0(): + """The same linearity must hold for the ISMIP6 non-local form.""" + tf = _thermal_forcing([0.5, 1.0, 2.0, -0.5]) + mean_tf = _thermal_forcing([1.0, 1.0, 1.0, 1.0]) + base = nonlocal_quadratic_melt(1.0, tf, mean_tf, constants=MALI) + + for factor in (0.25, 3.0, 17.0): + scaled = nonlocal_quadratic_melt(factor, tf, mean_tf, constants=MALI) + np.testing.assert_allclose(scaled.values, factor * base.values, + rtol=1.0e-14, atol=0.0) + + +def test_local_melt_is_quadratic_in_thermal_forcing(): + """Doubling the thermal forcing must quadruple the melt.""" + single = local_quadratic_melt(1.0, _thermal_forcing([1.0]), SALINITY, + SLOPE, constants=MALI).item() + double = local_quadratic_melt(1.0, _thermal_forcing([2.0]), SALINITY, + SLOPE, constants=MALI).item() + assert double == pytest.approx(4.0 * single) + + +def test_melt_keeps_the_sign_of_the_thermal_forcing(): + """ + The |TF| TF factor means warm water melts and cold water refreezes; a + sign error here would be invisible in a total but wrong per basin. + """ + tf = _thermal_forcing([2.0, -2.0]) + melt = local_quadratic_melt(1.0, tf, SALINITY, SLOPE, constants=MALI) + assert melt[0].item() > 0.0 + assert melt[1].item() < 0.0 + assert melt[0].item() == pytest.approx(-melt[1].item()) + + +def test_local_melt_is_linear_in_salinity(): + """ + Melt is linear in salinity, which is what bounds the effect of using a + constant value rather than a field. + """ + tf = _thermal_forcing([1.0]) + low = local_quadratic_melt(1.0, tf, 34.0, SLOPE, + constants=MALI).item() + high = local_quadratic_melt(1.0, tf, 34.0 * 2.0, SLOPE, + constants=MALI).item() + assert high == pytest.approx(2.0 * low) + + +def test_delta_t_shifts_both_thermal_forcing_factors(): + """ + Protocol Sect. 4.2.1 applies the basin correction wherever the thermal + forcing appears, so melt at TF with dT must equal melt at TF + dT with + no correction. + """ + tf = _thermal_forcing([1.0, 2.0]) + with_correction = local_quadratic_melt(1.0, tf, SALINITY, SLOPE, + constants=MALI, delta_t=0.5) + shifted = local_quadratic_melt(1.0, tf + 0.5, SALINITY, SLOPE, + constants=MALI) + np.testing.assert_allclose(with_correction.values, shifted.values, + rtol=1.0e-14) + + +def test_semi_local_is_degenerate_with_the_nonlocal_form(): + """ + With a constant salinity the Burgard semi-local form is algebraically + identical to the ISMIP6 non-local method MALI already has, so every K + has an exactly equivalent gamma0. That degeneracy is why the + semi-local form is calibrated through the ISMIP6 code path rather than + added as a third melt module, so it is worth asserting. + """ + tf = _thermal_forcing([0.5, 1.5, -0.5]) + mean_tf = _thermal_forcing([1.0, 1.0, 1.0]) + + melt_k = 8.5e-5 + # semi-local: the local form with the basin mean in the |TF| factor + semi_local = local_quadratic_melt(melt_k, tf, SALINITY, SLOPE, + thermal_forcing_avg=mean_tf, + constants=MALI) + + # the equivalent gamma0 follows from equating the two constants + cste = (MALI.rho_ocean * MALI.c_o / + (MALI.rho_ice * MALI.latent_heat))**2 + # gamma0 is already a per-year velocity scale, so the conversion + # carries the year length that the local form applies explicitly + coefficient = (np.sin(SLOPE) * (MALI.rho_ocean / MALI.rho_ice) * + (MALI.c_o / MALI.latent_heat)**2 * MALI.beta_s * + MALI.gravity / (2.0 * abs(MALI.coriolis)) * SALINITY * + MALI.seconds_per_year) + gamma0 = melt_k * coefficient / cste + + non_local = nonlocal_quadratic_melt(gamma0, tf, mean_tf, constants=MALI) + + np.testing.assert_allclose(semi_local.values, non_local.values, + rtol=1.0e-12) + + +def test_local_and_semi_local_differ_in_pattern(): + """ + Replacing the basin mean with the local forcing changes the spatial + pattern in a way no per-basin constant can undo. If these agreed, the + two forms would be indistinguishable and the comparison pointless. + """ + tf = _thermal_forcing([0.5, 2.0]) + mean_tf = _thermal_forcing([1.25, 1.25]) + + local = local_quadratic_melt(1.0, tf, SALINITY, SLOPE, constants=MALI) + semi_local = local_quadratic_melt(1.0, tf, SALINITY, SLOPE, + thermal_forcing_avg=mean_tf, + constants=MALI) + + ratio = (local / semi_local).values + assert not np.allclose(ratio, ratio[0]) + + +def test_year_length_differs_between_the_two_constant_sets(): + """ + MALI uses a 365-day year and the protocol's reference implementation + 365.2422 days. The 0.066% difference is small but showed up as + unexplained noise before it was tracked down, so it is pinned here. + """ + ratio = MULTIMELT.seconds_per_year / MALI.seconds_per_year + assert ratio == pytest.approx(1.000663562, rel=1.0e-8) + + +def test_melt_is_independent_of_ice_density_for_mali_constants(): + """ + The ice density appears once in the formula and once in the conversion + to a mass flux, so for MALI's constants the two cancel exactly. The + protocol's worked example uses 917 and 918, so they do not cancel there + -- which is why both are carried separately. + """ + assert MALI.rho_ice == MALI.rho_ice_flux + assert MULTIMELT.rho_ice != MULTIMELT.rho_ice_flux diff --git a/tests/landice/ismip7_calibration/test_remap_masks.py b/tests/landice/ismip7_calibration/test_remap_masks.py new file mode 100644 index 0000000000..0c05529f77 --- /dev/null +++ b/tests/landice/ismip7_calibration/test_remap_masks.py @@ -0,0 +1,79 @@ +""" +Tests for the ISMIP7 mask checks. + +The basin numbering is the most dangerous thing to get wrong here: ISMIP7 +is 0-based and MALI is 1-based, so confusing them mis-assigns every basin +while still producing plausible-looking aggregates. These test the checks +that catch it. +""" + +import logging + +import numpy as np +import pytest +import xarray as xr + +from compass.landice.tests.ismip7_calibration.ais.remap_masks import ( + EASTERN_AMUNDSEN_BASIN, + REGION_CODES, + _check_shelves_are_in_their_basin, +) + + +def _logger(): + logger = logging.getLogger('test_remap_masks') + logger.addHandler(logging.NullHandler()) + return logger + + +def _masks(basins, regions): + ds = xr.Dataset() + ds['ismip7BasinNumber'] = ('nCells', np.asarray(basins, dtype=np.int32)) + ds['ismip7ShelfRegion'] = ('nCells', np.asarray(regions, dtype=np.int32)) + return ds + + +def test_shelves_in_the_eastern_amundsen_pass(): + """PIG and Dotson both drain into ISMIP7 basin 9.""" + ds = _masks(basins=[EASTERN_AMUNDSEN_BASIN] * 4 + [3], + regions=[REGION_CODES['pig'], REGION_CODES['pig'], + REGION_CODES['dotson'], REGION_CODES['dotson'], + REGION_CODES['none']]) + + _check_shelves_are_in_their_basin(ds, _logger()) + + +def test_off_by_one_basin_numbering_is_caught(): + """ + Using MALI's 1-based numbering where ISMIP7's 0-based is expected puts + the shelves in basin 10 instead of 9. That is the mistake this check + exists for. + """ + ds = _masks(basins=[EASTERN_AMUNDSEN_BASIN + 1] * 4, + regions=[REGION_CODES['pig']] * 2 + + [REGION_CODES['dotson']] * 2) + + with pytest.raises(ValueError, match='numbering is probably off'): + _check_shelves_are_in_their_basin(ds, _logger()) + + +def test_a_few_stray_cells_are_tolerated(): + """ + Nearest-neighbour remapping can put a cell or two on a basin boundary + into a neighbour, which is not an error. + """ + basins = [EASTERN_AMUNDSEN_BASIN] * 19 + [8] + ds = _masks(basins=basins, regions=[REGION_CODES['pig']] * 20) + + _check_shelves_are_in_their_basin(ds, _logger()) + + +def test_a_mesh_with_no_shelf_cells_is_an_error(): + """ + If the shelf mask misses the mesh entirely, term J4 would have nothing + to aggregate over and would fail later and less clearly. + """ + ds = _masks(basins=[1, 2, 3], regions=[REGION_CODES['none']] * 3) + + with pytest.raises(ValueError, match='nothing to aggregate'): + _check_shelves_are_in_their_basin(ds, _logger()) diff --git a/tests/landice/ismip7_calibration/test_terms.py b/tests/landice/ismip7_calibration/test_terms.py new file mode 100644 index 0000000000..13cd285f2e --- /dev/null +++ b/tests/landice/ismip7_calibration/test_terms.py @@ -0,0 +1,188 @@ +""" +Tests for the area-weighted, mesh-agnostic objective-function terms. + +The point of these functions is that they give the *same* answer as the +upstream toolbox on a uniform grid, while being correct on a +variable-resolution mesh where the toolbox's single scalar cell area is not. +Both halves of that claim are tested here. +""" + +import numpy as np +import pytest +import xarray as xr + +from compass.landice.tests.ismip7_calibration.terms import ( + KG_PER_GT, + average_by_group, + integrate_by_group, + stack_cells, + uniform_area, +) + + +def _structured(melt, groups, resolution=8000.0): + """Build a small structured-grid case.""" + ny, nx = melt.shape + coords = {'y': np.arange(ny, dtype=float), + 'x': np.arange(nx, dtype=float)} + melt_da = xr.DataArray(melt, dims=('y', 'x'), coords=coords) + groups_da = xr.DataArray(groups, dims=('y', 'x'), coords=coords) + area = uniform_area(melt_da, resolution, ('y', 'x')) + mask = xr.ones_like(melt_da, dtype=bool) + return melt_da, area, mask, groups_da + + +def test_uniform_area_matches_the_toolbox_convention(): + """ + On a uniform grid the integral must equal the toolbox's + ``sum * reso**2 / 1e12``. + """ + melt = np.array([[1.0, 2.0], [3.0, 4.0]]) + groups = np.array([[0, 0], [1, 1]]) + resolution = 8000.0 + melt_da, area, mask, groups_da = _structured(melt, groups, resolution) + + result = integrate_by_group(melt_da, area, mask, groups_da, ('y', 'x')) + + expected_0 = (1.0 + 2.0) * resolution**2 / KG_PER_GT + expected_1 = (3.0 + 4.0) * resolution**2 / KG_PER_GT + assert float(result.sel(basins=0)) == pytest.approx(expected_0) + assert float(result.sel(basins=1)) == pytest.approx(expected_1) + + +def test_integral_is_area_weighted_on_a_variable_mesh(): + """ + With unequal cell areas the integral must weight by area. A plain sum + would give the wrong answer, which is the bug these functions exist to + avoid on a 4-20 km mesh. + """ + melt = xr.DataArray([1.0, 1.0, 1.0], dims='nCells') + area = xr.DataArray([1.0e9, 2.0e9, 4.0e9], dims='nCells') + groups = xr.DataArray([0, 0, 0], dims='nCells') + mask = xr.ones_like(melt, dtype=bool) + + result = integrate_by_group(melt, area, mask, groups, ('nCells',)) + + assert float(result.sel(basins=0)) == pytest.approx(7.0e9 / KG_PER_GT) + + +def test_mean_is_area_weighted_not_a_plain_mean(): + """A big cold cell must outweigh a small warm one.""" + melt = xr.DataArray([0.0, 100.0], dims='nCells') + area = xr.DataArray([9.0, 1.0], dims='nCells') + groups = xr.DataArray([0, 0], dims='nCells') + mask = xr.ones_like(melt, dtype=bool) + + result = average_by_group(melt, area, mask, groups, ('nCells',)) + + # area-weighted mean is 10, the plain mean would be 50 + assert float(result.sel(basins=0)) == pytest.approx(10.0) + + +def test_mask_excludes_cells_from_both_sums(): + """ + A masked cell must not contribute to the numerator or the denominator + of an area-weighted mean. + """ + melt = xr.DataArray([10.0, 1000.0], dims='nCells') + area = xr.DataArray([1.0, 1.0], dims='nCells') + groups = xr.DataArray([0, 0], dims='nCells') + mask = xr.DataArray([True, False], dims='nCells') + + result = average_by_group(melt, area, mask, groups, ('nCells',)) + + assert float(result.sel(basins=0)) == pytest.approx(10.0) + + +def test_nan_melt_does_not_bias_the_mean_denominator(): + """ + Where melt is NaN the cell's area must be left out of the weights, or + the mean is diluted toward zero. + """ + melt = xr.DataArray([10.0, np.nan], dims='nCells') + area = xr.DataArray([1.0, 99.0], dims='nCells') + groups = xr.DataArray([0, 0], dims='nCells') + mask = xr.ones_like(melt, dtype=bool) + + result = average_by_group(melt, area, mask, groups, ('nCells',)) + + assert float(result.sel(basins=0)) == pytest.approx(10.0) + + +def test_groups_are_kept_separate(): + """Cells in different groups must not be mixed.""" + melt = xr.DataArray([1.0, 2.0, 3.0, 4.0], dims='nCells') + area = xr.DataArray([1.0e12, 1.0e12, 1.0e12, 1.0e12], dims='nCells') + groups = xr.DataArray([0, 1, 0, 1], dims='nCells') + mask = xr.ones_like(melt, dtype=bool) + + result = integrate_by_group(melt, area, mask, groups, ('nCells',)) + + assert float(result.sel(basins=0)) == pytest.approx(4.0) + assert float(result.sel(basins=1)) == pytest.approx(6.0) + + +def test_empty_group_becomes_nan_not_zero(): + """ + A group with no contributing cells must drop out of the objective + rather than be treated as a real zero. + """ + melt = xr.DataArray([1.0, 2.0], dims='nCells') + area = xr.DataArray([1.0e12, 1.0e12], dims='nCells') + groups = xr.DataArray([0, 1], dims='nCells') + mask = xr.DataArray([True, False], dims='nCells') + + result = integrate_by_group(melt, area, mask, groups, ('nCells',)) + + assert np.isnan(float(result.sel(basins=1))) + + +def test_extra_dimensions_are_preserved(): + """ + Terms carry parameter, model and year dimensions through the + aggregation untouched. + """ + melt = xr.DataArray(np.ones((3, 4)), dims=('p1', 'nCells')) + area = xr.DataArray(np.full(4, 1.0e12), dims='nCells') + groups = xr.DataArray([0, 0, 1, 1], dims='nCells') + mask = xr.ones_like(area, dtype=bool) + + result = integrate_by_group(melt, area, mask, groups, ('nCells',)) + + assert result.sizes['p1'] == 3 + assert result.sizes['basins'] == 2 + + +def test_aggregation_is_linear_in_melt(): + """ + Melt is proportional to the calibration parameter, and the calibration + relies on the aggregates being proportional to it too. + """ + melt = xr.DataArray([1.0, 2.0, 3.0], dims='nCells') + area = xr.DataArray([1.0e12, 2.0e12, 3.0e12], dims='nCells') + groups = xr.DataArray([0, 0, 0], dims='nCells') + mask = xr.ones_like(melt, dtype=bool) + + once = integrate_by_group(melt, area, mask, groups, ('nCells',)) + twice = integrate_by_group(2.0 * melt, area, mask, groups, ('nCells',)) + + assert float(twice.sel(basins=0)) == pytest.approx( + 2.0 * float(once.sel(basins=0))) + + +def test_stack_cells_flattens_structured_dimensions(): + """A (y, x) field is flattened onto one cell dimension.""" + da = xr.DataArray(np.arange(6.0).reshape(2, 3), dims=('y', 'x')) + + stacked = stack_cells(da, ('y', 'x')) + + assert stacked.sizes == {'_cell': 6} + assert np.array_equal(np.sort(stacked.values), np.arange(6.0)) + + +def test_stack_cells_rejects_a_missing_dimension(): + """A typo in the cell dimensions should fail loudly.""" + da = xr.DataArray(np.zeros(3), dims='nCells') + + with pytest.raises(ValueError, match='not found in dimensions'): + stack_cells(da, ('ncells',)) diff --git a/tests/landice/ismip7_calibration/test_toolbox.py b/tests/landice/ismip7_calibration/test_toolbox.py new file mode 100644 index 0000000000..7bac2cb274 --- /dev/null +++ b/tests/landice/ismip7_calibration/test_toolbox.py @@ -0,0 +1,47 @@ +""" +Tests for the vendored ISMIP7 parameter-selection toolbox. + +The vendored file must stay byte-for-byte identical to upstream. The +published calibration numbers depend on it, so an accidental edit or a +partial update has to fail loudly rather than quietly change results. +""" + +import pytest + +from compass.landice.tests.ismip7_calibration import toolbox + + +def test_vendored_toolbox_matches_the_recorded_checksum(): + """The copy on disk must be the upstream revision PROVENANCE.md names.""" + assert toolbox.file_sha256() == toolbox.EXPECTED_SHA256 + + +def test_check_integrity_passes_for_the_shipped_copy(): + """The import-time check must not raise for an unmodified checkout.""" + toolbox.check_integrity() + + +def test_check_integrity_reports_a_mismatch(monkeypatch): + """A modified file must raise, naming both checksums.""" + monkeypatch.setattr(toolbox, 'EXPECTED_SHA256', '0' * 64) + + with pytest.raises(RuntimeError, match='does not match the recorded'): + toolbox.check_integrity() + + +def test_provenance_is_recorded(): + """The upstream commit and URL must be pinned in the module.""" + assert len(toolbox.UPSTREAM_COMMIT) == 40 + assert toolbox.UPSTREAM_URL.startswith('https://github.com/ismip/') + + +def test_the_functions_the_calibration_calls_are_present(): + """ + Guard against an upstream refactor silently removing something the + calibration drives. + """ + module = toolbox.parameter_selection_toolbox + for name in ('calculate_objective_function', 'optimise_deltaT', + 'select_optimal_deltaT', + 'select_subensemble_using_optimal_deltaT'): + assert callable(getattr(module, name))