diff --git a/MALI-Dev b/MALI-Dev
index 9ef2dcae2c..fab1754df3 160000
--- a/MALI-Dev
+++ b/MALI-Dev
@@ -1 +1 @@
-Subproject commit 9ef2dcae2cbd0dcd10d4823c0b590b1e820453a8
+Subproject commit fab1754df3178bc8a1439b3ba65e0a268060152e
diff --git a/compass/landice/__init__.py b/compass/landice/__init__.py
index 338f5237a2..e22cacec79 100644
--- a/compass/landice/__init__.py
+++ b/compass/landice/__init__.py
@@ -14,6 +14,7 @@
from compass.landice.tests.isunnguata_sermia import IsunnguataSermia
from compass.landice.tests.kangerlussuaq import Kangerlussuaq
from compass.landice.tests.koge_bugt_s import KogeBugtS
+from compass.landice.tests.mesh_convergence import MeshConvergence
from compass.landice.tests.mesh_modifications import MeshModifications
from compass.landice.tests.mismipplus import MISMIPplus
from compass.landice.tests.slm import Slm
@@ -48,6 +49,7 @@ def __init__(self):
self.add_test_group(IsunnguataSermia(mpas_core=self))
self.add_test_group(Kangerlussuaq(mpas_core=self))
self.add_test_group(KogeBugtS(mpas_core=self))
+ self.add_test_group(MeshConvergence(mpas_core=self))
self.add_test_group(MeshModifications(mpas_core=self))
self.add_test_group(MISMIPplus(mpas_core=self))
self.add_test_group(Slm(mpas_core=self))
diff --git a/compass/landice/tests/dome/setup_mesh.py b/compass/landice/tests/dome/setup_mesh.py
index 8eb3c2e34d..b5937891a7 100644
--- a/compass/landice/tests/dome/setup_mesh.py
+++ b/compass/landice/tests/dome/setup_mesh.py
@@ -80,11 +80,11 @@ def run(self):
make_graph_file(mesh_filename='landice_grid.nc',
graph_filename='graph.info')
- _setup_dome_initial_conditions(config, logger,
- filename='landice_grid.nc')
+ setup_dome_initial_conditions(config, logger,
+ filename='landice_grid.nc')
-def _setup_dome_initial_conditions(config, logger, filename):
+def setup_dome_initial_conditions(config, logger, filename):
"""
Add the initial condition to the given MPAS mesh file
@@ -104,7 +104,7 @@ def _setup_dome_initial_conditions(config, logger, filename):
dome_type = section.get('dome_type')
put_origin_on_a_cell = section.getboolean('put_origin_on_a_cell')
shelf = section.getboolean('shelf')
- hydro = section.getboolean('hyrdo')
+ hydro = section.getboolean('hydro')
# Open the file, get needed dimensions
gridfile = NetCDFFile(filename, 'r+')
diff --git a/compass/landice/tests/mesh_convergence/__init__.py b/compass/landice/tests/mesh_convergence/__init__.py
new file mode 100644
index 0000000000..86a34034e8
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/__init__.py
@@ -0,0 +1,24 @@
+from compass.landice.tests.mesh_convergence.halfar import Halfar
+from compass.landice.tests.mesh_convergence.horizontal_advection import (
+ HorizontalAdvection,
+)
+from compass.landice.tests.mesh_convergence.horizontal_advection_thickness import ( # noqa
+ HorizontalAdvectionThickness,
+)
+from compass.testgroup import TestGroup
+
+
+class MeshConvergence(TestGroup):
+ """
+ A test group for convergence tests with MALI
+ """
+ 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='mesh_convergence')
+
+ self.add_test_case(HorizontalAdvection(test_group=self))
+ self.add_test_case(HorizontalAdvectionThickness(test_group=self))
+ self.add_test_case(Halfar(test_group=self))
diff --git a/compass/landice/tests/mesh_convergence/conv_analysis.py b/compass/landice/tests/mesh_convergence/conv_analysis.py
new file mode 100644
index 0000000000..e58b1b4eef
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/conv_analysis.py
@@ -0,0 +1,34 @@
+from compass.step import Step
+
+
+class ConvAnalysis(Step):
+ """
+ A step for visualizing and/or analyzing the output from a convergence test
+ case
+
+ Attributes
+ ----------
+ resolutions : list of int
+ The resolutions of the meshes that have been run
+ """
+ def __init__(self, test_case, resolutions):
+ """
+ Create the step
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+
+ resolutions : list of int
+ The resolutions of the meshes that have been run
+ """
+ super().__init__(test_case=test_case, name='analysis')
+ self.resolutions = resolutions
+
+ # typically, the analysis will rely on the output from the forward
+ # steps
+ for resolution in resolutions:
+ self.add_input_file(
+ filename='{}km_output.nc'.format(resolution),
+ target='../{}km/forward/output.nc'.format(resolution))
diff --git a/compass/landice/tests/mesh_convergence/conv_init.py b/compass/landice/tests/mesh_convergence/conv_init.py
new file mode 100644
index 0000000000..38bbf9d761
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/conv_init.py
@@ -0,0 +1,67 @@
+from mpas_tools.io import write_netcdf
+from mpas_tools.mesh.conversion import convert, cull
+from mpas_tools.planar_hex import make_planar_hex_mesh
+from mpas_tools.translate import center
+
+from compass.model import make_graph_file
+from compass.step import Step
+
+
+class ConvInit(Step):
+ """
+ A step for creating a mesh for a given resolution in a mesh convergence
+ test case. A child class of this step should then create an appropriate
+ initial condition.
+
+ Attributes
+ ----------
+ resolution : int
+ The resolution of the test case
+ """
+ def __init__(self, test_case, resolution):
+ """
+ Create the step
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+
+ resolution : int
+ The resolution of the test case
+ """
+ super().__init__(test_case=test_case,
+ name='{}km_init'.format(resolution),
+ subdir='{}km/init'.format(resolution))
+
+ for file in ['mesh.nc', 'graph.info']:
+ self.add_output_file(file)
+
+ self.resolution = resolution
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+ logger = self.logger
+ config = self.config
+ resolution = float(self.resolution)
+
+ section = config['mesh_convergence']
+ nx_1km = section.getint('nx_1km')
+ ny_1km = section.getint('ny_1km')
+ nx = int(nx_1km / resolution)
+ ny = int(ny_1km / resolution)
+ dc = resolution * 1e3
+ nonperiodic = section.getboolean('nonperiodic')
+
+ ds_mesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc,
+ nonperiodic_x=nonperiodic,
+ nonperiodic_y=nonperiodic)
+
+ ds_mesh = cull(ds_mesh, logger=logger)
+ ds_mesh = convert(ds_mesh, logger=logger)
+ center(ds_mesh)
+
+ write_netcdf(ds_mesh, 'mesh.nc')
+ make_graph_file('mesh.nc', 'graph.info')
diff --git a/compass/landice/tests/mesh_convergence/conv_test_case.py b/compass/landice/tests/mesh_convergence/conv_test_case.py
new file mode 100644
index 0000000000..1d284371d7
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/conv_test_case.py
@@ -0,0 +1,143 @@
+from compass.config import CompassConfigParser
+from compass.landice.tests.mesh_convergence.forward import Forward
+from compass.testcase import TestCase
+
+
+class ConvTestCase(TestCase):
+ """
+ A test case for various convergence tests on in MALI with planar,
+ doubly periodic meshes
+
+ Attributes
+ ----------
+ resolutions : list of int
+ """
+ def __init__(self, test_group, name):
+ """
+ Create test case for creating a MALI mesh
+
+ Parameters
+ ----------
+ test_group : compass.ocean.tests.mesh_convergence.MeshConvergence
+ The test group that this test case belongs to
+
+ name : str
+ The name of the test case
+ """
+ super().__init__(test_group=test_group, name=name)
+ self.resolutions = None
+
+ # add the steps with default resolutions so they can be listed
+ config = CompassConfigParser()
+ module = 'compass.landice.tests.mesh_convergence'
+ config.add_from_package(module, 'mesh_convergence.cfg')
+ self._setup_steps(config)
+
+ def configure(self):
+ """
+ Set config options for the test case
+ """
+ config = self.config
+ # set up the steps again in case a user has provided new resolutions
+ self._setup_steps(config)
+
+ self.update_cores()
+
+ def update_cores(self):
+ """ Update the number of cores and min_tasks for each forward step """
+
+ config = self.config
+
+ goal_cells_per_core = config.getfloat('mesh_convergence',
+ 'goal_cells_per_core')
+ max_cells_per_core = config.getfloat('mesh_convergence',
+ 'max_cells_per_core')
+
+ section = config['mesh_convergence']
+ nx_1km = section.getint('nx_1km')
+ ny_1km = section.getint('ny_1km')
+
+ for resolution in self.resolutions:
+ nx = int(nx_1km / resolution)
+ ny = int(ny_1km / resolution)
+ # a heuristic based on
+ cell_count = nx * ny
+ # ideally, about 300 cells per core
+ # (make it a multiple of 4 because...it looks better?)
+ ntasks = max(1, 4 * round(cell_count / (4 * goal_cells_per_core)))
+ # In a pinch, about 3000 cells per core
+ min_tasks = max(1, round(cell_count / max_cells_per_core))
+ step = self.steps[f'{resolution}km_forward']
+ step.ntasks = ntasks
+ step.min_tasks = min_tasks
+
+ config.set('mesh_convergence', f'{resolution}km_ntasks',
+ str(ntasks))
+ config.set('mesh_convergence', f'{resolution}km_min_tasks',
+ str(min_tasks))
+
+ def _setup_steps(self, config):
+ """
+ setup steps given resolutions
+
+ Parameters
+ ----------
+ config : compass.config.CompassConfigParser
+ The config options containing the resolutions
+ """
+
+ resolutions = config.getlist('mesh_convergence', 'resolutions',
+ dtype=int)
+
+ if self.resolutions is not None and self.resolutions == resolutions:
+ return
+
+ # start fresh with no steps
+ self.steps = dict()
+ self.steps_to_run = list()
+
+ self.resolutions = resolutions
+
+ for resolution in resolutions:
+ self.add_step(self.create_init(resolution=resolution))
+ self.add_step(Forward(test_case=self, resolution=resolution))
+
+ self.add_step(self.create_analysis(resolutions=resolutions))
+
+ def create_init(self, resolution):
+ """
+
+ Child class must override this to return an instance of a
+ ConvergenceInit step
+
+ Parameters
+ ----------
+ resolution : int
+ The resolution of the step
+
+ Returns
+ -------
+ init : compass.landice.tests.mesh_convergence.convergence_init.ConvergenceInit # noqa
+ The init step object
+ """
+
+ pass
+
+ def create_analysis(self, resolutions):
+ """
+
+ Child class must override this to return an instance of a
+ ConvergenceInit step
+
+ Parameters
+ ----------
+ resolutions : list of int
+ The resolutions of the other steps in the test case
+
+ Returns
+ -------
+ analysis : compass.landice.tests.mesh_convergence.conv_analysis.ConvAnalysis # noqa
+ The init step object
+ """
+
+ pass
diff --git a/compass/landice/tests/mesh_convergence/forward.py b/compass/landice/tests/mesh_convergence/forward.py
new file mode 100644
index 0000000000..5f83e6ba7c
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/forward.py
@@ -0,0 +1,152 @@
+import math
+from importlib.resources import contents
+
+from compass.model import run_model
+from compass.step import Step
+
+
+class Forward(Step):
+ """
+ A step for performing forward MALI runs as part of a mesh
+ convergence test case
+
+ Attributes
+ ----------
+ resolution : int
+ The resolution of the (uniform) mesh in km
+ """
+
+ def __init__(self, test_case, resolution):
+ """
+ Create a new step
+
+ Parameters
+ ----------
+ test_case : compass.landice.tests.mesh_convergence.convergence_test_case.ConvergenceTestCase
+ The test case this step belongs to
+
+ resolution : int
+ The resolution of the (uniform) mesh in km
+ """ # noqa: E501
+ super().__init__(test_case=test_case,
+ name='{}km_forward'.format(resolution),
+ subdir='{}km/forward'.format(resolution))
+
+ self.resolution = resolution
+
+ self.add_namelist_file(
+ 'compass.landice.tests.mesh_convergence', 'namelist.forward')
+
+ self.add_streams_file('compass.landice.tests.mesh_convergence',
+ 'streams.forward')
+
+ module = self.test_case.__module__
+ mesh_package_contents = list(contents(module))
+ if 'namelist.forward' in mesh_package_contents:
+ self.add_namelist_file(module, 'namelist.forward')
+ if 'streams.forward' in mesh_package_contents:
+ self.add_streams_file(module, 'streams.forward')
+
+ self.add_input_file(filename='init.nc',
+ target='../init/initial_state.nc')
+ self.add_input_file(filename='graph.info',
+ target='../init/graph.info')
+
+ self.add_model_as_input()
+
+ self.add_output_file(filename='output.nc')
+
+ def setup(self):
+ """
+ Set namelist options based on config options
+ """
+
+ namelist_options, stream_replacements = self.get_dt_duration()
+ self.add_namelist_options(namelist_options)
+
+ self.add_streams_file('compass.landice.tests.mesh_convergence',
+ 'streams.template',
+ template_replacements=stream_replacements)
+ self._get_resources()
+
+ def constrain_resources(self, available_resources):
+ """
+ Update resources at runtime from config options
+ """
+ self._get_resources()
+ super().constrain_resources(available_resources)
+
+ def run(self):
+ """
+ Run this step of the testcase
+ """
+ namelist_options, stream_replacements = self.get_dt_duration()
+ self.update_namelist_at_runtime(
+ options=namelist_options,
+ out_name='namelist.landice')
+
+ self.update_streams_at_runtime(
+ 'compass.landice.tests.mesh_convergence',
+ 'streams.template', template_replacements=stream_replacements,
+ out_name='streams.landice')
+
+ run_model(self)
+
+ def get_dt_duration(self):
+ """
+ Get the time step and run duration as namelist options from config
+ options
+
+ Returns
+ -------
+ options : dict
+ Namelist options to update
+ """
+
+ sec_in_yr = 3600.0 * 24.0 * 365.0
+
+ config = self.config
+ duration = config.getint('mesh_convergence', 'duration')
+ dur_sec = duration * sec_in_yr
+ target_velocity = config.getfloat('mesh_convergence',
+ 'target_velocity')
+ resolutions = config.getlist('mesh_convergence', 'resolutions',
+ dtype=int)
+ max_res = max(resolutions)
+ # calculate dt in s for the resolution in km and velo in m/yr
+ # apply ceil to ensure dt * nt actually exceeds duration
+ dt_max_res_cfl = math.ceil(max_res * 1000.0 / target_velocity *
+ sec_in_yr)
+ nt_max_res = math.ceil(dur_sec / dt_max_res_cfl)
+ dt_max_res = math.ceil(dur_sec / nt_max_res)
+ print(f'max res: nt={nt_max_res}, dt={dt_max_res}, '
+ f'err={abs(dur_sec - nt_max_res * dt_max_res) / dur_sec}')
+
+ dt = dt_max_res * self.resolution / max_res
+ nt = math.ceil(dur_sec / dt)
+ print(f'{self.resolution}km res: nt={nt}, dt={dt}, '
+ f'err={abs(dur_sec - nt * dt) / dur_sec}')
+
+ # check that dt*nt is acceptably close to duration
+ time_err = abs(dur_sec - nt * dt) / dur_sec
+ if time_err >= 1e-4:
+ raise ValueError(
+ f'nt*dt differs from duration by {time_err}')
+
+ # the duration (years) of the run
+ duration = f'{duration:05d}-00-00_00:00:00'
+
+ namelist_replacements = {'config_dt': f"'{dt}'",
+ 'config_run_duration': f"'{duration}'"}
+
+ stream_replacements = {'output_interval': duration}
+
+ return namelist_replacements, stream_replacements
+
+ def _get_resources(self):
+ config = self.config
+ resolution = self.resolution
+ self.ntasks = config.getint('mesh_convergence',
+ f'{resolution}km_ntasks')
+ self.min_tasks = config.getint('mesh_convergence',
+ f'{resolution}km_min_tasks')
diff --git a/compass/landice/tests/mesh_convergence/halfar/__init__.py b/compass/landice/tests/mesh_convergence/halfar/__init__.py
new file mode 100644
index 0000000000..215e118fa8
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/halfar/__init__.py
@@ -0,0 +1,60 @@
+from compass.landice.tests.mesh_convergence.conv_test_case import ConvTestCase
+from compass.landice.tests.mesh_convergence.halfar.analysis import ( # noqa
+ Analysis,
+)
+from compass.landice.tests.mesh_convergence.halfar.init import Init
+
+
+class Halfar(ConvTestCase):
+ """
+ A test case for testing mesh convergence with the Halfar analytic test
+ """
+ def __init__(self, test_group):
+ """
+ Create test case for creating a MALI mesh
+
+ Parameters
+ ----------
+ test_group : compass.landice.tests.mesh_convergence.MeshConvergence
+ The landice test group that this test case belongs to
+ """
+ super().__init__(test_group=test_group, name='halfar')
+
+ self.add_step(Analysis(test_case=self, resolutions=self.resolutions))
+
+ def create_init(self, resolution):
+ """
+ Child class must override this to return an instance of a
+ ConvInit step
+
+ Parameters
+ ----------
+ resolution : int
+ The resolution of the test case
+
+ Returns
+ -------
+ init : compass.landice.tests.mesh_convergence.conv_init.ConvInit
+ The init step object
+ """
+
+ return Init(test_case=self, resolution=resolution)
+
+ def create_analysis(self, resolutions):
+ """
+
+ Child class must override this to return an instance of a
+ ConvergenceInit step
+
+ Parameters
+ ----------
+ resolutions : list of int
+ The resolutions of the other steps in the test case
+
+ Returns
+ -------
+ analysis : compass.landice.tests.mesh_convergence.conv_analysis.ConvAnalysis # noqa
+ The init step object
+ """
+
+ return Analysis(test_case=self, resolutions=resolutions)
diff --git a/compass/landice/tests/mesh_convergence/halfar/analysis.py b/compass/landice/tests/mesh_convergence/halfar/analysis.py
new file mode 100644
index 0000000000..5c7ac81c49
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/halfar/analysis.py
@@ -0,0 +1,190 @@
+import warnings
+
+import matplotlib.pyplot as plt
+import numpy as np
+import xarray as xr
+
+from compass.landice.tests.mesh_convergence.conv_analysis import ConvAnalysis
+
+
+class Analysis(ConvAnalysis):
+ """
+ A step for visualizing the output from the advection convergence test case
+ """
+ def __init__(self, test_case, resolutions):
+ """
+ Create the step
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+
+ resolutions : list of int
+ The resolutions of the meshes that have been run
+ """
+ super().__init__(test_case=test_case, resolutions=resolutions)
+ self.resolutions = resolutions
+ self.add_output_file('convergence_rmse.png')
+ self.add_output_file('convergence_dome.png')
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+ plt.switch_backend('Agg')
+ resolutions = self.resolutions
+ ncells_list = list()
+ rmse_errors = list()
+ center_errors = list()
+ for res in resolutions:
+ rms_error, center_error, ncells = self.rmse(res)
+ ncells_list.append(ncells)
+ rmse_errors.append(rms_error)
+ center_errors.append(center_error)
+
+ ncells = np.array(ncells_list)
+ rmse_errors = np.array(rmse_errors)
+ center_errors = np.array(center_errors)
+
+ # plot rmse errors
+ p = np.polyfit(np.log10(ncells), np.log10(rmse_errors), 1)
+ conv = abs(p[0]) * 2.0
+ error_fit = ncells**p[0] * 10**p[1]
+
+ plt.figure()
+ plt.loglog(ncells, error_fit, 'k')
+ plt.loglog(ncells, rmse_errors, 'or')
+ plt.annotate('Order of Convergence = {}'.format(np.round(conv, 3)),
+ xycoords='axes fraction', xy=(0.3, 0.95), fontsize=14)
+ plt.xlabel('Number of Grid Cells', fontsize=14)
+ plt.ylabel('L2 Norm', fontsize=14)
+ section = self.config['mesh_convergence']
+ duration = section.getfloat('duration')
+ plt.title(f'Halfar convergence test, {duration} yrs')
+ plt.savefig('convergence_rmse.png', bbox_inches='tight',
+ pad_inches=0.1)
+
+ # now repeat for center errors
+ plt.figure()
+ p = np.polyfit(np.log10(ncells), np.log10(center_errors), 1)
+ conv2 = abs(p[0]) * 2.0
+ error_fit = ncells**p[0] * 10**p[1]
+
+ plt.loglog(ncells, error_fit, 'k')
+ plt.loglog(ncells, center_errors, 'or')
+ plt.annotate('Order of Convergence = {}'.format(np.round(conv2, 3)),
+ xycoords='axes fraction', xy=(0.3, 0.95), fontsize=14)
+ plt.xlabel('Number of Grid Cells', fontsize=14)
+ plt.ylabel('Dome center error (m)', fontsize=14)
+ plt.title(f'Halfar convergence test, {duration} yrs')
+ plt.savefig('convergence_dome.png', bbox_inches='tight',
+ pad_inches=0.1)
+
+ section = self.config['halfar']
+ conv_thresh = section.getfloat('conv_thresh')
+ conv_max = section.getfloat('conv_max')
+
+ if conv < conv_thresh:
+ raise ValueError(f'order of convergence '
+ f'{conv} < min tolerance {conv_thresh}')
+
+ if conv > conv_max:
+ warnings.warn(f'order of convergence '
+ f'{conv} > max tolerance {conv_max}')
+
+ def rmse(self, resolution):
+ """
+ Compute the RMSE for a given resolution
+
+ Parameters
+ ----------
+ resolution : int
+ The resolution of the (uniform) mesh in km
+
+ Returns
+ -------
+ rms_error : float
+ The root-mean-squared error
+
+ ncells : int
+ The number of cells in the mesh
+ """
+ res_tag = '{}km'.format(resolution)
+
+ timelev = -1 # todo: determine a specified time level and err check
+
+ ds = xr.open_dataset('{}_output.nc'.format(res_tag),
+ decode_cf=False,
+ decode_timedelta=False)
+ # Find out what the ice density and flowA values for this run were.
+ print('Collecting parameter values from the output file.')
+ flowA = float(ds.config_default_flowParamA)
+ print(f'Using a flowParamA value of: {flowA}')
+ flow_n = float(ds.config_flowLawExponent)
+ print(f'Using a flowLawExponent value of: {flow_n}')
+ if flow_n != 3:
+ raise ValueError('Error: The Halfar script currently only '
+ 'supports a flow law exponent of 3.')
+ rhoi = ds.config_ice_density
+ print(f'Using an ice density value of: {rhoi}')
+ dynamicThickness = float(ds.config_dynamic_thickness)
+ print(f'Dynamic thickness for this run = {dynamicThickness}')
+ daysSinceStart = ds.daysSinceStart
+ print(f'Using model time of {daysSinceStart / 365.0} years')
+ if ds.config_calendar_type != "noleap":
+ raise ValueError('Error: The Halfar script currently assumes a '
+ 'noleap calendar.')
+
+ ncells = ds.sizes['nCells']
+ thk = ds['thickness'].isel(Time=timelev)
+ xCell = ds['xCell'].values
+ yCell = ds['yCell'].values
+ areaCell = ds['areaCell'].values
+
+ dt = (daysSinceStart[timelev] - daysSinceStart[0]) * 3600.0 * 24.0
+
+ thkHalfar = _halfar(dt, xCell, yCell, flowA, flow_n, rhoi)
+ mask = np.where(np.logical_or(thk > 0.0, thkHalfar > 0.0))
+ diff = thk - thkHalfar
+ rms_error = ((diff[mask]**2 * areaCell[mask] /
+ areaCell[mask].sum()).sum())**0.5
+
+ center_idx = np.where((xCell == 0.0) * (yCell == 0.0))
+ center_error = float(np.asarray(np.absolute(diff[center_idx])).item())
+
+ return rms_error, center_error, ncells
+
+
+def _halfar(t, x, y, A, n, rho):
+ # A # s^{-1} Pa^{-3}
+ # n # Glen flow law exponent
+ # rho # ice density kg m^{-3}
+
+ # These constants should come from setup_dome_initial_conditions.py.
+ # For now they are being hardcoded.
+ R0 = 60000.0 * np.sqrt(0.125) # initial dome radius
+ H0 = 2000.0 * np.sqrt(0.125) # initial dome thickness at center
+ g = 9.80616 # gravity m/s/s -- value used by mpas_constants
+ alpha = 1.0 / 9.0
+ beta = 1.0 / 18.0
+ Gamma = 2.0 / (n + 2.0) * A * (rho * g)**n
+
+ # The IC file puts the center of the dome and the domain at (0,0)
+ x0 = 0.0
+ y0 = 0.0
+
+ # NOTE: These constants assume n=3
+ # they need to be generalized to allow other n's
+ t0 = (beta / Gamma) * (7.0 / 4.0)**3 * (R0**4 / H0**7)
+ t = t + t0
+ t = t / t0
+
+ H = np.zeros(len(x))
+ for i in range(len(x)):
+ r = np.sqrt((x[i] - x0)**2 + (y[i] - y0)**2)
+ r = r / R0
+ inside = max(0.0, 1.0 - (r / t**beta)**((n + 1.0) / n))
+
+ H[i] = H0 * inside**(n / (2.0 * n + 1.0)) / t**alpha
+ return H
diff --git a/compass/landice/tests/mesh_convergence/halfar/halfar.cfg b/compass/landice/tests/mesh_convergence/halfar/halfar.cfg
new file mode 100644
index 0000000000..05f02b253e
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/halfar/halfar.cfg
@@ -0,0 +1,45 @@
+# options for halfar mesh convergence test
+[halfar]
+
+# Number of vertical levels
+vert_levels = 10
+
+# convergence threshold below which the test fails
+conv_thresh = 0.0
+
+# Convergence rate above which a warning is issued
+conv_max = 3.0
+
+# config options for dome test cases
+[dome]
+
+# the dome type ('halfar' or 'cism')
+dome_type = halfar
+
+# Whether to center the dome in the center of the cell that is closest to the
+# center of the domain
+put_origin_on_a_cell = True
+
+# whether to add a small shelf to the test
+shelf = False
+
+# whether to add hydrology to the initial condition
+hydro = False
+
+[mesh_convergence]
+
+# number of mesh cells in x and y for 1 km resolution. Other resolutions have
+# the same physical size.
+# must be divisible by the ratio of the other meshes resolutions to 1 km
+nx_1km = 128
+ny_1km = 128
+nonperiodic = True
+
+# target velocity (m/yr) to use for setting a time step based on the mesh resolution
+# and an advective CFL condition.
+# Note that Halfar is subject to a more restrictive diffusive CFL so this must be
+# artificially inflated
+target_velocity = 30000.0
+
+# the duration (years) of the run
+duration = 200
diff --git a/compass/landice/tests/mesh_convergence/halfar/init.py b/compass/landice/tests/mesh_convergence/halfar/init.py
new file mode 100644
index 0000000000..89217897a0
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/halfar/init.py
@@ -0,0 +1,57 @@
+import numpy
+import xarray
+from mpas_tools.io import write_netcdf
+
+from compass.landice.tests.dome.setup_mesh import setup_dome_initial_conditions
+from compass.landice.tests.mesh_convergence.conv_init import ConvInit
+
+
+class Init(ConvInit):
+ """
+ A step for creating an initial_condition for advection convergence test
+ case
+ """
+ def __init__(self, test_case, resolution):
+ """
+ Create the step
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+
+ resolution : int
+ The resolution of the test case
+ """
+ super().__init__(test_case=test_case, resolution=resolution)
+ self.add_output_file('initial_state.nc')
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+ # create the mesh and graph.info
+ super().run()
+
+ config = self.config
+ logger = self.logger
+ filename = 'initial_state.nc'
+
+ section = config['halfar']
+ nVertLevels = section.getint('vert_levels')
+
+ ds = xarray.open_dataset('mesh.nc')
+ xCell = ds.xCell
+
+ layerThicknessFractions = xarray.DataArray(
+ data=1.0 / nVertLevels * numpy.ones((nVertLevels,)),
+ dims=['nVertLevels'])
+ ds['layerThicknessFractions'] = layerThicknessFractions
+ thickness = xarray.zeros_like(xCell)
+ thickness = thickness.expand_dims(dim='Time', axis=0)
+ ds['thickness'] = thickness
+ ds['bedTopography'] = xarray.zeros_like(thickness)
+ ds['sfcMassBal'] = xarray.zeros_like(thickness)
+ write_netcdf(ds, filename)
+
+ setup_dome_initial_conditions(config, logger, filename)
diff --git a/compass/landice/tests/mesh_convergence/halfar/namelist.forward b/compass/landice/tests/mesh_convergence/halfar/namelist.forward
new file mode 100644
index 0000000000..a4e25ac67f
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/halfar/namelist.forward
@@ -0,0 +1,9 @@
+config_velocity_solver = 'sia'
+config_thickness_advection = 'fct'
+config_tracer_advection = 'fct'
+config_horiz_tracer_adv_order = 3
+config_advection_coef_3rd_order = 1.0
+config_start_time = '0001-01-01_00:00:00'
+config_time_integration = 'runge_kutta'
+config_rk_order = 3
+config_rk3_stages = 3
diff --git a/compass/landice/tests/mesh_convergence/halfar/streams.forward b/compass/landice/tests/mesh_convergence/halfar/streams.forward
new file mode 100644
index 0000000000..25123a88f7
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/halfar/streams.forward
@@ -0,0 +1,7 @@
+
+
+
+
+
+
+
diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/__init__.py b/compass/landice/tests/mesh_convergence/horizontal_advection/__init__.py
new file mode 100644
index 0000000000..094c426fd0
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection/__init__.py
@@ -0,0 +1,63 @@
+from compass.landice.tests.mesh_convergence.conv_test_case import ConvTestCase
+from compass.landice.tests.mesh_convergence.horizontal_advection.analysis import ( # noqa
+ Analysis,
+)
+from compass.landice.tests.mesh_convergence.horizontal_advection.init import (
+ Init,
+)
+
+
+class HorizontalAdvection(ConvTestCase):
+ """
+ A test case for testing horizontal advection in MALI with planar,
+ doubly periodic meshes
+ """
+ def __init__(self, test_group):
+ """
+ Create test case for creating a MALI mesh
+
+ Parameters
+ ----------
+ test_group : compass.landice.tests.mesh_convergence.MeshConvergence
+ The landice test group that this test case belongs to
+ """
+ super().__init__(test_group=test_group, name='horizontal_advection')
+
+ self.add_step(Analysis(test_case=self, resolutions=self.resolutions))
+
+ def create_init(self, resolution):
+ """
+ Child class must override this to return an instance of a
+ ConvInit step
+
+ Parameters
+ ----------
+ resolution : int
+ The resolution of the test case
+
+ Returns
+ -------
+ init : compass.landice.tests.mesh_convergence.conv_init.ConvInit
+ The init step object
+ """
+
+ return Init(test_case=self, resolution=resolution)
+
+ def create_analysis(self, resolutions):
+ """
+
+ Child class must override this to return an instance of a
+ ConvergenceInit step
+
+ Parameters
+ ----------
+ resolutions : list of int
+ The resolutions of the other steps in the test case
+
+ Returns
+ -------
+ analysis : compass.landice.tests.mesh_convergence.conv_analysis.ConvAnalysis # noqa
+ The init step object
+ """
+
+ return Analysis(test_case=self, resolutions=resolutions)
diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py
new file mode 100644
index 0000000000..377cbcef7e
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py
@@ -0,0 +1,104 @@
+import warnings
+
+import matplotlib.pyplot as plt
+import numpy as np
+import xarray as xr
+
+from compass.landice.tests.mesh_convergence.conv_analysis import ConvAnalysis
+
+
+class Analysis(ConvAnalysis):
+ """
+ A step for visualizing the output from the advection convergence test case
+ """
+ def __init__(self, test_case, resolutions):
+ """
+ Create the step
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+
+ resolutions : list of int
+ The resolutions of the meshes that have been run
+ """
+ super().__init__(test_case=test_case, resolutions=resolutions)
+ self.resolutions = resolutions
+ self.add_output_file('convergence.png')
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+ plt.switch_backend('Agg')
+ resolutions = self.resolutions
+ ncells_list = list()
+ errors = list()
+ for res in resolutions:
+ rms_error, ncells = self.rmse(res, variable='passiveTracer2d')
+ ncells_list.append(ncells)
+ errors.append(rms_error)
+
+ ncells = np.array(ncells_list)
+ errors = np.array(errors)
+
+ p = np.polyfit(np.log10(ncells), np.log10(errors), 1)
+ conv = abs(p[0]) * 2.0
+
+ error_fit = ncells**p[0] * 10**p[1]
+
+ plt.figure()
+ plt.loglog(ncells, error_fit, 'k')
+ plt.loglog(ncells, errors, 'or')
+ plt.annotate('Order of Convergence = {}'.format(np.round(conv, 3)),
+ xycoords='axes fraction', xy=(0.3, 0.95), fontsize=14)
+ plt.xlabel('Number of Grid Cells', fontsize=14)
+ plt.ylabel('L2 Norm', fontsize=14)
+ section = self.config['mesh_convergence']
+ duration = section.getfloat('duration')
+ plt.title(f'Horizontal advection convergence test, {duration} yrs')
+ plt.savefig('convergence.png', bbox_inches='tight', pad_inches=0.1)
+
+ section = self.config['horizontal_advection']
+ conv_thresh = section.getfloat('conv_thresh')
+ conv_max = section.getfloat('conv_max')
+
+ if conv < conv_thresh:
+ raise ValueError(f'order of convergence '
+ f'{conv} < min tolerance {conv_thresh}')
+
+ if conv > conv_max:
+ warnings.warn(f'order of convergence '
+ f'{conv} > max tolerance {conv_max}')
+
+ def rmse(self, resolution, variable):
+ """
+ Compute the RMSE for a given resolution
+
+ Parameters
+ ----------
+ resolution : int
+ The resolution of the (uniform) mesh in km
+
+ variable : str
+ The name of a variable in the output file to analyze.
+
+ Returns
+ -------
+ rms_error : float
+ The root-mean-squared error
+
+ ncells : int
+ The number of cells in the mesh
+ """
+ res_tag = '{}km'.format(resolution)
+
+ ds = xr.open_dataset('{}_output.nc'.format(res_tag),
+ decode_timedelta=False)
+ init = ds[variable].isel(Time=0)
+ final = ds[variable].isel(Time=-1)
+ diff = final - init
+ rms_error = np.sqrt((diff**2).mean()).values
+ ncells = ds.sizes['nCells']
+ return rms_error, ncells
diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/horizontal_advection.cfg b/compass/landice/tests/mesh_convergence/horizontal_advection/horizontal_advection.cfg
new file mode 100644
index 0000000000..22bc39cbe4
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection/horizontal_advection.cfg
@@ -0,0 +1,34 @@
+# options for planar horizontal advection test case
+[horizontal_advection]
+
+# Number of vertical levels
+vert_levels = 3
+
+# ice thickness (m)
+ice_thickness = 1000.0
+
+# bed elevation (m)
+bed_elevation = 0.0
+
+# center of the tracer gaussian (km)
+x_center = 0.
+y_center = 0.
+
+# width of gaussian tracer "blob" (km)
+gaussian_width = 50
+
+# whether to advect in x, y, or both
+advect_x = True
+advect_y = True
+
+# convergence threshold below which the test fails
+conv_thresh = 0.6
+
+# Convergence rate above which a warning is issued
+conv_max = 3.0
+
+# options for mesh convergence test cases
+[mesh_convergence]
+
+# a list of resolutions (km) to test
+resolutions = 16, 8, 4, 2
diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/init.py b/compass/landice/tests/mesh_convergence/horizontal_advection/init.py
new file mode 100644
index 0000000000..7f335a15e5
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection/init.py
@@ -0,0 +1,97 @@
+import numpy
+import xarray
+from mpas_tools.io import write_netcdf
+
+from compass.landice.tests.mesh_convergence.conv_init import ConvInit
+
+
+class Init(ConvInit):
+ """
+ A step for creating an initial_condition for advection convergence test
+ case
+ """
+ def __init__(self, test_case, resolution):
+ """
+ Create the step
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+
+ resolution : int
+ The resolution of the test case
+ """
+ super().__init__(test_case=test_case, resolution=resolution)
+ self.add_output_file('initial_state.nc')
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+ # create the mesh and graph.info
+ super().run()
+
+ config = self.config
+
+ section = config['horizontal_advection']
+ ice_thickness = section.getfloat('ice_thickness')
+ bed_elevation = section.getfloat('bed_elevation')
+ x_center = 1e3 * section.getfloat('x_center')
+ y_center = 1e3 * section.getfloat('y_center')
+ advect_x = section.getboolean('advect_x')
+ advect_y = section.getboolean('advect_y')
+ gaussian_width = 1e3 * section.getfloat('gaussian_width')
+ nVertLevels = section.getint('vert_levels')
+
+ section = config['mesh_convergence']
+ duration = section.getfloat('duration') * 3600.0 * 24.0 * 365.0
+
+ ds = xarray.open_dataset('mesh.nc')
+ xCell = ds.xCell
+ yCell = ds.yCell
+
+ if advect_x:
+ x_vel = ds.attrs['x_period'] / duration
+ else:
+ x_vel = 0.
+
+ if advect_y:
+ y_vel = ds.attrs['y_period'] / duration
+ else:
+ y_vel = 0.
+
+ layerThicknessFractions = xarray.DataArray(
+ data=1.0 / nVertLevels * numpy.ones((nVertLevels,)),
+ dims=['nVertLevels'])
+
+ thickness = ice_thickness * xarray.ones_like(xCell)
+ thickness = thickness.expand_dims(dim='Time', axis=0)
+
+ bedTopography = bed_elevation * xarray.ones_like(thickness)
+
+ uReconstructX = x_vel * xarray.ones_like(xCell)
+ uReconstructX = uReconstructX.expand_dims(dim={"nVertInterfaces":
+ nVertLevels + 1},
+ axis=1)
+ uReconstructX = uReconstructX.expand_dims(dim='Time', axis=0)
+
+ uReconstructY = y_vel * xarray.ones_like(xCell)
+ uReconstructY = uReconstructY.expand_dims(dim={"nVertInterfaces":
+ nVertLevels + 1},
+ axis=1)
+ uReconstructY = uReconstructY.expand_dims(dim='Time', axis=0)
+
+ dist_sq = (xCell - x_center)**2 + (yCell - y_center)**2
+
+ passiveTracer = numpy.exp(-0.5 * dist_sq / gaussian_width**2)
+ passiveTracer = passiveTracer.expand_dims(dim='Time', axis=0)
+
+ ds['layerThicknessFractions'] = layerThicknessFractions
+ ds['thickness'] = thickness
+ ds['bedTopography'] = bedTopography
+ ds['uReconstructX'] = uReconstructX
+ ds['uReconstructY'] = uReconstructY
+ ds['passiveTracer2d'] = passiveTracer
+
+ write_netcdf(ds, 'initial_state.nc')
diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/namelist.forward b/compass/landice/tests/mesh_convergence/horizontal_advection/namelist.forward
new file mode 100644
index 0000000000..66a39e998c
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection/namelist.forward
@@ -0,0 +1,10 @@
+config_velocity_solver = 'none'
+config_unrealistic_velocity = 1.0e36
+config_thickness_advection = 'fct'
+config_horiz_tracer_adv_order = 3
+config_advection_coef_3rd_order = 1.0
+config_tracer_advection = 'fct'
+config_start_time = '0001-01-01_00:00:00'
+config_time_integration = 'runge_kutta'
+config_rk_order = 3
+config_rk3_stages = 3
diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/streams.forward b/compass/landice/tests/mesh_convergence/horizontal_advection/streams.forward
new file mode 100644
index 0000000000..9f5d29f952
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection/streams.forward
@@ -0,0 +1,8 @@
+
+
+
+
+
+
+
+
diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/__init__.py b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/__init__.py
new file mode 100644
index 0000000000..2007bd986a
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/__init__.py
@@ -0,0 +1,64 @@
+from compass.landice.tests.mesh_convergence.conv_test_case import ConvTestCase
+from compass.landice.tests.mesh_convergence.horizontal_advection_thickness.analysis import ( # noqa
+ Analysis,
+)
+from compass.landice.tests.mesh_convergence.horizontal_advection_thickness.init import ( # noqa
+ Init,
+)
+
+
+class HorizontalAdvectionThickness(ConvTestCase):
+ """
+ A test case for testing horizontal advection of thickness in MALI with
+ planar, doubly periodic meshes
+ """
+ def __init__(self, test_group):
+ """
+ Create test case for creating a MALI mesh
+
+ Parameters
+ ----------
+ test_group : compass.landice.tests.mesh_convergence.MeshConvergence
+ The landice test group that this test case belongs to
+ """
+ super().__init__(test_group=test_group,
+ name='horizontal_advection_thickness')
+
+ self.add_step(Analysis(test_case=self, resolutions=self.resolutions))
+
+ def create_init(self, resolution):
+ """
+ Child class must override this to return an instance of a
+ ConvInit step
+
+ Parameters
+ ----------
+ resolution : int
+ The resolution of the test case
+
+ Returns
+ -------
+ init : compass.landice.tests.mesh_convergence.conv_init.ConvInit
+ The init step object
+ """
+
+ return Init(test_case=self, resolution=resolution)
+
+ def create_analysis(self, resolutions):
+ """
+
+ Child class must override this to return an instance of a
+ ConvergenceInit step
+
+ Parameters
+ ----------
+ resolutions : list of int
+ The resolutions of the other steps in the test case
+
+ Returns
+ -------
+ analysis : compass.landice.tests.mesh_convergence.conv_analysis.ConvAnalysis # noqa
+ The init step object
+ """
+
+ return Analysis(test_case=self, resolutions=resolutions)
diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py
new file mode 100644
index 0000000000..44de4b0e2e
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py
@@ -0,0 +1,105 @@
+import warnings
+
+import matplotlib.pyplot as plt
+import numpy as np
+import xarray as xr
+
+from compass.landice.tests.mesh_convergence.conv_analysis import ConvAnalysis
+
+
+class Analysis(ConvAnalysis):
+ """
+ A step for visualizing the output from the advection convergence test case
+ """
+ def __init__(self, test_case, resolutions):
+ """
+ Create the step
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+
+ resolutions : list of int
+ The resolutions of the meshes that have been run
+ """
+ super().__init__(test_case=test_case, resolutions=resolutions)
+ self.resolutions = resolutions
+ self.add_output_file('convergence.png')
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+ plt.switch_backend('Agg')
+ resolutions = self.resolutions
+ ncells_list = list()
+ errors = list()
+ for res in resolutions:
+ rms_error, ncells = self.rmse(res, variable='thickness')
+ ncells_list.append(ncells)
+ errors.append(rms_error)
+
+ ncells = np.array(ncells_list)
+ errors = np.array(errors)
+
+ p = np.polyfit(np.log10(ncells), np.log10(errors), 1)
+ conv = abs(p[0]) * 2.0
+
+ error_fit = ncells**p[0] * 10**p[1]
+
+ plt.figure()
+ plt.loglog(ncells, error_fit, 'k')
+ plt.loglog(ncells, errors, 'or')
+ plt.annotate('Order of Convergence = {}'.format(np.round(conv, 3)),
+ xycoords='axes fraction', xy=(0.3, 0.95), fontsize=14)
+ plt.xlabel('Number of Grid Cells', fontsize=14)
+ plt.ylabel('L2 Norm', fontsize=14)
+ section = self.config['mesh_convergence']
+ duration = section.getfloat('duration')
+ plt.title(f'Thickness horizontal advection convergence test, '
+ f'{duration} yrs')
+ plt.savefig('convergence.png', bbox_inches='tight', pad_inches=0.1)
+
+ section = self.config['horizontal_advection']
+ conv_thresh = section.getfloat('conv_thresh')
+ conv_max = section.getfloat('conv_max')
+
+ if conv < conv_thresh:
+ raise ValueError(f'order of convergence '
+ f'{conv} < min tolerance {conv_thresh}')
+
+ if conv > conv_max:
+ warnings.warn(f'order of convergence '
+ f'{conv} > max tolerance {conv_max}')
+
+ def rmse(self, resolution, variable):
+ """
+ Compute the RMSE for a given resolution
+
+ Parameters
+ ----------
+ resolution : int
+ The resolution of the (uniform) mesh in km
+
+ variable : str
+ The name of a variable in the output file to analyze.
+
+ Returns
+ -------
+ rms_error : float
+ The root-mean-squared error
+
+ ncells : int
+ The number of cells in the mesh
+ """
+ res_tag = '{}km'.format(resolution)
+
+ ds = xr.open_dataset('{}_output.nc'.format(res_tag),
+ decode_timedelta=False)
+ init = ds[variable].isel(Time=0)
+ final = ds[variable].isel(Time=-1)
+ diff = final - init
+ rms_error = np.sqrt((diff**2).mean()).values
+ ncells = ds.sizes['nCells']
+ return rms_error, ncells
diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/horizontal_advection_thickness.cfg b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/horizontal_advection_thickness.cfg
new file mode 100644
index 0000000000..22bc39cbe4
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/horizontal_advection_thickness.cfg
@@ -0,0 +1,34 @@
+# options for planar horizontal advection test case
+[horizontal_advection]
+
+# Number of vertical levels
+vert_levels = 3
+
+# ice thickness (m)
+ice_thickness = 1000.0
+
+# bed elevation (m)
+bed_elevation = 0.0
+
+# center of the tracer gaussian (km)
+x_center = 0.
+y_center = 0.
+
+# width of gaussian tracer "blob" (km)
+gaussian_width = 50
+
+# whether to advect in x, y, or both
+advect_x = True
+advect_y = True
+
+# convergence threshold below which the test fails
+conv_thresh = 0.6
+
+# Convergence rate above which a warning is issued
+conv_max = 3.0
+
+# options for mesh convergence test cases
+[mesh_convergence]
+
+# a list of resolutions (km) to test
+resolutions = 16, 8, 4, 2
diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/init.py b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/init.py
new file mode 100644
index 0000000000..bb86645cf3
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/init.py
@@ -0,0 +1,98 @@
+import numpy
+import xarray
+from mpas_tools.io import write_netcdf
+
+from compass.landice.tests.mesh_convergence.conv_init import ConvInit
+
+
+class Init(ConvInit):
+ """
+ A step for creating an initial_condition for advection convergence test
+ case
+ """
+ def __init__(self, test_case, resolution):
+ """
+ Create the step
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+
+ resolution : int
+ The resolution of the test case
+ """
+ super().__init__(test_case=test_case, resolution=resolution)
+ self.add_output_file('initial_state.nc')
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+ # create the mesh and graph.info
+ super().run()
+
+ config = self.config
+
+ section = config['horizontal_advection']
+ ice_thickness = section.getfloat('ice_thickness')
+ bed_elevation = section.getfloat('bed_elevation')
+ x_center = 1e3 * section.getfloat('x_center')
+ y_center = 1e3 * section.getfloat('y_center')
+ advect_x = section.getboolean('advect_x')
+ advect_y = section.getboolean('advect_y')
+ gaussian_width = 1e3 * section.getfloat('gaussian_width')
+ nVertLevels = section.getint('vert_levels')
+
+ section = config['mesh_convergence']
+ duration = section.getfloat('duration') * 3600.0 * 24.0 * 365.0
+
+ ds = xarray.open_dataset('mesh.nc')
+ xCell = ds.xCell
+ yCell = ds.yCell
+
+ if advect_x:
+ x_vel = ds.attrs['x_period'] / duration
+ else:
+ x_vel = 0.
+
+ if advect_y:
+ y_vel = ds.attrs['y_period'] / duration
+ else:
+ y_vel = 0.
+
+ layerThicknessFractions = xarray.DataArray(
+ data=1.0 / nVertLevels * numpy.ones((nVertLevels,)),
+ dims=['nVertLevels'])
+
+ # create base layer of uniform thickness over entire domain
+ # to ensure no margins
+ thickness = ice_thickness * xarray.ones_like(xCell)
+ # now add Gaussian bump on top
+ # using bump height equal to base height
+ dist_sq = (xCell - x_center)**2 + (yCell - y_center)**2
+ thickness += ice_thickness * numpy.exp(-0.5 * dist_sq /
+ gaussian_width**2)
+ thickness = thickness.expand_dims(dim='Time', axis=0)
+
+ bedTopography = bed_elevation * xarray.ones_like(thickness)
+
+ uReconstructX = x_vel * xarray.ones_like(xCell)
+ uReconstructX = uReconstructX.expand_dims(dim={"nVertInterfaces":
+ nVertLevels + 1},
+ axis=1)
+ uReconstructX = uReconstructX.expand_dims(dim='Time', axis=0)
+
+ uReconstructY = y_vel * xarray.ones_like(xCell)
+ uReconstructY = uReconstructY.expand_dims(dim={"nVertInterfaces":
+ nVertLevels + 1},
+ axis=1)
+ uReconstructY = uReconstructY.expand_dims(dim='Time', axis=0)
+
+ ds['layerThicknessFractions'] = layerThicknessFractions
+ ds['thickness'] = thickness
+ ds['bedTopography'] = bedTopography
+ ds['uReconstructX'] = uReconstructX
+ ds['uReconstructY'] = uReconstructY
+
+ write_netcdf(ds, 'initial_state.nc')
diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/namelist.forward b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/namelist.forward
new file mode 100644
index 0000000000..82767ef85a
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/namelist.forward
@@ -0,0 +1,10 @@
+config_velocity_solver = 'none'
+config_unrealistic_velocity = 1.0e36
+config_thickness_advection = 'fct'
+config_horiz_tracer_adv_order = 3
+config_advection_coef_3rd_order = 1.0
+config_tracer_advection = 'none'
+config_start_time = '0001-01-01_00:00:00'
+config_time_integration = 'runge_kutta'
+config_rk_order = 3
+config_rk3_stages = 3
diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/streams.forward b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/streams.forward
new file mode 100644
index 0000000000..25123a88f7
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/streams.forward
@@ -0,0 +1,7 @@
+
+
+
+
+
+
+
diff --git a/compass/landice/tests/mesh_convergence/mesh_convergence.cfg b/compass/landice/tests/mesh_convergence/mesh_convergence.cfg
new file mode 100644
index 0000000000..f07310f108
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/mesh_convergence.cfg
@@ -0,0 +1,30 @@
+# options for mesh convergence test cases
+[mesh_convergence]
+
+# a list of resolutions (km) to test
+resolutions = 8, 4, 2, 1
+
+# number of mesh cells in x and y for 1 km resolution. Other resolutions have
+# the same physical size. The default is approximately square, because of the
+# staggering of the hex mesh.
+# The integers used need to be divisible by the ratio of lowest to highest resolutions
+nx_1km = 512
+ny_1km = 640
+
+# whether to make mesh nonperiodic or not
+nonperiodic = False
+
+# 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 = 3000
+
+# target velocity (m/yr) to be used for defining the timestep
+# the actual velocity will be determined by the test case
+# This should be slightly larger than the maximum velocity expected
+target_velocity = 2000.0
+
+# the duration (years) of the run
+duration = 1000
diff --git a/compass/landice/tests/mesh_convergence/namelist.forward b/compass/landice/tests/mesh_convergence/namelist.forward
new file mode 100644
index 0000000000..a7970f6389
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/namelist.forward
@@ -0,0 +1,5 @@
+config_velocity_solver = 'none'
+config_unrealistic_velocity = 1.0e36
+config_thickness_advection = 'fo'
+config_tracer_advection = 'fo'
+config_start_time = '0001-01-01_00:00:00'
diff --git a/compass/landice/tests/mesh_convergence/streams.forward b/compass/landice/tests/mesh_convergence/streams.forward
new file mode 100644
index 0000000000..78bf23f620
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/streams.forward
@@ -0,0 +1,20 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/compass/landice/tests/mesh_convergence/streams.template b/compass/landice/tests/mesh_convergence/streams.template
new file mode 100644
index 0000000000..16d4905de5
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/streams.template
@@ -0,0 +1,6 @@
+
+
+
+
+
diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst
index 6be0af6e49..ab032fc7ea 100644
--- a/docs/developers_guide/landice/api.rst
+++ b/docs/developers_guide/landice/api.rst
@@ -417,6 +417,64 @@ koge_bugt_s
mesh.Mesh
mesh.Mesh.run
+mesh_convergence
+~~~~~~~~~~~~~~~~
+
+.. currentmodule:: compass.landice.tests.mesh_convergence
+
+.. autosummary::
+ :toctree: generated/
+
+ MeshConvergence
+
+ conv_test_case.ConvTestCase
+ conv_test_case.ConvTestCase.configure
+ conv_test_case.ConvTestCase.update_cores
+
+ conv_init.ConvInit
+ conv_init.ConvInit.run
+
+ conv_analysis.ConvAnalysis
+
+ forward.Forward
+ forward.Forward.setup
+ forward.Forward.constrain_resources
+ forward.Forward.run
+ forward.Forward.get_dt_duration
+
+ halfar.Halfar
+ halfar.Halfar.create_init
+ halfar.Halfar.create_analysis
+
+ halfar.init.Init
+ halfar.init.Init.run
+
+ halfar.analysis.Analysis
+ halfar.analysis.Analysis.run
+ halfar.analysis.Analysis.rmse
+
+ horizontal_advection.HorizontalAdvection
+ horizontal_advection.HorizontalAdvection.create_init
+ horizontal_advection.HorizontalAdvection.create_analysis
+
+ horizontal_advection.init.Init
+ horizontal_advection.init.Init.run
+
+ horizontal_advection.analysis.Analysis
+ horizontal_advection.analysis.Analysis.run
+ horizontal_advection.analysis.Analysis.rmse
+
+ horizontal_advection_thickness.HorizontalAdvectionThickness
+ horizontal_advection_thickness.HorizontalAdvectionThickness.create_init
+ horizontal_advection_thickness.HorizontalAdvectionThickness.create_analysis
+
+ horizontal_advection_thickness.init.Init
+ horizontal_advection_thickness.init.Init.run
+
+ horizontal_advection_thickness.analysis.Analysis
+ horizontal_advection_thickness.analysis.Analysis.run
+ horizontal_advection_thickness.analysis.Analysis.rmse
+
mesh_modifications
~~~~~~~~~~~~~~~~~~
diff --git a/docs/developers_guide/landice/test_groups/index.rst b/docs/developers_guide/landice/test_groups/index.rst
index a5ee4cc659..6d9db2558b 100644
--- a/docs/developers_guide/landice/test_groups/index.rst
+++ b/docs/developers_guide/landice/test_groups/index.rst
@@ -23,6 +23,7 @@ Test groups
isunnguata_sermia
kangerlussuaq
koge_bugt_s
+ mesh_convergence
mesh_modifications
mismipplus
slm
diff --git a/docs/developers_guide/landice/test_groups/mesh_convergence.rst b/docs/developers_guide/landice/test_groups/mesh_convergence.rst
new file mode 100644
index 0000000000..bb00d09927
--- /dev/null
+++ b/docs/developers_guide/landice/test_groups/mesh_convergence.rst
@@ -0,0 +1,145 @@
+.. _dev_landice_mesh_convergence:
+
+mesh_convergence
+================
+
+The ``mesh_convergence`` test group
+(:py:class:`compass.landice.tests.mesh_convergence.MeshConvergence`)
+implements a suite of spatial convergence tests for MALI on planar, doubly
+periodic hexagonal meshes (see :ref:`landice_mesh_convergence`). Here we
+describe the shared framework and the individual test cases.
+
+.. _dev_landice_mesh_convergence_framework:
+
+framework
+---------
+
+The test group shares several base classes and config/namelist/streams files
+that are used by all three test cases.
+
+ConvTestCase
+~~~~~~~~~~~~
+
+:py:class:`compass.landice.tests.mesh_convergence.conv_test_case.ConvTestCase`
+is the shared base class for all convergence test cases. It reads the list of
+resolutions from ``[mesh_convergence] resolutions`` and calls
+:py:meth:`~compass.landice.tests.mesh_convergence.conv_test_case.ConvTestCase._setup_steps`
+to construct the appropriate ``init`` and ``forward`` steps for each
+resolution. Each child class must implement ``create_init`` and
+``create_analysis`` to supply the test-case-specific init and analysis step
+objects.
+
+The ``configure`` method re-runs ``_setup_steps`` in case the user has changed
+the resolution list in a config file, and then calls ``update_cores`` to set
+the number of MPI tasks for each forward step based on the mesh size and the
+``goal_cells_per_core`` / ``max_cells_per_core`` config options.
+
+ConvInit
+~~~~~~~~
+
+:py:class:`compass.landice.tests.mesh_convergence.conv_init.ConvInit`
+is the shared base class for all init steps. Its ``run`` method creates a
+planar hexagonal mesh of the requested resolution using
+``mpas_tools.planar_hex.make_planar_hex_mesh``, culls and converts the mesh,
+centers it, and writes ``mesh.nc`` and ``graph.info``. Each test case
+provides a child class of ``ConvInit`` that calls ``super().run()`` and then
+adds the appropriate initial conditions for that test.
+
+ConvAnalysis
+~~~~~~~~~~~~
+
+:py:class:`compass.landice.tests.mesh_convergence.conv_analysis.ConvAnalysis`
+is the shared base class for all analysis steps. Its ``__init__`` method
+registers the forward-step output files as analysis inputs using the naming
+convention ``{resolution}km_output.nc → ../{resolution}km/forward/output.nc``.
+Each test case provides a child class that implements the ``run`` method and,
+typically, an ``rmse`` helper method.
+
+Forward
+~~~~~~~
+
+:py:class:`compass.landice.tests.mesh_convergence.forward.Forward`
+is shared by all test cases. It reads the common ``namelist.forward`` and
+``streams.forward`` files from the ``mesh_convergence`` package, and
+additionally reads test-case-specific namelist/streams files if they exist in
+the child test-case package. The time step for each resolution is set via the
+``get_dt_duration`` method to satisfy an advective CFL condition based on
+``target_velocity`` and scaled proportionally to each mesh resolution.
+
+Config / namelist / streams files
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+``mesh_convergence.cfg`` sets default values for options shared across all
+test cases (resolutions, domain size, core counts, etc.). Each test case
+additionally ships its own ``.cfg`` with test-specific defaults.
+
+``namelist.forward`` sets common MALI namelist options (e.g. third-order
+spatial and temporal advection schemes) and is applied first; test-case
+namelist files are applied on top to allow overrides.
+
+``streams.forward`` and ``streams.template`` define the ``input``
+and ``output`` streams; the output interval is filled in by the ``Forward``
+step from the computed run duration.
+
+.. _dev_landice_mesh_convergence_test_cases:
+
+test cases
+----------
+
+horizontal_advection
+~~~~~~~~~~~~~~~~~~~~
+
+:py:class:`compass.landice.tests.mesh_convergence.horizontal_advection.HorizontalAdvection`
+tests convergence of passive-tracer advection.
+
+Init
+(:py:class:`compass.landice.tests.mesh_convergence.horizontal_advection.init.Init`)
+sets up a uniform ice slab with a Gaussian blob of ``passiveTracer2d`` at the
+domain center plus a spatially uniform prescribed velocity
+(``uReconstructX``, ``uReconstructY``) chosen so that the tracer completes
+exactly one circuit of the periodic domain in the configured run duration.
+
+Analysis
+(:py:class:`compass.landice.tests.mesh_convergence.horizontal_advection.analysis.Analysis`)
+computes the RMSE of ``passiveTracer2d`` between the initial and final states
+at each resolution, fits a power law, and saves ``convergence.png``.
+
+horizontal_advection_thickness
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+:py:class:`compass.landice.tests.mesh_convergence.horizontal_advection_thickness.HorizontalAdvectionThickness`
+is structurally identical to ``horizontal_advection`` but tests convergence of
+ice-thickness advection instead.
+
+Init
+(:py:class:`compass.landice.tests.mesh_convergence.horizontal_advection_thickness.init.Init`)
+places a Gaussian bump of ice thickness on top of a uniform background layer,
+with the same prescribed velocity approach as the tracer test.
+
+Analysis
+(:py:class:`compass.landice.tests.mesh_convergence.horizontal_advection_thickness.analysis.Analysis`)
+computes the RMSE of ``thickness`` between the initial and final states.
+
+halfar
+~~~~~~
+
+:py:class:`compass.landice.tests.mesh_convergence.halfar.Halfar`
+tests convergence of the ice-dynamics solver by comparing against the Halfar
+analytic similarity solution.
+
+Init
+(:py:class:`compass.landice.tests.mesh_convergence.halfar.init.Init`)
+calls
+:py:func:`compass.landice.tests.dome.setup_mesh.setup_dome_initial_conditions`
+to create the Halfar dome initial condition on the planar hexagonal mesh.
+
+Analysis
+(:py:class:`compass.landice.tests.mesh_convergence.halfar.analysis.Analysis`)
+computes two error metrics:
+
+* The RMSE between simulated and analytic Halfar thickness over all cells
+ that contain ice in either solution (``convergence_rmse.png``).
+* The absolute error at the dome center (``convergence_dome.png``).
+
+The order of convergence is derived from the RMSE metric and compared against
+``conv_thresh`` (fail below) and ``conv_max`` (warn above) config options.
diff --git a/docs/users_guide/landice/test_groups/index.rst b/docs/users_guide/landice/test_groups/index.rst
index 438d2dbb09..1e514ef951 100644
--- a/docs/users_guide/landice/test_groups/index.rst
+++ b/docs/users_guide/landice/test_groups/index.rst
@@ -28,6 +28,7 @@ physics but that are not run routinely.
isunnguata_sermia
kangerlussuaq
koge_bugt_s
+ mesh_convergence
mesh_modifications
mismipplus
slm
diff --git a/docs/users_guide/landice/test_groups/mesh_convergence.rst b/docs/users_guide/landice/test_groups/mesh_convergence.rst
new file mode 100644
index 0000000000..209d95d926
--- /dev/null
+++ b/docs/users_guide/landice/test_groups/mesh_convergence.rst
@@ -0,0 +1,188 @@
+.. _landice_mesh_convergence:
+
+mesh_convergence
+================
+
+The ``landice/mesh_convergence`` test group includes tests for assessing the
+spatial (mesh) convergence of MALI on planar, doubly periodic hexagonal
+meshes. Currently, three test cases are available:
+:ref:`landice_mesh_convergence_horizontal_advection`,
+:ref:`landice_mesh_convergence_horizontal_advection_thickness`, and
+:ref:`landice_mesh_convergence_halfar`.
+
+Each test case runs MALI at a series of resolutions, then computes and plots
+the order of convergence of an error metric as a function of the number of
+cells in the mesh. The analysis step raises an error if the computed order of
+convergence falls below a configurable threshold.
+
+config options
+--------------
+
+All test cases in this test group share the following ``[mesh_convergence]``
+config options, which may be overridden on a per-test-case basis:
+
+.. code-block:: cfg
+
+ # options for mesh convergence test cases
+ [mesh_convergence]
+
+ # a list of resolutions (km) to test
+ resolutions = 8, 4, 2, 1
+
+ # number of mesh cells in x and y for 1 km resolution. Other resolutions
+ # have the same physical size.
+ nx_1km = 512
+ ny_1km = 640
+
+ # whether to make mesh nonperiodic or not
+ nonperiodic = False
+
+ # the number of cells per core to aim for
+ goal_cells_per_core = 300
+
+ # the approximate maximum number of cells per core
+ max_cells_per_core = 3000
+
+ # target velocity (m/yr) used to define the time step via an advective CFL
+ # condition
+ target_velocity = 2000.0
+
+ # the duration (years) of the run
+ duration = 1000
+
+.. _landice_mesh_convergence_horizontal_advection:
+
+horizontal_advection
+--------------------
+
+``landice/mesh_convergence/horizontal_advection`` tests the spatial convergence
+of passive tracer advection in MALI. A Gaussian blob of a 2-D passive tracer
+(``passiveTracer2d``) is placed on a doubly periodic mesh and advected with a
+prescribed, spatially uniform velocity exactly once around the domain. The
+root-mean-squared difference between the final and initial tracer fields is
+then used as the error metric.
+
+The test case config options are:
+
+.. code-block:: cfg
+
+ # options for planar horizontal advection test case
+ [horizontal_advection]
+
+ # Number of vertical levels
+ vert_levels = 3
+
+ # ice thickness (m)
+ ice_thickness = 1000.0
+
+ # bed elevation (m)
+ bed_elevation = 0.0
+
+ # center of the tracer gaussian (km)
+ x_center = 0.
+ y_center = 0.
+
+ # width of gaussian tracer "blob" (km)
+ gaussian_width = 50
+
+ # whether to advect in x, y, or both
+ advect_x = True
+ advect_y = True
+
+ # convergence threshold below which the test fails
+ conv_thresh = 0.6
+
+ # Convergence rate above which a warning is issued
+ conv_max = 3.0
+
+ # options for mesh convergence test cases
+ [mesh_convergence]
+
+ # a list of resolutions (km) to test
+ resolutions = 16, 8, 4, 2
+
+.. _landice_mesh_convergence_horizontal_advection_thickness:
+
+horizontal_advection_thickness
+-------------------------------
+
+``landice/mesh_convergence/horizontal_advection_thickness`` tests the spatial
+convergence of ice-thickness advection in MALI. A Gaussian bump of ice
+thickness (on top of a uniform background layer) is placed on a doubly
+periodic mesh and advected with a prescribed, spatially uniform velocity
+exactly once around the domain. The root-mean-squared difference between the
+final and initial thickness fields is used as the error metric.
+
+The test case uses the same config options as
+:ref:`landice_mesh_convergence_horizontal_advection`.
+
+.. _landice_mesh_convergence_halfar:
+
+halfar
+------
+
+``landice/mesh_convergence/halfar`` tests the spatial convergence of MALI's
+ice dynamics by comparing the simulated ice-sheet geometry to the Halfar
+analytic similarity solution (Halfar 1983). The test begins from the
+analytic Halfar ice dome (using the same initial condition as the
+:ref:`landice_dome` test case with ``dome_type = halfar``) and runs the model
+for a configurable duration. Two error metrics are computed at the end of the
+run:
+
+* **RMSE** — the root-mean-squared difference between the simulated and
+ analytic thickness fields over all cells that contain ice in either solution,
+ plotted in ``convergence_rmse.png``.
+* **dome center error** — the absolute difference between the simulated and
+ analytic thickness at the dome center, plotted in ``convergence_dome.png``.
+
+The order of convergence is estimated from the RMSE metric. The test fails if
+the order of convergence is below ``conv_thresh`` and issues a warning if it
+exceeds ``conv_max``.
+
+The test case config options are:
+
+.. code-block:: cfg
+
+ # options for halfar mesh convergence test
+ [halfar]
+
+ # Number of vertical levels
+ vert_levels = 10
+
+ # convergence threshold below which the test fails
+ conv_thresh = 0.0
+
+ # Convergence rate above which a warning is issued
+ conv_max = 3.0
+
+ # config options for dome test cases
+ [dome]
+
+ # the dome type ('halfar' or 'cism')
+ dome_type = halfar
+
+ # Whether to center the dome in the center of the cell that is closest to
+ # the center of the domain
+ put_origin_on_a_cell = True
+
+ # whether to add a small shelf to the test
+ shelf = False
+
+ # whether to add hydrology to the initial condition
+ hydro = False
+
+ [mesh_convergence]
+
+ # a list of resolutions (km) to test
+ resolutions = 8, 4, 2, 1
+
+ nx_1km = 128
+ ny_1km = 128
+ nonperiodic = True
+
+ # target velocity (m/yr); artificially large to satisfy the more
+ # restrictive diffusive CFL condition for Halfar
+ target_velocity = 30000.0
+
+ # the duration (years) of the run
+ duration = 200