From e51e3bf361b8969d55b3d4705b1ceb77efbb332e Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 18 Aug 2023 22:07:01 -0700 Subject: [PATCH 01/26] Add landice mesh_convergence test case This is copied from the ocean's planar_convergence test case with minor adjustments to work for the landice core. --- compass/landice/__init__.py | 2 + .../tests/mesh_convergence/__init__.py | 18 +++ .../tests/mesh_convergence/conv_analysis.py | 34 +++++ .../tests/mesh_convergence/conv_init.py | 62 ++++++++ .../tests/mesh_convergence/conv_test_case.py | 143 ++++++++++++++++++ .../landice/tests/mesh_convergence/forward.py | 136 +++++++++++++++++ .../horizontal_advection/__init__.py | 63 ++++++++ .../horizontal_advection/analysis.py | 99 ++++++++++++ .../horizontal_advection.cfg | 28 ++++ .../horizontal_advection/init.py | 107 +++++++++++++ .../mesh_convergence/mesh_convergence.cfg | 24 +++ .../tests/mesh_convergence/namelist.forward | 3 + .../tests/mesh_convergence/streams.forward | 24 +++ .../tests/mesh_convergence/streams.template | 6 + 14 files changed, 749 insertions(+) create mode 100644 compass/landice/tests/mesh_convergence/__init__.py create mode 100644 compass/landice/tests/mesh_convergence/conv_analysis.py create mode 100644 compass/landice/tests/mesh_convergence/conv_init.py create mode 100644 compass/landice/tests/mesh_convergence/conv_test_case.py create mode 100644 compass/landice/tests/mesh_convergence/forward.py create mode 100644 compass/landice/tests/mesh_convergence/horizontal_advection/__init__.py create mode 100644 compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py create mode 100644 compass/landice/tests/mesh_convergence/horizontal_advection/horizontal_advection.cfg create mode 100644 compass/landice/tests/mesh_convergence/horizontal_advection/init.py create mode 100644 compass/landice/tests/mesh_convergence/mesh_convergence.cfg create mode 100644 compass/landice/tests/mesh_convergence/namelist.forward create mode 100644 compass/landice/tests/mesh_convergence/streams.forward create mode 100644 compass/landice/tests/mesh_convergence/streams.template diff --git a/compass/landice/__init__.py b/compass/landice/__init__.py index 338f5237a2..317520a130 100644 --- a/compass/landice/__init__.py +++ b/compass/landice/__init__.py @@ -15,6 +15,7 @@ from compass.landice.tests.kangerlussuaq import Kangerlussuaq from compass.landice.tests.koge_bugt_s import KogeBugtS from compass.landice.tests.mesh_modifications import MeshModifications +from compass.landice.tests.mesh_convergence import MeshConvergence from compass.landice.tests.mismipplus import MISMIPplus from compass.landice.tests.slm import Slm from compass.landice.tests.thwaites import Thwaites @@ -49,6 +50,7 @@ def __init__(self): self.add_test_group(Kangerlussuaq(mpas_core=self)) self.add_test_group(KogeBugtS(mpas_core=self)) self.add_test_group(MeshModifications(mpas_core=self)) + self.add_test_group(MeshConvergence(mpas_core=self)) self.add_test_group(MISMIPplus(mpas_core=self)) self.add_test_group(Slm(mpas_core=self)) self.add_test_group(Thwaites(mpas_core=self)) diff --git a/compass/landice/tests/mesh_convergence/__init__.py b/compass/landice/tests/mesh_convergence/__init__.py new file mode 100644 index 0000000000..83bfa57b6d --- /dev/null +++ b/compass/landice/tests/mesh_convergence/__init__.py @@ -0,0 +1,18 @@ +from compass.landice.tests.mesh_convergence.horizontal_advection import ( + HorizontalAdvection, +) +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)) 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..aa9978963b --- /dev/null +++ b/compass/landice/tests/mesh_convergence/conv_init.py @@ -0,0 +1,62 @@ +from mpas_tools.io import write_netcdf +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 + """ + 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 + + ds_mesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, + nonperiodic_x=False, + nonperiodic_y=False) + + 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..1c819c7a7f --- /dev/null +++ b/compass/landice/tests/mesh_convergence/forward.py @@ -0,0 +1,136 @@ +import time +from datetime import timedelta +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 base 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 + """ + config = self.config + # dt is proportional to resolution: default 30 seconds per km + dt_1km = config.getint('mesh_convergence', 'dt_1km') + + dt = dt_1km * self.resolution + # https://stackoverflow.com/a/1384565/7728169 + dt = time.strftime('%H:%M:%S', time.gmtime(dt)) + + # the duration (hours) of the run + duration = \ + int(3600 * config.getfloat('mesh_convergence', 'duration')) + delta = timedelta(seconds=duration) + hours = delta.seconds // 3600 + minutes = delta.seconds // 60 % 60 + seconds = delta.seconds % 60 + duration = f'{delta.days:03d}_{hours:02d}:{minutes:02d}:{seconds:02d}' + + 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/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..34d5913132 --- /dev/null +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py @@ -0,0 +1,99 @@ +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='passiveTracer') + 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) + 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..f5d407ec4b --- /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 = 1.9 + +# Convergence rate above which a warning is issued +conv_max = 2.3 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..ff53bcb464 --- /dev/null +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/init.py @@ -0,0 +1,107 @@ +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 = 3600. * section.getfloat('duration') + dt_1km = section.getint('dt_1km') + resolution = float(self.resolution) + dt = dt_1km * resolution + dc = resolution * 1e3 + + ds = xarray.open_dataset('mesh.nc') + xCell = ds.xCell + yCell = ds.yCell + + if advect_x: + x_vel = ds.attrs['x_period'] / duration + x_cfl = x_vel * dt / dc + print(f'x_cfl: {x_cfl}') + else: + x_vel = 0. + + if advect_y: + y_vel = ds.attrs['y_period'] / duration + y_cfl = y_vel * dt / dc + print(f'y_cfl: {y_cfl}') + 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 + + tracer1 = numpy.exp(-0.5 * dist_sq / gaussian_width**2) + tracer1 = tracer1.expand_dims(dim={'nVertLevels': nVertLevels}, + axis=1) + tracer1 = tracer1.expand_dims(dim='Time', axis=0) + + ds['layerThicknessFractions'] = layerThicknessFractions + ds['thickness'] = thickness + ds['bedTopography'] = bedTopography + ds['uReconstructX'] = uReconstructX + ds['uReconstructY'] = uReconstructY + ds['passiveTracer'] = tracer1 + + write_netcdf(ds, 'initial_state.nc') 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..e71280706f --- /dev/null +++ b/compass/landice/tests/mesh_convergence/mesh_convergence.cfg @@ -0,0 +1,24 @@ +# options for mesh convergence test cases +[mesh_convergence] + +# a list of resolutions (km) to test +resolutions = 2, 4, 8, 16, 32 + +# 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. +nx_1km = 512 +ny_1km = 640 + +# 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 + +# time step at 1 km. dt at other resolutions is proportional to the resolution +dt_1km = 15 + +# the duration (hours) of the run +duration = 24 diff --git a/compass/landice/tests/mesh_convergence/namelist.forward b/compass/landice/tests/mesh_convergence/namelist.forward new file mode 100644 index 0000000000..6d917b45b6 --- /dev/null +++ b/compass/landice/tests/mesh_convergence/namelist.forward @@ -0,0 +1,3 @@ +config_velocity_solver = 'none' +config_thickness_advection = 'fo' +config_tracer_advection = 'fo' diff --git a/compass/landice/tests/mesh_convergence/streams.forward b/compass/landice/tests/mesh_convergence/streams.forward new file mode 100644 index 0000000000..5c457e680d --- /dev/null +++ b/compass/landice/tests/mesh_convergence/streams.forward @@ -0,0 +1,24 @@ + + + + + + + + + + + + + + + + + + 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 @@ + + + + + From a812a87d2591ee806014f3b6b5f2c5e19f56108d Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 18 Aug 2023 23:18:56 -0700 Subject: [PATCH 02/26] Correct passiceTracer to have no vertical dim --- .../tests/mesh_convergence/horizontal_advection/init.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/init.py b/compass/landice/tests/mesh_convergence/horizontal_advection/init.py index ff53bcb464..08b711bb41 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection/init.py +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/init.py @@ -92,16 +92,14 @@ def run(self): dist_sq = (xCell - x_center)**2 + (yCell - y_center)**2 - tracer1 = numpy.exp(-0.5 * dist_sq / gaussian_width**2) - tracer1 = tracer1.expand_dims(dim={'nVertLevels': nVertLevels}, - axis=1) - tracer1 = tracer1.expand_dims(dim='Time', axis=0) + 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['passiveTracer'] = tracer1 + ds['passiveTracer'] = passiveTracer write_netcdf(ds, 'initial_state.nc') From fcc862af79b45d62b19e89649f9e81cb301ef935 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Sun, 20 Aug 2023 19:51:35 -0700 Subject: [PATCH 03/26] Remove velo cap to allow prescribed velocity to be used --- compass/landice/tests/mesh_convergence/namelist.forward | 1 + 1 file changed, 1 insertion(+) diff --git a/compass/landice/tests/mesh_convergence/namelist.forward b/compass/landice/tests/mesh_convergence/namelist.forward index 6d917b45b6..326d8c5239 100644 --- a/compass/landice/tests/mesh_convergence/namelist.forward +++ b/compass/landice/tests/mesh_convergence/namelist.forward @@ -1,3 +1,4 @@ config_velocity_solver = 'none' +config_unrealistic_velocity = 1.0e36 config_thickness_advection = 'fo' config_tracer_advection = 'fo' From 28828509dbf5dc78c497d5b9bcc43f4980378878 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 21 Aug 2023 12:42:09 -0700 Subject: [PATCH 04/26] Add Halfar mesh convergence test Also modify the mesh_convergence test-group-level architecture to be flexible enough to support these two different test cases. --- compass/landice/tests/dome/setup_mesh.py | 6 +- .../tests/mesh_convergence/__init__.py | 2 + .../tests/mesh_convergence/conv_init.py | 9 +- .../landice/tests/mesh_convergence/forward.py | 22 +-- .../tests/mesh_convergence/halfar/__init__.py | 60 +++++++ .../tests/mesh_convergence/halfar/analysis.py | 161 ++++++++++++++++++ .../tests/mesh_convergence/halfar/halfar.cfg | 43 +++++ .../tests/mesh_convergence/halfar/init.py | 57 +++++++ .../mesh_convergence/halfar/namelist.forward | 4 + .../mesh_convergence/halfar/streams.forward | 24 +++ 10 files changed, 367 insertions(+), 21 deletions(-) create mode 100644 compass/landice/tests/mesh_convergence/halfar/__init__.py create mode 100644 compass/landice/tests/mesh_convergence/halfar/analysis.py create mode 100644 compass/landice/tests/mesh_convergence/halfar/halfar.cfg create mode 100644 compass/landice/tests/mesh_convergence/halfar/init.py create mode 100644 compass/landice/tests/mesh_convergence/halfar/namelist.forward create mode 100644 compass/landice/tests/mesh_convergence/halfar/streams.forward diff --git a/compass/landice/tests/dome/setup_mesh.py b/compass/landice/tests/dome/setup_mesh.py index 8eb3c2e34d..09d5f340a0 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 diff --git a/compass/landice/tests/mesh_convergence/__init__.py b/compass/landice/tests/mesh_convergence/__init__.py index 83bfa57b6d..5c60a077ea 100644 --- a/compass/landice/tests/mesh_convergence/__init__.py +++ b/compass/landice/tests/mesh_convergence/__init__.py @@ -1,3 +1,4 @@ +from compass.landice.tests.mesh_convergence.halfar import Halfar from compass.landice.tests.mesh_convergence.horizontal_advection import ( HorizontalAdvection, ) @@ -16,3 +17,4 @@ def __init__(self, mpas_core): super().__init__(mpas_core=mpas_core, name='mesh_convergence') self.add_test_case(HorizontalAdvection(test_group=self)) + self.add_test_case(Halfar(test_group=self)) diff --git a/compass/landice/tests/mesh_convergence/conv_init.py b/compass/landice/tests/mesh_convergence/conv_init.py index aa9978963b..4e41c3d3e8 100644 --- a/compass/landice/tests/mesh_convergence/conv_init.py +++ b/compass/landice/tests/mesh_convergence/conv_init.py @@ -1,4 +1,5 @@ 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 @@ -42,6 +43,7 @@ def run(self): """ Run this step of the test case """ + logger = self.logger config = self.config resolution = float(self.resolution) @@ -51,12 +53,15 @@ def run(self): 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=False, - nonperiodic_y=False) + nonperiodic_x=nonperiodic, + nonperiodic_y=nonperiodic) center(ds_mesh) + ds_mesh = cull(ds_mesh, logger=logger) + ds_mesh = convert(ds_mesh, logger=logger) write_netcdf(ds_mesh, 'mesh.nc') make_graph_file('mesh.nc', 'graph.info') diff --git a/compass/landice/tests/mesh_convergence/forward.py b/compass/landice/tests/mesh_convergence/forward.py index 1c819c7a7f..a92c958dd1 100644 --- a/compass/landice/tests/mesh_convergence/forward.py +++ b/compass/landice/tests/mesh_convergence/forward.py @@ -1,5 +1,3 @@ -import time -from datetime import timedelta from importlib.resources import contents from compass.model import run_model @@ -59,7 +57,7 @@ def __init__(self, test_case, resolution): def setup(self): """ - Set namelist options base on config options + Set namelist options based on config options """ namelist_options, stream_replacements = self.get_dt_duration() @@ -104,21 +102,13 @@ def get_dt_duration(self): Namelist options to update """ config = self.config - # dt is proportional to resolution: default 30 seconds per km + # dt is proportional to resolution dt_1km = config.getint('mesh_convergence', 'dt_1km') + dt = dt_1km * self.resolution * 3600.0 * 24.0 - dt = dt_1km * self.resolution - # https://stackoverflow.com/a/1384565/7728169 - dt = time.strftime('%H:%M:%S', time.gmtime(dt)) - - # the duration (hours) of the run - duration = \ - int(3600 * config.getfloat('mesh_convergence', 'duration')) - delta = timedelta(seconds=duration) - hours = delta.seconds // 3600 - minutes = delta.seconds // 60 % 60 - seconds = delta.seconds % 60 - duration = f'{delta.days:03d}_{hours:02d}:{minutes:02d}:{seconds:02d}' + # the duration (days) of the run + duration = config.getint('mesh_convergence', 'duration') + duration = f'{duration:05d}_00:00:00' namelist_replacements = {'config_dt': f"'{dt}'", 'config_run_duration': f"'{duration}'"} 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..02cec73f1d --- /dev/null +++ b/compass/landice/tests/mesh_convergence/halfar/analysis.py @@ -0,0 +1,161 @@ +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.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) + 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) + plt.savefig('convergence.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 + + return rms_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..75b07c777f --- /dev/null +++ b/compass/landice/tests/mesh_convergence/halfar/halfar.cfg @@ -0,0 +1,43 @@ +# 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. The default is approximately square, because of the +# staggering of the hex mesh. +nx_1km = 64 +ny_1km = 64 +nonperiodic = True + +# time step at 1 km in days. dt at other resolutions is proportional to the resolution +dt_1km = 73 + +# the duration (days) of the run. Needs to be a multiple of the ratio +# of the coarsest to finest resolutions +duration = 58400 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..e2564e21b2 --- /dev/null +++ b/compass/landice/tests/mesh_convergence/halfar/namelist.forward @@ -0,0 +1,4 @@ +config_velocity_solver = 'sia' +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/halfar/streams.forward b/compass/landice/tests/mesh_convergence/halfar/streams.forward new file mode 100644 index 0000000000..226250ed14 --- /dev/null +++ b/compass/landice/tests/mesh_convergence/halfar/streams.forward @@ -0,0 +1,24 @@ + + + + + + + + + + + + + + + + + + From 8335a2c97a3a3f6b9b3840cab985e73105fc095b Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 21 Aug 2023 14:38:11 -0700 Subject: [PATCH 05/26] Make each test case in the group define its own output stream contents Adding both a test group stream and a test case stream will result in the union of their contents, so I've moved test-case-specific variables to the test-case streams file. --- .../mesh_convergence/halfar/streams.forward | 19 +------------------ .../horizontal_advection/streams.forward | 8 ++++++++ .../tests/mesh_convergence/streams.forward | 6 +----- 3 files changed, 10 insertions(+), 23 deletions(-) create mode 100644 compass/landice/tests/mesh_convergence/horizontal_advection/streams.forward diff --git a/compass/landice/tests/mesh_convergence/halfar/streams.forward b/compass/landice/tests/mesh_convergence/halfar/streams.forward index 226250ed14..25123a88f7 100644 --- a/compass/landice/tests/mesh_convergence/halfar/streams.forward +++ b/compass/landice/tests/mesh_convergence/halfar/streams.forward @@ -1,24 +1,7 @@ - - - - - - - - - - - + - 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..14724c16c9 --- /dev/null +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/streams.forward @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/compass/landice/tests/mesh_convergence/streams.forward b/compass/landice/tests/mesh_convergence/streams.forward index 5c457e680d..2894cb36e6 100644 --- a/compass/landice/tests/mesh_convergence/streams.forward +++ b/compass/landice/tests/mesh_convergence/streams.forward @@ -14,11 +14,7 @@ - - - - - + From 4e27b35cb2cdf5328109876092a0c69311e071ad Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 22 Aug 2023 15:14:58 -0700 Subject: [PATCH 06/26] Get both convergence tests working together This commit changes a lot of details about the convergence tests to get them both working at the same time after the halfar one broke the horizontal_advection one. Changes include: * replacing dt_1km cfg option with target_velocity that is then used to figure out an appropriate timestep that can be scaled for all resolutions. Note this change may not be desirable in the end, because there are a number of reasons one might want to directly control the dt. e.g. Halfar has a diffusive CFL constraint, which the current logic to set the dt from the target velocity is unable to deal with, and for FCT, one may want to take substantially smaller dt. * change horizontal_advection to use glacier-like velocities and times * change the duration cfg option to use years * minor cfg and plotting updates --- .../tests/mesh_convergence/conv_init.py | 2 +- .../landice/tests/mesh_convergence/forward.py | 36 +++++++++++++++---- .../tests/mesh_convergence/halfar/analysis.py | 3 ++ .../tests/mesh_convergence/halfar/halfar.cfg | 20 ++++++----- .../mesh_convergence/halfar/namelist.forward | 2 +- .../horizontal_advection/analysis.py | 3 ++ .../horizontal_advection.cfg | 4 +-- .../horizontal_advection/init.py | 10 +----- .../mesh_convergence/mesh_convergence.cfg | 16 ++++++--- 9 files changed, 63 insertions(+), 33 deletions(-) diff --git a/compass/landice/tests/mesh_convergence/conv_init.py b/compass/landice/tests/mesh_convergence/conv_init.py index 4e41c3d3e8..38bbf9d761 100644 --- a/compass/landice/tests/mesh_convergence/conv_init.py +++ b/compass/landice/tests/mesh_convergence/conv_init.py @@ -59,9 +59,9 @@ def run(self): nonperiodic_x=nonperiodic, nonperiodic_y=nonperiodic) - center(ds_mesh) 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/forward.py b/compass/landice/tests/mesh_convergence/forward.py index a92c958dd1..31b27704af 100644 --- a/compass/landice/tests/mesh_convergence/forward.py +++ b/compass/landice/tests/mesh_convergence/forward.py @@ -1,3 +1,4 @@ +import math from importlib.resources import contents from compass.model import run_model @@ -101,14 +102,37 @@ def get_dt_duration(self): options : dict Namelist options to update """ - config = self.config - # dt is proportional to resolution - dt_1km = config.getint('mesh_convergence', 'dt_1km') - dt = dt_1km * self.resolution * 3600.0 * 24.0 - # the duration (days) of the run + sec_in_yr = 3600.0 * 24.0 * 365.0 + + config = self.config duration = config.getint('mesh_convergence', 'duration') - duration = f'{duration:05d}_00:00:00' + 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}'"} diff --git a/compass/landice/tests/mesh_convergence/halfar/analysis.py b/compass/landice/tests/mesh_convergence/halfar/analysis.py index 02cec73f1d..1c7f67032e 100644 --- a/compass/landice/tests/mesh_convergence/halfar/analysis.py +++ b/compass/landice/tests/mesh_convergence/halfar/analysis.py @@ -55,6 +55,9 @@ def run(self): 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.png', bbox_inches='tight', pad_inches=0.1) section = self.config['halfar'] diff --git a/compass/landice/tests/mesh_convergence/halfar/halfar.cfg b/compass/landice/tests/mesh_convergence/halfar/halfar.cfg index 75b07c777f..0f256b99ce 100644 --- a/compass/landice/tests/mesh_convergence/halfar/halfar.cfg +++ b/compass/landice/tests/mesh_convergence/halfar/halfar.cfg @@ -29,15 +29,17 @@ hydro = False [mesh_convergence] # 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. -nx_1km = 64 -ny_1km = 64 +# 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 -# time step at 1 km in days. dt at other resolutions is proportional to the resolution -dt_1km = 73 +# 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 = 3000.0 -# the duration (days) of the run. Needs to be a multiple of the ratio -# of the coarsest to finest resolutions -duration = 58400 +# the duration (years) of the run +duration = 200 diff --git a/compass/landice/tests/mesh_convergence/halfar/namelist.forward b/compass/landice/tests/mesh_convergence/halfar/namelist.forward index e2564e21b2..de24182c1a 100644 --- a/compass/landice/tests/mesh_convergence/halfar/namelist.forward +++ b/compass/landice/tests/mesh_convergence/halfar/namelist.forward @@ -1,4 +1,4 @@ config_velocity_solver = 'sia' config_thickness_advection = 'fo' -config_tracer_advection = 'none' +config_tracer_advection = 'fo' config_start_time = '0001-01-01_00:00:00' diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py index 34d5913132..6f1698a7c8 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py @@ -54,6 +54,9 @@ def run(self): 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'] diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/horizontal_advection.cfg b/compass/landice/tests/mesh_convergence/horizontal_advection/horizontal_advection.cfg index f5d407ec4b..b19e2aa314 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection/horizontal_advection.cfg +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/horizontal_advection.cfg @@ -22,7 +22,7 @@ advect_x = True advect_y = True # convergence threshold below which the test fails -conv_thresh = 1.9 +conv_thresh = 0.6 # Convergence rate above which a warning is issued -conv_max = 2.3 +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 index 08b711bb41..42a74f13f2 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection/init.py +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/init.py @@ -45,11 +45,7 @@ def run(self): nVertLevels = section.getint('vert_levels') section = config['mesh_convergence'] - duration = 3600. * section.getfloat('duration') - dt_1km = section.getint('dt_1km') - resolution = float(self.resolution) - dt = dt_1km * resolution - dc = resolution * 1e3 + duration = section.getfloat('duration') * 3600.0 * 24.0 * 365.0 ds = xarray.open_dataset('mesh.nc') xCell = ds.xCell @@ -57,15 +53,11 @@ def run(self): if advect_x: x_vel = ds.attrs['x_period'] / duration - x_cfl = x_vel * dt / dc - print(f'x_cfl: {x_cfl}') else: x_vel = 0. if advect_y: y_vel = ds.attrs['y_period'] / duration - y_cfl = y_vel * dt / dc - print(f'y_cfl: {y_cfl}') else: y_vel = 0. diff --git a/compass/landice/tests/mesh_convergence/mesh_convergence.cfg b/compass/landice/tests/mesh_convergence/mesh_convergence.cfg index e71280706f..16a5a21e05 100644 --- a/compass/landice/tests/mesh_convergence/mesh_convergence.cfg +++ b/compass/landice/tests/mesh_convergence/mesh_convergence.cfg @@ -2,14 +2,18 @@ [mesh_convergence] # a list of resolutions (km) to test -resolutions = 2, 4, 8, 16, 32 +resolutions = 16, 8, 4, 2 # 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 @@ -17,8 +21,10 @@ goal_cells_per_core = 300 # few cores are available) max_cells_per_core = 3000 -# time step at 1 km. dt at other resolutions is proportional to the resolution -dt_1km = 15 +# 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 (hours) of the run -duration = 24 +# the duration (years) of the run +duration = 1000 From 9224dacb9cc88334d1a777b01f4f01e06328f8c6 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 22 Aug 2023 15:45:38 -0700 Subject: [PATCH 07/26] Add dome center error convergence in addition to rmse All advection will reduce to first order at margin, so the dome center error provides a metric of error "away" from the margin. --- .../tests/mesh_convergence/halfar/analysis.py | 42 +++++++++++++++---- 1 file changed, 33 insertions(+), 9 deletions(-) diff --git a/compass/landice/tests/mesh_convergence/halfar/analysis.py b/compass/landice/tests/mesh_convergence/halfar/analysis.py index 1c7f67032e..e41aaf1bf9 100644 --- a/compass/landice/tests/mesh_convergence/halfar/analysis.py +++ b/compass/landice/tests/mesh_convergence/halfar/analysis.py @@ -35,22 +35,26 @@ def run(self): plt.switch_backend('Agg') resolutions = self.resolutions ncells_list = list() - errors = list() + rmse_errors = list() + center_errors = list() for res in resolutions: - rms_error, ncells = self.rmse(res) + rms_error, center_error, ncells = self.rmse(res) ncells_list.append(ncells) - errors.append(rms_error) + rmse_errors.append(rms_error) + center_errors.append(center_error) ncells = np.array(ncells_list) - errors = np.array(errors) + rmse_errors = np.array(rmse_errors) + center_errors = np.array(center_errors) - p = np.polyfit(np.log10(ncells), np.log10(errors), 1) + # 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, errors, 'or') + 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) @@ -58,7 +62,24 @@ def run(self): section = self.config['mesh_convergence'] duration = section.getfloat('duration') plt.title(f'Halfar convergence test, {duration} yrs') - plt.savefig('convergence.png', bbox_inches='tight', pad_inches=0.1) + 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') @@ -127,7 +148,10 @@ def rmse(self, resolution): rms_error = ((diff[mask]**2 * areaCell[mask] / areaCell[mask].sum()).sum())**0.5 - return rms_error, ncells + 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): From ef32615cbf0c37b0c648ef4bd1c7ee13adf0c3f7 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 23 Aug 2023 08:10:00 -0700 Subject: [PATCH 08/26] Minor tweaks to configurations * drop 16 km and add 1 km to ensemble. 16 km appears too coarse to be very useful, especially for halfar, where it makes a dome of only 7 grid cells * adjust halfar target_velocity accordingly (accounting for diffusive cfl) * move start time of year 1 to the test group level so it applies to both tests --- compass/landice/tests/mesh_convergence/halfar/halfar.cfg | 2 +- compass/landice/tests/mesh_convergence/halfar/namelist.forward | 1 - compass/landice/tests/mesh_convergence/mesh_convergence.cfg | 2 +- compass/landice/tests/mesh_convergence/namelist.forward | 1 + 4 files changed, 3 insertions(+), 3 deletions(-) diff --git a/compass/landice/tests/mesh_convergence/halfar/halfar.cfg b/compass/landice/tests/mesh_convergence/halfar/halfar.cfg index 0f256b99ce..45f3e63019 100644 --- a/compass/landice/tests/mesh_convergence/halfar/halfar.cfg +++ b/compass/landice/tests/mesh_convergence/halfar/halfar.cfg @@ -39,7 +39,7 @@ nonperiodic = True # and an advective CFL condition. # Note that Halfar is subject to a more restrictive diffusive CFL so this must be # artificially inflated -target_velocity = 3000.0 +target_velocity = 30000.0 # the duration (years) of the run duration = 200 diff --git a/compass/landice/tests/mesh_convergence/halfar/namelist.forward b/compass/landice/tests/mesh_convergence/halfar/namelist.forward index de24182c1a..78761f9b49 100644 --- a/compass/landice/tests/mesh_convergence/halfar/namelist.forward +++ b/compass/landice/tests/mesh_convergence/halfar/namelist.forward @@ -1,4 +1,3 @@ config_velocity_solver = 'sia' 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/mesh_convergence.cfg b/compass/landice/tests/mesh_convergence/mesh_convergence.cfg index 16a5a21e05..f07310f108 100644 --- a/compass/landice/tests/mesh_convergence/mesh_convergence.cfg +++ b/compass/landice/tests/mesh_convergence/mesh_convergence.cfg @@ -2,7 +2,7 @@ [mesh_convergence] # a list of resolutions (km) to test -resolutions = 16, 8, 4, 2 +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 diff --git a/compass/landice/tests/mesh_convergence/namelist.forward b/compass/landice/tests/mesh_convergence/namelist.forward index 326d8c5239..a7970f6389 100644 --- a/compass/landice/tests/mesh_convergence/namelist.forward +++ b/compass/landice/tests/mesh_convergence/namelist.forward @@ -2,3 +2,4 @@ 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' From 42f1207dd45ad10f0eb0cace6a28110d4d0f8273 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Sat, 26 Aug 2023 09:49:25 -0700 Subject: [PATCH 09/26] Add thickness advection version of horizontal_advection test This is identical to the horizontal_advection test of a passive tracer, but Gaussian bump is applied to thickness field and tracer advection is not considered. This isolates the analysis to a single advection call; if there are issues with the thickness advection, we presumably cannot expect tracer advection to converge properly. It might have been possible to create an argument to toggle between versions, but it was easier (and arguably less confusing) to just make a copy. This uses a Gaussian bump of 1000 m sitting on top of a 1000 m base layer in a doubly periodic mesh. This ensures a margin-less domain, which would avoid the issue of advection degrading to 1st order at margins. --- .../tests/mesh_convergence/__init__.py | 4 + .../__init__.py | 64 +++++++++++ .../analysis.py | 103 ++++++++++++++++++ .../horizontal_advection_thickness.cfg | 34 ++++++ .../horizontal_advection_thickness/init.py | 98 +++++++++++++++++ .../namelist.forward | 5 + .../streams.forward | 7 ++ 7 files changed, 315 insertions(+) create mode 100644 compass/landice/tests/mesh_convergence/horizontal_advection_thickness/__init__.py create mode 100644 compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py create mode 100644 compass/landice/tests/mesh_convergence/horizontal_advection_thickness/horizontal_advection_thickness.cfg create mode 100644 compass/landice/tests/mesh_convergence/horizontal_advection_thickness/init.py create mode 100644 compass/landice/tests/mesh_convergence/horizontal_advection_thickness/namelist.forward create mode 100644 compass/landice/tests/mesh_convergence/horizontal_advection_thickness/streams.forward diff --git a/compass/landice/tests/mesh_convergence/__init__.py b/compass/landice/tests/mesh_convergence/__init__.py index 5c60a077ea..86a34034e8 100644 --- a/compass/landice/tests/mesh_convergence/__init__.py +++ b/compass/landice/tests/mesh_convergence/__init__.py @@ -2,6 +2,9 @@ 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 @@ -17,4 +20,5 @@ def __init__(self, mpas_core): 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/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 @@ + + + + + + + From 00591007b31bb8433ae616642adb53ba56a2eba0 Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Tue, 7 Nov 2023 11:46:46 -0800 Subject: [PATCH 10/26] Fix name of passive tracer variable --- .../landice/tests/mesh_convergence/horizontal_advection/init.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/init.py b/compass/landice/tests/mesh_convergence/horizontal_advection/init.py index 42a74f13f2..7f335a15e5 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection/init.py +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/init.py @@ -92,6 +92,6 @@ def run(self): ds['bedTopography'] = bedTopography ds['uReconstructX'] = uReconstructX ds['uReconstructY'] = uReconstructY - ds['passiveTracer'] = passiveTracer + ds['passiveTracer2d'] = passiveTracer write_netcdf(ds, 'initial_state.nc') From 9e6536f47776bb4ed5a10258dfe4bc601a040815 Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Tue, 7 Nov 2023 12:11:11 -0800 Subject: [PATCH 11/26] Update analysis output files have compass expect seperate files for the dome and rmse convergence. --- compass/landice/tests/mesh_convergence/halfar/analysis.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/compass/landice/tests/mesh_convergence/halfar/analysis.py b/compass/landice/tests/mesh_convergence/halfar/analysis.py index e41aaf1bf9..fbeaf5aae0 100644 --- a/compass/landice/tests/mesh_convergence/halfar/analysis.py +++ b/compass/landice/tests/mesh_convergence/halfar/analysis.py @@ -26,7 +26,8 @@ def __init__(self, test_case, resolutions): """ super().__init__(test_case=test_case, resolutions=resolutions) self.resolutions = resolutions - self.add_output_file('convergence.png') + self.add_output_file('convergence_rmse.png') + self.add_output_file('convergence_dome.png') def run(self): """ From 10b9f0661114ab1ae501ae98c5949a956dbc587f Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Tue, 7 Nov 2023 13:09:22 -0800 Subject: [PATCH 12/26] Upate tracer field in the streams file --- .../tests/mesh_convergence/horizontal_advection/streams.forward | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/streams.forward b/compass/landice/tests/mesh_convergence/horizontal_advection/streams.forward index 14724c16c9..9f5d29f952 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection/streams.forward +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/streams.forward @@ -3,6 +3,6 @@ - + From 1a70fa139cb737c5b17aa9009af0c220158a10c6 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 8 Nov 2023 20:34:44 -0700 Subject: [PATCH 13/26] Change restart output interval from 30 days to 30 years. Change restart output interval from 30 days to 30 years. --- compass/landice/tests/mesh_convergence/streams.forward | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/landice/tests/mesh_convergence/streams.forward b/compass/landice/tests/mesh_convergence/streams.forward index 2894cb36e6..87c53d6818 100644 --- a/compass/landice/tests/mesh_convergence/streams.forward +++ b/compass/landice/tests/mesh_convergence/streams.forward @@ -4,7 +4,7 @@ filename_template="init.nc"/> + output_interval="0030-00-00_00:00:00"/> Date: Thu, 9 Nov 2023 07:59:14 -0700 Subject: [PATCH 14/26] Update passiveTracer to passiveTracer2d in analysis.py Update passiveTracer to passiveTracer2d in analysis.py for horizontal_advection. --- .../tests/mesh_convergence/horizontal_advection/analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py index 6f1698a7c8..649c3a017d 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py @@ -36,7 +36,7 @@ def run(self): ncells_list = list() errors = list() for res in resolutions: - rms_error, ncells = self.rmse(res, variable='passiveTracer') + rms_error, ncells = self.rmse(res, variable='passiveTracer2d') ncells_list.append(ncells) errors.append(rms_error) From fdb701635e3c8932feee33659a9c9e41dd6a59e0 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 19 Aug 2024 13:53:43 -0700 Subject: [PATCH 15/26] Change clobber mode from truncate to overwrite Change clobber mode from truncate to overwrite because we want to be able to use restarts without losing the initial time level in the output file. --- compass/landice/tests/mesh_convergence/streams.forward | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/landice/tests/mesh_convergence/streams.forward b/compass/landice/tests/mesh_convergence/streams.forward index 87c53d6818..78bf23f620 100644 --- a/compass/landice/tests/mesh_convergence/streams.forward +++ b/compass/landice/tests/mesh_convergence/streams.forward @@ -9,7 +9,7 @@ From f9c2f07895d3db866e5664292836c0eb4ebf9be3 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 19 Aug 2024 13:57:26 -0700 Subject: [PATCH 16/26] Fix isort issue from linter --- compass/landice/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/compass/landice/__init__.py b/compass/landice/__init__.py index 317520a130..e22cacec79 100644 --- a/compass/landice/__init__.py +++ b/compass/landice/__init__.py @@ -14,8 +14,8 @@ 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_modifications import MeshModifications 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 from compass.landice.tests.thwaites import Thwaites @@ -49,8 +49,8 @@ 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(MeshModifications(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)) self.add_test_group(Thwaites(mpas_core=self)) From 96e08070402ef5a804c81d2565a3ce4e8ac4dfea Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 8 Apr 2026 15:33:20 -0700 Subject: [PATCH 17/26] Suppress warnings with decode_timedelta=False Suppress warnings with decode_timedelta=False when using xr.open_dataset() --- compass/landice/tests/mesh_convergence/halfar/analysis.py | 6 ++++-- .../tests/mesh_convergence/horizontal_advection/analysis.py | 3 ++- .../horizontal_advection_thickness/analysis.py | 3 ++- 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/compass/landice/tests/mesh_convergence/halfar/analysis.py b/compass/landice/tests/mesh_convergence/halfar/analysis.py index fbeaf5aae0..d29d90a6c4 100644 --- a/compass/landice/tests/mesh_convergence/halfar/analysis.py +++ b/compass/landice/tests/mesh_convergence/halfar/analysis.py @@ -115,7 +115,9 @@ def rmse(self, resolution): timelev = -1 # todo: determine a specified time level and err check - ds = xr.open_dataset('{}_output.nc'.format(res_tag), decode_cf=False) + 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) @@ -130,7 +132,7 @@ def rmse(self, resolution): 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') + 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.') diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py index 649c3a017d..ebd8668ce7 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py @@ -93,7 +93,8 @@ def rmse(self, resolution, variable): """ res_tag = '{}km'.format(resolution) - ds = xr.open_dataset('{}_output.nc'.format(res_tag)) + 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 diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py index 7bb25a60e8..357df7f394 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py +++ b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py @@ -94,7 +94,8 @@ def rmse(self, resolution, variable): """ res_tag = '{}km'.format(resolution) - ds = xr.open_dataset('{}_output.nc'.format(res_tag)) + 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 From fb2de3914e63d5dd9891c8164e84a923d9d82f65 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 8 Apr 2026 15:34:12 -0700 Subject: [PATCH 18/26] Make 2, 4, 8, 16km the default resolution for horizontal_advection Make 2, 4, 8, 16km the default resolution for horizontal_advection. This gives the proper order of convergence (2.746 for 3rd-order) and is much faster than going to 1km. --- .../horizontal_advection/horizontal_advection.cfg | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/horizontal_advection.cfg b/compass/landice/tests/mesh_convergence/horizontal_advection/horizontal_advection.cfg index b19e2aa314..22bc39cbe4 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection/horizontal_advection.cfg +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/horizontal_advection.cfg @@ -26,3 +26,9 @@ 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 From 8dfb101af804888670f692804df8e545fc09cfed Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 8 Apr 2026 15:42:07 -0700 Subject: [PATCH 19/26] Make third-order space and time default for convergence tests Make third-order space and time default for convergence tests, since these aren't very useful for first-order. --- .../tests/mesh_convergence/halfar/namelist.forward | 10 ++++++++-- .../horizontal_advection/namelist.forward | 10 ++++++++++ .../horizontal_advection_thickness/namelist.forward | 7 ++++++- 3 files changed, 24 insertions(+), 3 deletions(-) create mode 100644 compass/landice/tests/mesh_convergence/horizontal_advection/namelist.forward diff --git a/compass/landice/tests/mesh_convergence/halfar/namelist.forward b/compass/landice/tests/mesh_convergence/halfar/namelist.forward index 78761f9b49..a4e25ac67f 100644 --- a/compass/landice/tests/mesh_convergence/halfar/namelist.forward +++ b/compass/landice/tests/mesh_convergence/halfar/namelist.forward @@ -1,3 +1,9 @@ config_velocity_solver = 'sia' -config_thickness_advection = 'fo' -config_tracer_advection = 'fo' +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/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_thickness/namelist.forward b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/namelist.forward index 9fe876ca03..82767ef85a 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/namelist.forward +++ b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/namelist.forward @@ -1,5 +1,10 @@ config_velocity_solver = 'none' config_unrealistic_velocity = 1.0e36 -config_thickness_advection = 'fo' +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 From c814164c0075e3ce7777f81aee32fc72fb0a4b3a Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 9 Apr 2026 07:53:35 -0700 Subject: [PATCH 20/26] Create documentation for landice/mesh_convergence cases --- docs/developers_guide/landice/api.rst | 58 ++++++ .../landice/test_groups/index.rst | 1 + .../landice/test_groups/mesh_convergence.rst | 145 ++++++++++++++ .../users_guide/landice/test_groups/index.rst | 1 + .../landice/test_groups/mesh_convergence.rst | 188 ++++++++++++++++++ 5 files changed, 393 insertions(+) create mode 100644 docs/developers_guide/landice/test_groups/mesh_convergence.rst create mode 100644 docs/users_guide/landice/test_groups/mesh_convergence.rst 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..5e2294af13 --- /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 :cite:`Halfar1983`. 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 From 688e4bb7181459fe0f80406ccf0cc93fbe85825e Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 9 Apr 2026 08:10:14 -0700 Subject: [PATCH 21/26] Update MALI-Dev submodule to include recent bug fixes to RK Update MALI-Dev submodule to include recent bug fixes to RK time integration. --- MALI-Dev | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 10006b91c9be688492e1e2923de8c031ef994faf Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 9 Apr 2026 08:59:30 -0700 Subject: [PATCH 22/26] Fix error when building docs Remove use of `:cite:`, which is not available in compass docs. --- docs/users_guide/landice/test_groups/mesh_convergence.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/users_guide/landice/test_groups/mesh_convergence.rst b/docs/users_guide/landice/test_groups/mesh_convergence.rst index 5e2294af13..209d95d926 100644 --- a/docs/users_guide/landice/test_groups/mesh_convergence.rst +++ b/docs/users_guide/landice/test_groups/mesh_convergence.rst @@ -123,7 +123,7 @@ 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 :cite:`Halfar1983`. The test begins from the +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 From c3476831d44e0411ee88952f1e887fa73295e414 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 9 Apr 2026 11:59:14 -0600 Subject: [PATCH 23/26] Apply suggestions from code review Fix various typos. Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- compass/landice/tests/mesh_convergence/forward.py | 4 +++- compass/landice/tests/mesh_convergence/halfar/analysis.py | 6 +++--- compass/landice/tests/mesh_convergence/halfar/halfar.cfg | 2 +- .../tests/mesh_convergence/horizontal_advection/analysis.py | 4 ++-- .../horizontal_advection_thickness/analysis.py | 6 +++--- 5 files changed, 12 insertions(+), 10 deletions(-) diff --git a/compass/landice/tests/mesh_convergence/forward.py b/compass/landice/tests/mesh_convergence/forward.py index 31b27704af..5f83e6ba7c 100644 --- a/compass/landice/tests/mesh_convergence/forward.py +++ b/compass/landice/tests/mesh_convergence/forward.py @@ -129,7 +129,9 @@ def get_dt_duration(self): # 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}' + 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' diff --git a/compass/landice/tests/mesh_convergence/halfar/analysis.py b/compass/landice/tests/mesh_convergence/halfar/analysis.py index d29d90a6c4..955aa86df8 100644 --- a/compass/landice/tests/mesh_convergence/halfar/analysis.py +++ b/compass/landice/tests/mesh_convergence/halfar/analysis.py @@ -88,11 +88,11 @@ def run(self): if conv < conv_thresh: raise ValueError(f'order of convergence ' - f' {conv} < min tolerence {conv_thresh}') + f'{conv} < min tolerance {conv_thresh}') if conv > conv_max: warnings.warn(f'order of convergence ' - f'{conv} > max tolerence {conv_max}') + f'{conv} > max tolerance {conv_max}') def rmse(self, resolution): """ @@ -152,7 +152,7 @@ def rmse(self, resolution): areaCell[mask].sum()).sum())**0.5 center_idx = np.where((xCell == 0.0) * (yCell == 0.0)) - center_error = np.absolute(diff[center_idx]) + center_error = float(np.asarray(np.absolute(diff[center_idx])).item()) return rms_error, center_error, ncells diff --git a/compass/landice/tests/mesh_convergence/halfar/halfar.cfg b/compass/landice/tests/mesh_convergence/halfar/halfar.cfg index 45f3e63019..05f02b253e 100644 --- a/compass/landice/tests/mesh_convergence/halfar/halfar.cfg +++ b/compass/landice/tests/mesh_convergence/halfar/halfar.cfg @@ -30,7 +30,7 @@ hydro = False # 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 +# must be divisible by the ratio of the other meshes resolutions to 1 km nx_1km = 128 ny_1km = 128 nonperiodic = True diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py index ebd8668ce7..51934a08f9 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py @@ -65,11 +65,11 @@ def run(self): if conv < conv_thresh: raise ValueError(f'order of convergence ' - f' {conv} < min tolerence {conv_thresh}') + f'{conv} < min tolerance {conv_thresh}') if conv > conv_max: warnings.warn(f'order of convergence ' - f'{conv} > max tolerence {conv_max}') + f'{conv} > max tolerance {conv_max}') def rmse(self, resolution, variable): """ diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py index 357df7f394..b7d0d5462e 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py +++ b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py @@ -56,7 +56,7 @@ def run(self): plt.ylabel('L2 Norm', fontsize=14) section = self.config['mesh_convergence'] duration = section.getfloat('duration') - plt.title(f'Thickness horizontal advection convergence test,' + plt.title(f'Thickness horizontal advection convergence test, ' f'{duration} yrs') plt.savefig('convergence.png', bbox_inches='tight', pad_inches=0.1) @@ -66,11 +66,11 @@ def run(self): if conv < conv_thresh: raise ValueError(f'order of convergence ' - f' {conv} < min tolerence {conv_thresh}') + f'{conv} < min tolerance {conv_thresh}') if conv > conv_max: warnings.warn(f'order of convergence ' - f'{conv} > max tolerence {conv_max}') + f'{conv} > max tolerance {conv_max}') def rmse(self, resolution, variable): """ From a8ad8461c3a921a6c9e5445af1aa366f0dccb03f Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 9 Apr 2026 11:04:53 -0700 Subject: [PATCH 24/26] Fix typo in landice dome setup Change "hyrdo" typo to "hydro" --- compass/landice/tests/dome/setup_mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/landice/tests/dome/setup_mesh.py b/compass/landice/tests/dome/setup_mesh.py index 09d5f340a0..b5937891a7 100644 --- a/compass/landice/tests/dome/setup_mesh.py +++ b/compass/landice/tests/dome/setup_mesh.py @@ -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+') From 09c2e82623dad70ab8404e63d0b0ad14ae88ffea Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 9 Apr 2026 11:13:13 -0700 Subject: [PATCH 25/26] Raise ValueError instead of sys.exit Raise ValueError instead of sys.exit in halfar analysis for clean error reporting. --- .../landice/tests/mesh_convergence/halfar/analysis.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/compass/landice/tests/mesh_convergence/halfar/analysis.py b/compass/landice/tests/mesh_convergence/halfar/analysis.py index 955aa86df8..e0971ab527 100644 --- a/compass/landice/tests/mesh_convergence/halfar/analysis.py +++ b/compass/landice/tests/mesh_convergence/halfar/analysis.py @@ -1,4 +1,3 @@ -import sys import warnings import matplotlib.pyplot as plt @@ -125,8 +124,8 @@ def rmse(self, resolution): 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.') + 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) @@ -134,8 +133,8 @@ def rmse(self, resolution): 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.') + raise ValueError('Error: The Halfar script currently assumes a ' + 'noleap calendar.') ncells = ds.sizes['nCells'] thk = ds['thickness'].isel(Time=timelev) From 46607bc9048636a54bba224cd340655b4d0b2b79 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 9 Apr 2026 14:16:28 -0700 Subject: [PATCH 26/26] Fix bug in landice mesh convergence plots Fix bug in landice mesh convergence plots that led to all results being plotted on the same axes when multiple test cases are run in the same job. --- compass/landice/tests/mesh_convergence/halfar/analysis.py | 4 ++-- .../tests/mesh_convergence/horizontal_advection/analysis.py | 1 + .../horizontal_advection_thickness/analysis.py | 1 + 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/compass/landice/tests/mesh_convergence/halfar/analysis.py b/compass/landice/tests/mesh_convergence/halfar/analysis.py index e0971ab527..5c7ac81c49 100644 --- a/compass/landice/tests/mesh_convergence/halfar/analysis.py +++ b/compass/landice/tests/mesh_convergence/halfar/analysis.py @@ -52,7 +52,7 @@ def run(self): conv = abs(p[0]) * 2.0 error_fit = ncells**p[0] * 10**p[1] - plt.figure(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)), @@ -66,7 +66,7 @@ def run(self): pad_inches=0.1) # now repeat for center errors - plt.figure(2) + 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] diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py index 51934a08f9..377cbcef7e 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py +++ b/compass/landice/tests/mesh_convergence/horizontal_advection/analysis.py @@ -48,6 +48,7 @@ def run(self): 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)), diff --git a/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py index b7d0d5462e..44de4b0e2e 100644 --- a/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py +++ b/compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py @@ -48,6 +48,7 @@ def run(self): 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)),