diff --git a/compass/landice/__init__.py b/compass/landice/__init__.py
index 97b33a010c..e0e38ee71e 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.thwaites import Thwaites
@@ -47,6 +48,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(Thwaites(mpas_core=self))
diff --git a/compass/landice/tests/dome/setup_mesh.py b/compass/landice/tests/dome/setup_mesh.py
index 912cf3dac8..9436d61c85 100644
--- a/compass/landice/tests/dome/setup_mesh.py
+++ b/compass/landice/tests/dome/setup_mesh.py
@@ -1,10 +1,9 @@
import numpy
-from netCDF4 import Dataset as NetCDFFile
-
-from mpas_tools.planar_hex import make_planar_hex_mesh
from mpas_tools.io import write_netcdf
-from mpas_tools.mesh.conversion import convert, cull
from mpas_tools.logging import check_call
+from mpas_tools.mesh.conversion import convert, cull
+from mpas_tools.planar_hex import make_planar_hex_mesh
+from netCDF4 import Dataset as NetCDFFile
from compass.model import make_graph_file
from compass.step import Step
@@ -81,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
@@ -158,7 +157,7 @@ def _setup_dome_initial_conditions(config, logger, filename):
thickness_field[r < r0] = h0 * (1.0 - (r[r < r0] / r0) ** 2) ** 0.5
elif dome_type == 'halfar':
thickness_field[r < r0] = h0 * (
- 1.0 - (r[r < r0] / r0) ** (4.0 / 3.0)) ** (3.0 / 7.0)
+ 1.0 - (r[r < r0] / r0) ** (4.0 / 3.0)) ** (3.0 / 7.0)
else:
raise ValueError('Unexpected dome_type: {}'.format(dome_type))
thickness[0, :] = thickness_field
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..31b27704af
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/forward.py
@@ -0,0 +1,150 @@
+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
+ assert time_err < 1e-4, 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..fbeaf5aae0
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/halfar/analysis.py
@@ -0,0 +1,189 @@
+import sys
+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(1)
+ 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(2)
+ 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 tolerence {conv_thresh}')
+
+ if conv > conv_max:
+ warnings.warn(f'order of convergence '
+ f'{conv} > max tolerence {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)
+ # 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:
+ sys.exit('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":
+ sys.exit('Error: The Halfar script currently assumes a '
+ 'gregorian_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 = np.absolute(diff[center_idx])
+
+ 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..45f3e63019
--- /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 divisble 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..78761f9b49
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/halfar/namelist.forward
@@ -0,0 +1,3 @@
+config_velocity_solver = 'sia'
+config_thickness_advection = 'fo'
+config_tracer_advection = 'fo'
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..649c3a017d
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py
@@ -0,0 +1,102 @@
+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.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 tolerence {conv_thresh}')
+
+ if conv > conv_max:
+ warnings.warn(f'order of convergence '
+ f'{conv} > max tolerence {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))
+ 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..b19e2aa314
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection/horizontal_advection.cfg
@@ -0,0 +1,28 @@
+# 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
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/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..7bb25a60e8
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py
@@ -0,0 +1,103 @@
+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.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 tolerence {conv_thresh}')
+
+ if conv > conv_max:
+ warnings.warn(f'order of convergence '
+ f'{conv} > max tolerence {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))
+ 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..9fe876ca03
--- /dev/null
+++ b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/namelist.forward
@@ -0,0 +1,5 @@
+config_velocity_solver = 'none'
+config_unrealistic_velocity = 1.0e36
+config_thickness_advection = 'fo'
+config_tracer_advection = 'none'
+config_start_time = '0001-01-01_00:00:00'
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 @@
+
+
+
+
+