Merge pull request #737 from andrewdnolan/mismip+_spinup

Port the MISMIP+ spinup testcase from the old version of compass 

Currently the only test case within the MISMIP+ test group is a SmokeTest, which uses an archived "spun up" initial condition. This PR adds a SpinUp test case that generates the initial conditions and other input files (e.g. namelist and streams) needed to run MALI for a ~20,000 year spin up. The SpinUp testcase supports an arbitrary resolution, which needs to be specified in the config file at the time of running compass setup. Furthermore, the user can specify a gutterLength in the config file to extend the eastern boundary of the MISMIP+ domain, which was added to support simulations using a dynamic claving law that will results in an irregular grounding line. To this end, we have added a SetupMesh step, used by the SpinUp testcase that creates a MALI mesh (and initial conditions) following the MISMIP+ requirements (Asay-Davis et al. 2016).

We've also slightly altered the SmokeTest to explicitly list the supported resolutions (e.g. compass list | grep mismipplus). Currently this is only 2000m, but in the future the SpinUp testcase will be used to generated the spun up conditions needed to support additional resolutions of the SmokeTest.

compass/landice/tests/mismipplus/__init__.py

@@ -1,11 +1,11 @@
-from compass.testgroup import TestGroup
 from compass.landice.tests.mismipplus.smoke_test import SmokeTest
+from compass.landice.tests.mismipplus.spin_up import SpinUp
+from compass.testgroup import TestGroup
 
 
 class MISMIPplus(TestGroup):
     """
     A test group for MISMIP+ test cases.
-    This test group uses a pre-made mesh file.
     """
     def __init__(self, mpas_core):
         """
@@ -14,4 +14,8 @@ def __init__(self, mpas_core):
         """
         super().__init__(mpas_core=mpas_core, name='mismipplus')
 
-        self.add_test_case(SmokeTest(test_group=self))
+        self.add_test_case(SpinUp(test_group=self))
+
+        for resolution in [2000]:
+            self.add_test_case(SmokeTest(test_group=self,
+                                         resolution=resolution))

compass/landice/tests/mismipplus/albany_input.yaml

@@ -3,6 +3,24 @@
 ANONYMOUS:
   Build Type: Tpetra
 
+  Problem:
+    Dirichlet BCs:
+      SDBC on NS dirichlet for DOF U1 prescribe Field: dirichlet_field
+    LandIce BCs:
+      BC 0:
+        Side Set Name: basalside
+        Type: Basal Friction
+        Basal Friction Coefficient:
+          Type: Power Law
+          Mu Type: Field
+          Power Exponent: 0.3333333333
+          Effective Pressure Type: Field
+          Zero Beta On Floating Ice: true
+    LandIce Viscosity:
+      Type: 'Glen''s Law'
+      Flow Rate Type: Uniform
+      Glen's Law A: 6.338e-25
+
 # Discretization Description
 #  Discretization:
 #    Exodus Output File Name: albany_output.exo

compass/landice/tests/mismipplus/mismipplus.cfg

@@ -0,0 +1,8 @@
+[mismipplus]
+
+# the number of cells per core to aim for
+goal_cells_per_core = 300
+
+# the approximate maximum number of cells per core (the test will fail if too
+# few cores are available)
+max_cells_per_core = 5000

compass/landice/tests/mismipplus/namelist.landice

@@ -1,24 +1,19 @@
-config_start_time = '0000-01-01_00:00:00'
-config_run_duration = '0005-00-00_00:00:00'
+config_start_time = '0001-01-01_00:00:00'
+config_stop_time = '0006-01-01_00:00:00'
 config_adaptive_timestep = .true.
 config_adaptive_timestep_include_calving = .false.
-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_dynamic_thickness = 0.001
+config_dynamic_thickness = 10.0
+config_adaptive_timestep_CFL_fraction = 0.8
 
 config_calving = 'none'
-config_floating_von_Mises_threshold_stress_source = 'scalar'
-config_floating_von_Mises_threshold_stress = 20000.0
-config_calving_specified_source = 'const'
-config_calving_velocity_const = 0.000127388
-config_distribute_unablatedVolumeDynCell = .false.
-config_calving_error_threshold = 1000.0
 
+config_year_digits = 5
 config_write_output_on_startup = .true.
 config_AM_globalStats_enable = .true.
 config_AM_globalStats_compute_interval = '0000-00-00_01:00:00'

compass/landice/tests/mismipplus/run_model.py

@@ -1,3 +1,10 @@
+import os
+
+from compass.landice.tests.mismipplus.tasks import (
+    approx_cell_count,
+    exact_cell_count,
+    get_ntasks_from_cell_count,
+)
 from compass.model import make_graph_file, run_model
 from compass.step import Step
 
@@ -15,8 +22,11 @@ class RunModel(Step):
         the ``restart_run`` step of the ``restart_test`` runs the model
         twice, the second time with ``namelist.landice.rst`` and
         ``streams.landice.rst``
