Merge pull request #386 from matthewhoffman/landice/mismipplus_calvin…

…gdt_test

MALI calving dt convergence tests

This merge adds a new set of tests for MALI that explore convergence with timestep of the calving physics. This is added in a new test group that runs a configuration multiple times with the adaptive timestepper restricted by different fractions of the calculated calving CFL limit. The results are then analyzed and a summary plot is produced. This test can be applied to a number of different calving laws and meshes and using either data velocity or prognostic velocity. These set of tests are not intended to be run routinely as part of integration, but instead are a tool for evaluating development of calving physics. A test suite is also included that runs a number of these tests that have been evaluated as being useful. There are some known ways this test group could be improved, but given its limited intended usage, they are not worth developing further at this time.

A MISMIP+ smoke test is also added as part of this PR. The calving tests were first developed to use the MISMIP+ configuration, which was not yet in COMPASS. Later, things were reorganized to make the calving tests more general and support multiple meshes. Having a MISMIP+ smoke test is generally useful, so that test was retained in this branch after reorganizing things.

Note that this PR includes an update to the MALI-Dev submodule to include the merged MALI PR here: MALI-Dev/E3SM#33.

compass/landice/__init__.py

@@ -1,5 +1,6 @@
 from compass.mpas_core import MpasCore
 from compass.landice.tests.antarctica import Antarctica
+from compass.landice.tests.calving_dt_convergence import CalvingDtConvergence
 from compass.landice.tests.circular_shelf import CircularShelf
 from compass.landice.tests.dome import Dome
 from compass.landice.tests.eismint2 import Eismint2
@@ -9,6 +10,7 @@
 from compass.landice.tests.hydro_radial import HydroRadial
 from compass.landice.tests.kangerlussuaq import Kangerlussuaq
 from compass.landice.tests.koge_bugt_s import KogeBugtS
+from compass.landice.tests.mismipplus import MISMIPplus
 from compass.landice.tests.thwaites import Thwaites
 
 
@@ -24,6 +26,7 @@ def __init__(self):
         super().__init__(name='landice')
 
         self.add_test_group(Antarctica(mpas_core=self))
+        self.add_test_group(CalvingDtConvergence(mpas_core=self))
         self.add_test_group(CircularShelf(mpas_core=self))
         self.add_test_group(Dome(mpas_core=self))
         self.add_test_group(Eismint2(mpas_core=self))
@@ -33,4 +36,5 @@ def __init__(self):
         self.add_test_group(HydroRadial(mpas_core=self))
         self.add_test_group(Kangerlussuaq(mpas_core=self))
         self.add_test_group(KogeBugtS(mpas_core=self))
+        self.add_test_group(MISMIPplus(mpas_core=self))
         self.add_test_group(Thwaites(mpas_core=self))

compass/landice/suites/calving_dt_convergence.txt

@@ -0,0 +1,9 @@
+landice/calving_dt_convergence/mismip+.specified_calving_velocity.none
+landice/calving_dt_convergence/mismip+.von_Mises_stress.none
+landice/calving_dt_convergence/humboldt.specified_calving_velocity.none
+landice/calving_dt_convergence/humboldt.von_Mises_stress.none
+landice/calving_dt_convergence/thwaites.specified_calving_velocity.none
+landice/calving_dt_convergence/thwaites.von_Mises_stress.none
+landice/calving_dt_convergence/mismip+.von_Mises_stress.FO
+landice/calving_dt_convergence/humboldt.von_Mises_stress.FO
+landice/calving_dt_convergence/thwaites.von_Mises_stress.FO

compass/landice/tests/calving_dt_convergence/__init__.py

