Add debris friction test with MISMIP+

compass/landice/suites/full_integration.txt

@@ -18,6 +18,8 @@ landice/dome/2000m/fo_restart_test
 landice/dome/variable_resolution/fo_decomposition_test
 landice/dome/variable_resolution/fo_restart_test
 landice/circular_shelf/decomposition_test
+landice/mismipplus/smoke_test/2000m
+landice/mismipplus/smoke_test/debris_friction/2000m
 landice/greenland/fo_decomposition_test
 landice/greenland/fo_restart_test
 landice/thwaites/fo_decomposition_test

compass/landice/tests/mismipplus/__init__.py

@@ -19,3 +19,6 @@ def __init__(self, mpas_core):
         for resolution in [2000]:
             self.add_test_case(SmokeTest(test_group=self,
                                          resolution=resolution))
+            self.add_test_case(SmokeTest(test_group=self,
+                                         resolution=resolution,
+                                         debris_friction=True))

compass/landice/tests/mismipplus/albany_input_debrisfriction.yaml

@@ -0,0 +1,241 @@
+%YAML 1.1
+---
+ANONYMOUS:
+  Problem:
+    LandIce Field Norm:
+      sliding_velocity_basalside:
+        Regularization Type: Given Value
+        Regularization Value: 1.0e-4
+    LandIce BCs:
+      BC 0:
+        Cubature Degree: 4
+        Basal Friction Coefficient:
+          Type: Debris Friction
+          Power Exponent: 0.3333333333
+          Effective Pressure Type: Hydrostatic Computed at Nodes
+          Effective Pressure Regularization: 1.0e3
+          Sliding Velocity Regularization: 5.0e2
+          Flow Rate Type: Constant
+          Flow Rate: 4.537e-24      # Pa^-n s^-1
+  # # Discretization Description
+  # Discretization:
+  #   Exodus Output File Name: albany_output.exo
+  #   Side Set Discretizations:
+  #     basalside:
+  #       Exodus Output File Name: albany_output_basal.exo
+
+  Piro: 
+#   Nonlinear Solver Information
+    NOX: 
+      Nonlinear Solver: Line Search Based
+      Line Search: 
+        Full Step: 
+          Full Step: 1.0e+00
+        Method: Backtrack
+      Solver Options: 
+        Status Test Check Type: Minimal
+      Status Tests: 
+        Test Type: Combo
+        Combo Type: OR
+        Number of Tests: 2
+        Test 0: 
+          Test Type: NormF
+          Norm Type: Two Norm
+          Scale Type: Scaled
+          Tolerance: 1.0e-05
+        Test 1: 
+          Test Type: MaxIters
+          Maximum Iterations: 50
+      Printing: 
+        Output Precision: 3
+        Output Processor: 0
+        Output Information: 
+          Error: true
+          Warning: true
+          Outer Iteration: true
+          Parameters: false
+          Details: false
+          Linear Solver Details: false
+          Stepper Iteration: true
+          Stepper Details: true
+          Stepper Parameters: true
+
+      Direction: 
+        Method: Newton
+        Newton: 
+          Forcing Term Method: Type 2
+          Rescue Bad Newton Solve: true
+          Linear Solver:
+            Write Linear System: false
+            Tolerance: 1.0e-8
+
+          Stratimikos Linear Solver: 
+            Stratimikos:
+
+#             Linear Solver Information
+              Linear Solver Type: Belos
+              Linear Solver Types: 
+                Belos:
+                  Solver Type: Block GMRES
+                  Solver Types:
+                    Block GMRES:
+                      Output Frequency: 20
+                      Output Style: 1
+                      Verbosity: 33
+                      Maximum Iterations: 200
+                      Block Size: 1
+                      Num Blocks: 200
+                      Flexible Gmres: false
+                  VerboseObject:
+                    Output File: none
+                    Verbosity Level: low
+
+#             Preconditioner Information
+              Preconditioner Type: MueLu
+              Preconditioner Types: 
+
+                Ifpack2: 
+                  Overlap: 1
+                  Prec Type: ILUT                
+
+                MueLu: 
+                  Matrix: 
+                    PDE equations: 2
+                  Factories: 
+                    myLineDetectionFact: 
+                      factory: LineDetectionFactory
+                      'linedetection: orientation': coordinates
+                    mySemiCoarsenPFact1: 
+                      factory: SemiCoarsenPFactory
+                      'semicoarsen: coarsen rate': 14
+                    UncoupledAggregationFact2: 
+                      factory: UncoupledAggregationFactory
+                      'aggregation: ordering': graph
+                      'aggregation: max selected neighbors': 0
+                      'aggregation: min agg size': 3
+                      'aggregation: phase3 avoid singletons': true
+                    MyCoarseMap2: 
+                      factory: CoarseMapFactory
+                      Aggregates: UncoupledAggregationFact2
+                    myTentativePFact2: 
+                      'tentative: calculate qr': true
+                      factory: TentativePFactory
+                      Aggregates: UncoupledAggregationFact2
+                      CoarseMap: MyCoarseMap2
+                    mySaPFact2: 
+                      'sa: eigenvalue estimate num iterations': 10
+                      'sa: damping factor': 1.33333e+00
+                      factory: SaPFactory
+                      P: myTentativePFact2
+                    myTransferCoordinatesFact: 
+                      factory: CoordinatesTransferFactory
+                      CoarseMap: MyCoarseMap2
+                      Aggregates: UncoupledAggregationFact2
+                    myTogglePFact: 
+                      factory: TogglePFactory
+                      'semicoarsen: number of levels': 2
+                      TransferFactories: 
+                        P1: mySemiCoarsenPFact1
+                        P2: mySaPFact2
+                        Ptent1: mySemiCoarsenPFact1
+                        Ptent2: myTentativePFact2
+                        Nullspace1: mySemiCoarsenPFact1
+                        Nullspace2: myTentativePFact2
+                    myRestrictorFact: 
+                      factory: TransPFactory
+                      P: myTogglePFact
+                    myToggleTransferCoordinatesFact: 
+                      factory: ToggleCoordinatesTransferFactory
+                      Chosen P: myTogglePFact
+                      TransferFactories: 
+                        Coordinates1: mySemiCoarsenPFact1
+                        Coordinates2: myTransferCoordinatesFact
+                    myRAPFact: 
+                      factory: RAPFactory
+                      P: myTogglePFact
+                      R: myRestrictorFact
+                      TransferFactories: 
+                        For Coordinates: myToggleTransferCoordinatesFact
+                    myRepartitionHeuristicFact: 
+                      factory: RepartitionHeuristicFactory
+                      A: myRAPFact
+                      'repartition: min rows per proc': 3000
+                      'repartition: max imbalance': 1.327e+00
+                      'repartition: start level': 1
+                    myZoltanInterface: 
+                      factory: ZoltanInterface
+                      A: myRAPFact
+                      Coordinates: myToggleTransferCoordinatesFact
+                      number of partitions: myRepartitionHeuristicFact
+                    myRepartitionFact: 
+                      factory: RepartitionFactory
+                      A: myRAPFact
+                      Partition: myZoltanInterface
+                      'repartition: remap parts': true
+                      number of partitions: myRepartitionHeuristicFact
+                    myRebalanceProlongatorFact: 
+                      factory: RebalanceTransferFactory
+                      type: Interpolation
+                      P: myTogglePFact
+                      Coordinates: myToggleTransferCoordinatesFact
+                      Nullspace: myTogglePFact
+                    myRebalanceRestrictionFact: 
+                      factory: RebalanceTransferFactory
+                      type: Restriction
+                      R: myRestrictorFact
+                    myRebalanceAFact: 
+                      factory: RebalanceAcFactory
+                      A: myRAPFact
+                      TransferFactories: { }
+                    mySmoother1: 
+                      factory: TrilinosSmoother
+                      type: LINESMOOTHING_BANDEDRELAXATION
+                      'smoother: pre or post': both
+                      ParameterList: 
+                        'relaxation: type': Gauss-Seidel
+                        'relaxation: sweeps': 1
+                        'relaxation: damping factor': 1.0
+                    mySmoother3: 
+                      factory: TrilinosSmoother
+                      type: RELAXATION
+                      'smoother: pre or post': both
+                      ParameterList: 
+                        'relaxation: type': Gauss-Seidel
+                        'relaxation: sweeps': 1
+                        'relaxation: damping factor': 1.0
+                    mySmoother4: 
+                      factory: TrilinosSmoother
+                      type: RELAXATION
+                      'smoother: pre or post': pre
+                      ParameterList: 
+                        'relaxation: type': Gauss-Seidel
+                        'relaxation: sweeps': 4
+                        'relaxation: damping factor': 1.0
+                  Hierarchy: 
+                    max levels: 7
+                    'coarse: max size': 2000
+                    verbosity: None
+                    Finest: 
+                      Smoother: mySmoother1
+                      CoarseSolver: mySmoother4
+                      P: myRebalanceProlongatorFact
+                      Nullspace: myRebalanceProlongatorFact
+                      CoarseNumZLayers: myLineDetectionFact
+                      LineDetection_Layers: myLineDetectionFact
+                      LineDetection_VertLineIds: myLineDetectionFact
+                      A: myRebalanceAFact
+                      Coordinates: myRebalanceProlongatorFact
+                      Importer: myRepartitionFact
+                    All: 
+                      startLevel: 1
+                      Smoother: mySmoother4
+                      CoarseSolver: mySmoother4
+                      P: myRebalanceProlongatorFact
+                      Nullspace: myRebalanceProlongatorFact
+                      CoarseNumZLayers: myLineDetectionFact
+                      LineDetection_Layers: myLineDetectionFact
+                      LineDetection_VertLineIds: myLineDetectionFact
+                      A: myRebalanceAFact
+                      Coordinates: myRebalanceProlongatorFact
+                      Importer: myRepartitionFact
+