+
+    resolution : float
+        The nominal distance [m] between horizontal grid points (dcEdge).
     """
-    def __init__(self, test_case, name='run_model', subdir=None,
+    def __init__(self, test_case, name, subdir=None, resolution=None,
                  ntasks=1, min_tasks=None, openmp_threads=1, suffixes=None):
         """
         Create a new test case
@@ -26,12 +36,18 @@ def __init__(self, test_case, name='run_model', subdir=None,
         test_case : compass.TestCase
             The test case this step belongs to
 
-        name : str, optional
+        name : str
             the name of the test case
 
         subdir : str, optional
             the subdirectory for the step.  The default is ``name``
 
+        resolution : float, optional
+            The nominal distance [m] between horizontal grid points (dcEdge).
+            For the ``SpinUp`` testcase no value should be passed. Instead the
+            value will be parsed from the configuration options at the
+            time of ``compass setup``.
+
         ntasks : int, optional
             the number of tasks the step would ideally use.  If fewer tasks
             are available on the system, the step will run on all available
@@ -56,17 +72,29 @@ def __init__(self, test_case, name='run_model', subdir=None,
         if suffixes is None:
             suffixes = ['landice']
         self.suffixes = suffixes
+
+        # The condition below will only be true for the `SpinUp` testcase
+        # where resolution is not know at the time of object construction
+        # (e.g. `compass list`), but instead is parsed from the config file
+        # when `compass setup` is executed. For the `SpinUp` testcase the
+        # resolution will be parsed by the testcases `configure`` method,
+        # where the attribute value is also overwritten to properly reflect
+        # the requested value. A default value is initialized here to
+        # follow best python practices and have all the attributes be set
+        # in the constructor.
+        if resolution is None:
+            resolution = 0.0
+        self.resolution = resolution
+
         if min_tasks is None:
             min_tasks = ntasks
-        super().__init__(test_case=test_case, name=name, subdir=subdir,
+
+        super().__init__(test_case=test_case,
+                         name=name,
+                         subdir=subdir,
                          ntasks=ntasks, min_tasks=min_tasks,
                          openmp_threads=openmp_threads)
 
-        # download and link the mesh
-        self.mesh_file = 'MISMIP_2km_20220502.nc'
-        self.add_input_file(filename=self.mesh_file, target=self.mesh_file,
-                            database='')
-
         for suffix in suffixes:
             self.add_namelist_file(
                 'compass.landice.tests.mismipplus', 'namelist.landice',
@@ -85,14 +113,74 @@ def __init__(self, test_case, name='run_model', subdir=None,
         self.add_output_file(filename='output.nc')
         self.add_output_file(filename='globalStats.nc')
 
-    # no setup() is needed
+    def setup(self):
+        """
+        Set the number of MPI tasks based on a tentative scaling of
+        a the `ncells_at_1km_res` heuristic from the config file
+        based on the desired resolution.
+        """
+
+        config = self.config
+        # use the resolution at the time of ``compass setup``
+        resolution = self.resolution
+
+        # default gutter_length is 0. This encompasses `SmokeTest`, which does
+        # not use the gutter_length option.
+        gutter_length = 0.0
+
+        if self.test_case.name == "spin_up":
+            gutter_length = config.getfloat('mesh', 'gutter_length')
+
+        # approximate the number of cells in the mesh, based on the requested
+        # resolution and gutter_length.
+        cell_count = approx_cell_count(resolution, gutter_length)
+
+        # find optimal and minimum number of task for the desired resolution
+        ntasks, min_tasks = get_ntasks_from_cell_count(config, cell_count)
+
+        # set values as attributes
+        self.ntasks, self.min_tasks = (ntasks, min_tasks)
+
+        super().setup()
+
+    def constrain_resources(self, available_resources):
+        """
+        Update resources at runtime from config options
+        """
+
+        config = self.config
+
+        mesh_filename = os.path.join(self.work_dir, 'landice_grid.nc')
+
+        # find the number of cells from the mesh created as part of the
+        # `setup_mesh` step.
+        cell_count = exact_cell_count(mesh_filename)
+
+        # find optimal and minimum number of task for the cell count
+        ntasks, min_tasks = get_ntasks_from_cell_count(config, cell_count)
+
+        # set values as attributes
+        self.ntasks, self.min_tasks = (ntasks, min_tasks)
+
+        super().constrain_resources(available_resources)
 
     def run(self):
         """
         Run this step of the test case
         """
+
+        config = self.config
+
+        if self.test_case.name == "spin_up":
+            # read the density value from config file and update
+            # namelist. This will only work for the `SpinUp` test case.
+            ice_density = config['mesh'].getfloat('ice_density')
+            self.update_namelist_at_runtime(
+                {'config_ice_density': f'{ice_density}'})
+
         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))