@@ -0,0 +1,40 @@
+from compass.testgroup import TestGroup
+from compass.landice.tests.calving_dt_convergence.dt_convergence_test \
+     import DtConvergenceTest
+
+
+class CalvingDtConvergence(TestGroup):
+    """
+    A test group for a series of test cases for assessing timestep
+    convergence related to calving.
+    This test group uses a pre-made mesh files and various calving
+    and velocity options.
+    """
+    def __init__(self, mpas_core):
+        """
+        mpas_core : compass.landice.Landice
+            the MPAS core that this test group belongs to
+        """
+        super().__init__(mpas_core=mpas_core, name='calving_dt_convergence')
+
+        for mesh in [
+                     'mismip+',
+                     'humboldt',
+                     'thwaites'
+                    ]:
+            for calving in [
+                            'specified_calving_velocity',
+                            'von_Mises_stress',
+                            'eigencalving'
+                           ]:
+                for velo in [
+                             'none',
+                             'FO'
+                            ]:
+                    if (calving == 'specified_calving_velocity' and
+                       velo == 'FO'):
+                        continue  # This combination is not useful
+                    self.add_test_case(DtConvergenceTest(test_group=self,
+                                                         mesh=mesh,
+                                                         calving=calving,
+                                                         velo=velo))

compass/landice/tests/calving_dt_convergence/dt_convergence_test/__init__.py

@@ -0,0 +1,115 @@
+from compass.testcase import TestCase
+from compass.landice.tests.calving_dt_convergence.run_model import RunModel
+# from compass.validate import compare_variables  # not currently used
+import numpy
+import netCDF4
+import matplotlib.pyplot as plt
+import matplotlib.cm
+
+
+class DtConvergenceTest(TestCase):
+    """
+    A test case for running the same configuration with a series of
+    values for config_adaptive_timestep_calvingCFL_fraction
+    to check for convergence
+
+    Attributes
+    ----------
+    name : str
+        The name of the case
+    """
+
+    def __init__(self, test_group, mesh, calving, velo):
+        """
+        Create the test case
+
+        Parameters
+        ----------
+        test_group : compass.landice.tests.calving_dt_convergence.CalvingDtConvergence
+            The test group that this test case belongs to
+
+        """
+        self.name = f'calving_dt_convergence_test_{mesh}_{calving}_{velo}'
+        subdir = f'{mesh}.{calving}.{velo}'
+        super().__init__(test_group=test_group, name=self.name, subdir=subdir)
+
+        cores = 36
+        min_cores = 4
+
+        # Do fewer runs if FO solver
+        if velo == 'FO':
+            self.fractions = numpy.arange(0.25, 2.3, 0.25)
+        else:
+            self.fractions = numpy.arange(0.2, 3.1, 0.2)
+
+        for frac in self.fractions:
+            name = f'frac{frac:.2f}'
+            step = RunModel(test_case=self, name=name,
+                            mesh=mesh,
+                            calving=calving,
+                            velo=velo,
+                            calv_dt_frac=frac,
+                            cores=cores, min_cores=min_cores, threads=1)
+            self.add_step(step)
+
+    # no configure() method is needed
+
+    # no run() method is needed
+
+    def validate(self):
+        """
+        Test cases can override this method to perform validation of variables
+        and timers
+        """
+        # If variable comparison is added, need to uncomment line 3 as well
+        # variables = ['thickness', 'surfaceSpeed']
+        # compare_variables(test_case=self, variables=variables,
+        #                   filename1='full_run/output.nc',
+        #                   filename2='restart_run/output.nc')
+
+        # plot results
+        fig, ax = plt.subplots(4, figsize=(10, 7))
+        ax[0].set(xlabel='year', ylabel='calving flux (kg/yr)')
+        ax[1].set(xlabel='year', ylabel='cum. calving flux (kg)')
+        ax[2].set(xlabel='year', ylabel='actual dt to calving dt ratio')
+        ax[3].set(xlabel='fraction')
+        ax[3].set_ylabel('# warnings', color='c')
+        ax2 = ax[3].twinx()
+        ax2.set_ylabel('fraction with warnings', color='g')
+        colors = matplotlib.cm.jet(numpy.linspace(0, 1, len(self.fractions)))
+        nWarn = numpy.zeros([len(self.fractions)])
+        nTimesteps = numpy.zeros([len(self.fractions)])
+
+        i = 0
+        for frac in self.fractions:
+            name = f'frac{frac:.2f}'
+            f = netCDF4.Dataset(f'{name}/globalStats.nc', 'r')
+            yr = f.variables['daysSinceStart'][:] / 365.0
+            calv = f.variables['totalCalvingFlux'][:]
+            deltat = f.variables['deltat'][:]
+            ax[0].plot(yr[1:], calv[1:], '-', label=f'{frac:.2f}',
+                       color=colors[i])
+
+            ax[1].plot(yr[1:], (calv[1:]*deltat[1:]).cumsum(), '-',
+                       color=colors[i])
+
+            ratio = f.variables['dtCalvingCFLratio'][:]
+            ax[2].plot(yr[1:], numpy.ones(yr[1:].shape) * frac, 'k:',
+                       label=f'{frac:.2f}')
+            ax[2].plot(yr[1:], ratio[1:], '-', label=f'{frac:.2f}',
+                       color=colors[i])
+
+            # Now count errors
+            file = open(f"{name}/log.landice.0000.out", "r")
+            logcontents = file.read()
+            # get number of occurrences of the substring in the string
+            nWarn[i] = logcontents.count("WARNING: Failed to ablate")
+            nTimesteps[i] = logcontents.count("Starting timestep number")
+            ax[3].plot(frac, nWarn[i], 'co')
+            ax2.plot(frac, nWarn[i]/nTimesteps[i], 'gx')
+
+            f.close()
+            i += 1
+
+        ax[0].legend(loc='best', prop={'size': 5})
+        plt.savefig(f'{self.name}_comparison_plot.png', dpi=150)