compass/landice/tests/mismipplus/namelist.landice

@@ -1,5 +1,5 @@
 config_start_time = '0001-01-01_00:00:00'
-config_stop_time = '0006-01-01_00:00:00'
+config_stop_time = '0003-01-01_00:00:00'
 config_adaptive_timestep = .true.
 config_adaptive_timestep_include_calving = .false.
 config_block_decomp_file_prefix = 'graph.info.part.'

compass/landice/tests/mismipplus/run_model.py

@@ -1,5 +1,7 @@
 import os
 
+import xarray as xr
+
 from compass.landice.tests.mismipplus.tasks import (
     approx_cell_count,
     exact_cell_count,
@@ -27,7 +29,9 @@ class RunModel(Step):
         The nominal distance [m] between horizontal grid points (dcEdge).
     """
     def __init__(self, test_case, name, subdir=None, resolution=None,
-                 ntasks=1, min_tasks=None, openmp_threads=1, suffixes=None):
+                 ntasks=1, min_tasks=None, openmp_threads=1, suffixes=None,
+                 albany_input_yaml='albany_input.yaml',
+                 debris_friction=False):
         """
         Create a new test case
 
@@ -67,11 +71,20 @@ def __init__(self, test_case, name, subdir=None, resolution=None,
             the ``restart_run`` step of the ``restart_test`` runs the model
             twice, the second time with ``namelist.landice.rst`` and
             ``streams.landice.rst``
+
+        albany_input_yaml : str, optional
+            The name of the Albany input yaml file in the package to use for
+            this run.
+
+        debris_friction : bool, optional
+            Whether to apply mesh field updates for the debris-friction
+            variant.
         """
 
         if suffixes is None:
             suffixes = ['landice']
         self.suffixes = suffixes
+        self.debris_friction = debris_friction
 
         # The condition below will only be true for the `SpinUp` testcase
         # where resolution is not know at the time of object construction
@@ -105,8 +118,9 @@ def __init__(self, test_case, name, subdir=None, resolution=None,
                 out_name='streams.{}'.format(suffix))
 
         self.add_input_file(filename='albany_input.yaml',
-                            package='compass.landice.tests.mismipplus',
-                            copy=True)
+                    target=albany_input_yaml,
+                    package='compass.landice.tests.mismipplus',
+                    copy=True)
 
         self.add_model_as_input()
 
@@ -178,9 +192,21 @@ def run(self):
             self.update_namelist_at_runtime(
                 {'config_ice_density': f'{ice_density}'})
 
+        if self.debris_friction:
+            self._update_mesh_for_debris_friction()
+
         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))
+
+    def _update_mesh_for_debris_friction(self):
+        with xr.open_dataset(self.mesh_file) as ds_mesh:
+            ds_mesh['muFriction'].loc[:] = 0.4
+            ds_mesh['bedRoughnessBC'] = xr.full_like(ds_mesh['muFriction'],
+                                                     6000.0)
+            ds_mesh['basalDebrisFactor'] = xr.where(
+                ds_mesh['xCell'] < 200000.0, 3.2e-2, 0.0)
+            ds_mesh.to_netcdf(self.mesh_file, mode='a')

compass/landice/tests/mismipplus/smoke_test/__init__.py

@@ -13,7 +13,7 @@ class SmokeTest(TestCase):
     spin-up that has previously been run to steady state
     """
 
-    def __init__(self, test_group, resolution):
+    def __init__(self, test_group, resolution, debris_friction=False):
         """
         Create the test case
 
@@ -25,22 +25,37 @@ def __init__(self, test_group, resolution):
         resolution : float
             The resolution of the test case. Valid options are defined in the
             test group constructor.
+
+        debris_friction : bool, optional
+            Whether to configure a variant that uses debris-friction fields
+            and Albany input options.
         """
         name = 'smoke_test'
-        subdir = f"{name}/{resolution:4.0f}m"
+        if debris_friction:
+            subdir = f"{name}/debris_friction/{resolution:4.0f}m"
+        else:
+            subdir = f"{name}/{resolution:4.0f}m"
 
         super().__init__(test_group=test_group, name=name, subdir=subdir)
 
         step_name = 'run_model'
-        step = RunModel(test_case=self, name=step_name, resolution=resolution)
+        if debris_friction:
+            albany_input_yaml = 'albany_input_debrisfriction.yaml'
+        else:
+            albany_input_yaml = 'albany_input.yaml'
+
+        step = RunModel(test_case=self, name=step_name, resolution=resolution,
+                        albany_input_yaml=albany_input_yaml,
+                        debris_friction=debris_friction)
 
         # download and link the mesh, eventually this will need to be
         # resolution aware. ``configure`` method is probably a better place
         # for parsing and adding the correct IC file based on resolution.
         step.mesh_file = 'landice_grid.nc'
         step.add_input_file(filename=step.mesh_file,
                             target='MISMIP_2km_20220502.nc',
-                            database='')
+                            database='',
+                            copy=debris_friction)
 
         self.add_step(step)