MALI calving dt convergence tests

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))