compass/landice/tests/calving_dt_convergence/namelist.landice

@@ -0,0 +1,29 @@
+config_start_time = '0000-01-01_00:00:00'
+config_run_duration = '0025-00-00_00:00:00'
+config_adaptive_timestep = .true.
+config_adaptive_timestep_CFL_fraction = 0.95
+config_adaptive_timestep_include_calving = .true.
+config_adaptive_timestep_calvingCFL_fraction = 1.0
+config_block_decomp_file_prefix = 'graph.info.part.'
+config_velocity_solver = 'FO'
+config_thickness_advection = 'fo'
+config_tracer_advection = 'none'
+config_thermal_solver = 'none'
+config_temperature_init = 'file'
+
+config_calving = 'von_Mises_stress'
+config_floating_von_Mises_threshold_stress_source = 'scalar'
+config_floating_von_Mises_threshold_stress = 20000.0
+config_grounded_von_Mises_threshold_stress_source = 'scalar'
+config_grounded_von_Mises_threshold_stress = 20000.0
+config_calving_specified_source = 'const'
+config_calving_velocity_const = 0.000127388
+config_calving_eigencalving_parameter_scalar_value = 1.0e17
+config_distribute_unablatedVolumeDynCell = .false.
+config_calving_error_threshold = 1000.0
+
+config_write_output_on_startup = .true.
+config_AM_globalStats_enable = .true.
+config_AM_globalStats_compute_interval = '0000-00-00_01:00:00'
+config_AM_globalStats_compute_on_startup = .true.
+config_AM_globalStats_write_on_startup = .true.

compass/landice/tests/calving_dt_convergence/run_model.py

