From a19398d992926e12edf3993daca804e2a5aaa4f5 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 18 Aug 2023 22:07:01 -0700 Subject: [PATCH 01/16] 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 97b33a010c..58c097b462 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.thwaites import Thwaites from compass.mpas_core import MpasCore @@ -48,5 +49,6 @@ 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(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 6e328bf03299f61cb855657ce117cf21401774c8 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Fri, 18 Aug 2023 23:18:56 -0700 Subject: [PATCH 02/16] 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 df30d40959c0bc43518b0a92779ea1ceb250cb68 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Sun, 20 Aug 2023 19:51:35 -0700 Subject: [PATCH 03/16] 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 60a4e867f5a09528588bf9d7d0af4c0c35d0750a Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 21 Aug 2023 12:42:09 -0700 Subject: [PATCH 04/16] 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 | 15 +- .../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, 371 insertions(+), 26 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 912cf3dac8..9436d61c85 100644 --- a/compass/landice/tests/dome/setup_mesh.py +++ b/compass/landice/tests/dome/setup_mesh.py @@ -1,10 +1,9 @@ import numpy -from netCDF4 import Dataset as NetCDFFile - -from mpas_tools.planar_hex import make_planar_hex_mesh from mpas_tools.io import write_netcdf -from mpas_tools.mesh.conversion import convert, cull from mpas_tools.logging import check_call +from mpas_tools.mesh.conversion import convert, cull +from mpas_tools.planar_hex import make_planar_hex_mesh +from netCDF4 import Dataset as NetCDFFile from compass.model import make_graph_file from compass.step import Step @@ -81,11 +80,11 @@ def run(self): make_graph_file(mesh_filename='landice_grid.nc', graph_filename='graph.info') - _setup_dome_initial_conditions(config, logger, - filename='landice_grid.nc') + setup_dome_initial_conditions(config, logger, + filename='landice_grid.nc') -def _setup_dome_initial_conditions(config, logger, filename): +def setup_dome_initial_conditions(config, logger, filename): """ Add the initial condition to the given MPAS mesh file @@ -158,7 +157,7 @@ def _setup_dome_initial_conditions(config, logger, filename): thickness_field[r < r0] = h0 * (1.0 - (r[r < r0] / r0) ** 2) ** 0.5 elif dome_type == 'halfar': thickness_field[r < r0] = h0 * ( - 1.0 - (r[r < r0] / r0) ** (4.0 / 3.0)) ** (3.0 / 7.0) + 1.0 - (r[r < r0] / r0) ** (4.0 / 3.0)) ** (3.0 / 7.0) else: raise ValueError('Unexpected dome_type: {}'.format(dome_type)) thickness[0, :] = thickness_field diff --git a/compass/landice/tests/mesh_convergence/__init__.py b/compass/landice/tests/mesh_convergence/__init__.py 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 3978d2f03fb55c4e186e10d2fe80848b9439a650 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 21 Aug 2023 14:38:11 -0700 Subject: [PATCH 05/16] 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 52c2b6e61122c8ec8ba83f0fe97b9252dcce0cab Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 22 Aug 2023 15:14:58 -0700 Subject: [PATCH 06/16] 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 cb81780ce2594a830bf44180757705a8ac5658bd Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 22 Aug 2023 15:45:38 -0700 Subject: [PATCH 07/16] 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 8e01f1b707dc61aa5ae6bc273905684166473392 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 23 Aug 2023 08:10:00 -0700 Subject: [PATCH 08/16] 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 ddad348409f9d0623382cbd13f6c1b69e5ae5808 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Sat, 26 Aug 2023 09:49:25 -0700 Subject: [PATCH 09/16] 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 09fea5f4fbcb9d16a68d0ab5a21129c770587089 Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Tue, 7 Nov 2023 11:46:46 -0800 Subject: [PATCH 10/16] 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 71f5ef47d40d576fdeb0d08db4b454c7c29595ba Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Tue, 7 Nov 2023 12:11:11 -0800 Subject: [PATCH 11/16] 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 95004efabf96389a4d9441dc2bca0d26bae35d09 Mon Sep 17 00:00:00 2001 From: Andrew Nolan Date: Tue, 7 Nov 2023 13:09:22 -0800 Subject: [PATCH 12/16] 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 5e79cc005dda4722ed6273203e1fb83e3e06678d Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 8 Nov 2023 20:34:44 -0700 Subject: [PATCH 13/16] 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/16] 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 b017cf3621258fa0bd1bad6ed261b27bcec6e02d Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 19 Aug 2024 13:53:43 -0700 Subject: [PATCH 15/16] 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 a9ebafc2676306ac4d39d9d4e3f77e061deab769 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 19 Aug 2024 13:57:26 -0700 Subject: [PATCH 16/16] 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 58c097b462..e0e38ee71e 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.thwaites import Thwaites from compass.mpas_core import MpasCore @@ -48,7 +48,7 @@ def __init__(self): self.add_test_group(IsunnguataSermia(mpas_core=self)) self.add_test_group(Kangerlussuaq(mpas_core=self)) self.add_test_group(KogeBugtS(mpas_core=self)) - self.add_test_group(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(Thwaites(mpas_core=self))