@@ -0,0 +1,149 @@
+from compass.model import make_graph_file, run_model
+from compass.step import Step
+
+
+class RunModel(Step):
+    """
+    A step for performing forward MALI runs as part of calving dt convergence
+    test cases.
+
+    Attributes
+    ----------
+    mesh_file : str
+        name of mesh file used
+
+    suffixes : list of str, optional
+        a list of suffixes for namelist and streams files produced
+        for this step.  Most steps most runs will just have a
+        ``namelist.landice`` and a ``streams.landice`` (the default) but
+        the ``restart_run`` step of the ``restart_test`` runs the model
+        twice, the second time with ``namelist.landice.rst`` and
+        ``streams.landice.rst``
+    """
+    def __init__(self, test_case, name, mesh, calving, velo, calv_dt_frac,
+                 subdir=None, cores=1, min_cores=None, threads=1,
+                 suffixes=None):
+        """
+        Create a new test case
+
+        Parameters
+        ----------
+        test_case : compass.TestCase
+            The test case this step belongs to
+
+        name : str
+            the name of the test case
+
+        mesh : str
+            the name of the mesh to be used. Valid values are:
+            'mismip+', 'humboldt', 'thwaites'
+
+        calving : str
+            the name of the calving option to be used. Valid values are:
+            'specified_calving_velocity', 'eigencalving', 'von_Mises_stress'
+
+        velo : str
+            the name of the velocity setting to use.  Valid values are:
+            'none', 'FO'
+
+        calv_dt_frac : float
+            the value to use for calving dt fraction
+
+        subdir : str, optional
+            the subdirectory for the step.  The default is ``name``
+
+        cores : int, optional
+            the number of cores the step would ideally use.  If fewer cores
+            are available on the system, the step will run on all available
+            cores as long as this is not below ``min_cores``
+
+        min_cores : int, optional
+            the number of cores the step requires.  If the system has fewer
+            than this number of cores, the step will fail
+
+        threads : int, optional
+            the number of threads the step will use
+
+        suffixes : list of str, optional
+            a list of suffixes for namelist and streams files produced
+            for this step.  Most run steps will just have a
+            ``namelist.landice`` and a ``streams.landice`` (the default) but
+            the ``restart_run`` step of the ``restart_test`` runs the model
+            twice, the second time with ``namelist.landice.rst`` and
+            ``streams.landice.rst``
+        """
+
+        if suffixes is None:
+            suffixes = ['landice']
+        self.suffixes = suffixes
+        if min_cores is None:
+            min_cores = cores
+        super().__init__(test_case=test_case, name=name, subdir=subdir,
+                         cores=cores, min_cores=min_cores, threads=threads)
+
+        # download and link the mesh
+        if mesh == 'mismip+':
+            self.mesh_file = 'MISMIP_2km_20220502.nc'
+            self.add_input_file(filename='albany_input.yaml',
+                                package='compass.landice.tests.mismipplus',
+                                copy=True)
+        elif mesh == 'humboldt':
+            self.mesh_file = 'Humboldt_3to30km_r04_20210615.nc'
+            self.add_input_file(filename='albany_input.yaml',
+                                package='compass.landice.tests.humboldt',
+                                copy=True)
+        elif mesh == 'thwaites':
+            self.mesh_file = 'thwaites.4km.210608.nc'
+            self.add_input_file(filename='albany_input.yaml',
+                                package='compass.landice.tests.thwaites',
+                                copy=True)
+        self.add_input_file(filename=self.mesh_file, target=self.mesh_file,
+                            database='')
+
+        for suffix in suffixes:
+            self.add_namelist_file(
+                'compass.landice.tests.calving_dt_convergence',
+                'namelist.landice',
+                out_name='namelist.{}'.format(suffix))
+
+            # Modify nl & streams for the requested options
+            options = {'config_calving': f"'{calving}'",
+                       'config_velocity_solver': f"'{velo}'",
+                       'config_adaptive_timestep_calvingCFL_fraction':
+                       f"{calv_dt_frac}"}
+            if velo == 'FO':
+                # Make FO runs shorter
+                options['config_run_duration'] = "'0005-00-00_00:00:00'"
+            if mesh == 'thwaites':
+                # adjust calving params for Thwaites
+                options['config_floating_von_Mises_threshold_stress'] = "50000"
+                options['config_grounded_von_Mises_threshold_stress'] = "50000"
+                options['config_calving_velocity_const'] = "0.000477705"
+                # (15km/yr)
+
+            self.add_namelist_options(options=options,
+                                      out_name='namelist.{}'.format(suffix))
+
+            stream_replacements = {'INPUT_FILE': self.mesh_file}
+            self.add_streams_file(
+                'compass.landice.tests.calving_dt_convergence',
+                'streams.landice.template',
+                out_name='streams.{}'.format(suffix),
+                template_replacements=stream_replacements)
+
+        self.add_model_as_input()
+
+        self.add_output_file(filename='output.nc')
+        self.add_output_file(filename='globalStats.nc')
+
+    # no setup() is needed
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        make_graph_file(mesh_filename=self.mesh_file,
+                        graph_filename='graph.info')
+        for suffix in self.suffixes:
+            run_model(step=self, namelist='namelist.{}'.format(suffix),
+                      streams='streams.{}'.format(suffix))