diff --git a/compass/ocean/__init__.py b/compass/ocean/__init__.py
index f9961ae451..929720e50f 100644
--- a/compass/ocean/__init__.py
+++ b/compass/ocean/__init__.py
@@ -13,6 +13,7 @@
from compass.ocean.tests.merry_go_round import MerryGoRound
from compass.ocean.tests.nonhydro import Nonhydro
from compass.ocean.tests.overflow import Overflow
+from compass.ocean.tests.parabolic_bowl import ParabolicBowl
from compass.ocean.tests.planar_convergence import PlanarConvergence
from compass.ocean.tests.soma import Soma
from compass.ocean.tests.sphere_transport import SphereTransport
@@ -49,6 +50,7 @@ def __init__(self):
self.add_test_group(MerryGoRound(mpas_core=self))
self.add_test_group(Nonhydro(mpas_core=self))
self.add_test_group(Overflow(mpas_core=self))
+ self.add_test_group(ParabolicBowl(mpas_core=self))
self.add_test_group(PlanarConvergence(mpas_core=self))
self.add_test_group(Soma(mpas_core=self))
self.add_test_group(SphereTransport(mpas_core=self))
diff --git a/compass/ocean/tests/parabolic_bowl/__init__.py b/compass/ocean/tests/parabolic_bowl/__init__.py
new file mode 100644
index 0000000000..167b3842bb
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/__init__.py
@@ -0,0 +1,20 @@
+from compass.ocean.tests.parabolic_bowl.default import Default
+from compass.testgroup import TestGroup
+
+
+class ParabolicBowl(TestGroup):
+ """
+ A test group for parabolic bowl (wetting-and-drying) test cases
+ """
+
+ def __init__(self, mpas_core):
+ """
+ mpas_core : compass.MpasCore
+ the MPAS core that this test group belongs to
+ """
+ super().__init__(mpas_core=mpas_core, name='parabolic_bowl')
+ for wetdry in ['standard']: # , 'subgrid']:
+ for ramp_type in ['ramp', 'noramp']:
+ self.add_test_case(
+ Default(test_group=self, ramp_type=ramp_type,
+ wetdry=wetdry))
diff --git a/compass/ocean/tests/parabolic_bowl/default/__init__.py b/compass/ocean/tests/parabolic_bowl/default/__init__.py
new file mode 100644
index 0000000000..ccb7fd437f
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/default/__init__.py
@@ -0,0 +1,146 @@
+import numpy as np
+
+from compass.config import CompassConfigParser
+from compass.ocean.tests.parabolic_bowl.forward import Forward
+from compass.ocean.tests.parabolic_bowl.initial_state import InitialState
+from compass.ocean.tests.parabolic_bowl.viz import Viz
+from compass.testcase import TestCase
+from compass.validate import compare_variables
+
+
+class Default(TestCase):
+ """
+ The default parabolic_bowl test case
+
+ Attributes
+ ----------
+ ramp_type : str
+ The type of vertical coordinate (``ramp``, ``noramp``, etc.)
+ """
+
+ def __init__(self, test_group, ramp_type, wetdry):
+ """
+ Create the test case
+
+ Parameters
+ ----------
+ test_group : compass.ocean.tests.parabolic_bowl.ParabolicBowl
+ The test group that this test case belongs to
+
+ ramp_type : str
+ The type of vertical coordinate (``ramp``, ``noramp``)
+
+ wetdry : str
+ The type of wetting and drying used (``standard``, ``subgrid``)
+ """
+ name = f'{wetdry}_{ramp_type}'
+
+ subdir = f'{wetdry}/{ramp_type}'
+ super().__init__(test_group=test_group, name=name,
+ subdir=subdir)
+
+ self.resolutions = None
+ self.wetdry = wetdry
+ self.ramp_type = ramp_type
+
+ # add the steps with default resolutions so they can be listed
+ config = CompassConfigParser()
+ config.add_from_package('compass.ocean.tests.parabolic_bowl',
+ 'parabolic_bowl.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('parabolic_bowl',
+ 'goal_cells_per_core')
+ max_cells_per_core = config.getfloat('parabolic_bowl',
+ 'max_cells_per_core')
+
+ for resolution in self.resolutions:
+
+ lx = config.getfloat('parabolic_bowl', 'Lx')
+ ly = config.getfloat('parabolic_bowl', 'Ly')
+
+ nx = 2 * int(0.5 * lx / resolution + 0.5)
+ ny = 2 * int(0.5 * ly * (2. / np.sqrt(3)) / resolution + 0.5)
+
+ approx_cells = nx * ny
+ # ideally, about 300 cells per core
+ # (make it a multiple of 4 because...it looks better?)
+ ntasks = max(1,
+ 4 * round(approx_cells / (4 * goal_cells_per_core)))
+ # In a pinch, about 3000 cells per core
+ min_tasks = max(1,
+ round(approx_cells / max_cells_per_core))
+
+ res_name = f'{resolution}km'
+ step = self.steps[f'forward_{res_name}']
+ step.ntasks = ntasks
+ step.min_tasks = min_tasks
+
+ config.set('parabolic_bowl', f'{res_name}_ntasks', str(ntasks),
+ comment=f'Target core count for {res_name} mesh')
+ config.set('parabolic_bowl', f'{res_name}_min_tasks',
+ str(min_tasks),
+ comment=f'Minimum core count for {res_name} mesh')
+
+ def _setup_steps(self, config):
+ """ setup steps given resolutions """
+
+ default_resolutions = '20, 10, 5'
+
+ # set the default values that a user may change before setup
+ config.set('parabolic_bowl', 'resolutions', default_resolutions,
+ comment='a list of resolutions (km) to test')
+
+ # get the resolutions back, perhaps with values set in the user's
+ # config file
+ resolutions = config.getlist('parabolic_bowl',
+ '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 self.resolutions:
+
+ res_name = f'{resolution}km'
+
+ self.add_step(InitialState(test_case=self,
+ name=f'initial_state_{res_name}',
+ resolution=resolution,
+ wetdry=self.wetdry))
+ self.add_step(Forward(test_case=self,
+ name=f'forward_{res_name}',
+ resolution=resolution,
+ ramp_type=self.ramp_type,
+ wetdry=self.wetdry))
+ self.add_step(Viz(test_case=self, resolutions=resolutions))
+
+ def validate(self):
+ """
+ Validate variables against a baseline
+ """
+ super().validate()
+ variables = ['layerThickness', 'normalVelocity']
+ for res in self.resolutions:
+ compare_variables(test_case=self, variables=variables,
+ filename1=f'forward_{res}km/output.nc')
diff --git a/compass/ocean/tests/parabolic_bowl/forward.py b/compass/ocean/tests/parabolic_bowl/forward.py
new file mode 100644
index 0000000000..62fdd4d49d
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/forward.py
@@ -0,0 +1,123 @@
+import time
+
+from compass.model import run_model
+from compass.step import Step
+
+
+class Forward(Step):
+ """
+ A step for performing forward MPAS-Ocean runs as part of parabolic bowl
+ test cases.
+ """
+ def __init__(self, test_case, resolution, name,
+ ramp_type='ramp', coord_type='single_layer',
+ wetdry='standard'):
+ """
+ Create a new test case
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+
+ resolution : float
+ The resolution of the test case
+
+ name : str
+ the name of the test case
+
+ subdir : str, optional
+ the subdirectory for the step. The default is ``name``
+
+ coord_type : str, optional
+ vertical coordinate configuration
+
+ ramp_type : str, optional
+ vertical coordinate configuration
+ """
+
+ self.resolution = resolution
+
+ super().__init__(test_case=test_case, name=name)
+
+ self.add_namelist_file('compass.ocean.tests.parabolic_bowl',
+ 'namelist.forward')
+
+ res_name = f'{resolution}km'
+ self.add_namelist_file('compass.ocean.tests.parabolic_bowl',
+ f'namelist.{coord_type}.forward')
+
+ if ramp_type == 'ramp':
+ self.add_namelist_file('compass.ocean.tests.parabolic_bowl',
+ 'namelist.ramp.forward')
+
+ self.add_streams_file('compass.ocean.tests.parabolic_bowl',
+ 'streams.forward')
+
+ input_path = f'../initial_state_{res_name}'
+ self.add_input_file(filename='mesh.nc',
+ target=f'{input_path}/culled_mesh.nc')
+
+ self.add_input_file(filename='init.nc',
+ target=f'{input_path}/ocean.nc')
+
+ self.add_input_file(filename='graph.info',
+ target=f'{input_path}/culled_graph.info')
+
+ self.add_model_as_input()
+
+ self.add_output_file(filename='output.nc')
+
+ def setup(self):
+ """
+ Set namelist options based on config options
+ """
+ dt = self.get_dt()
+ self.add_namelist_options({'config_dt': dt})
+ self._get_resources()
+
+ def constrain_resources(self, available_cores):
+ """
+ Update resources at runtime from config options
+ """
+ self._get_resources()
+ super().constrain_resources(available_cores)
+
+ def run(self):
+ """
+ Run this step of the testcase
+ """
+ # update dt in case the user has changed dt_per_km
+ dt = self.get_dt()
+ self.update_namelist_at_runtime(options={'config_dt': dt},
+ out_name='namelist.ocean')
+
+ run_model(self)
+
+ def get_dt(self):
+ """
+ Get the time step
+
+ Returns
+ -------
+ dt : str
+ the time step in HH:MM:SS
+ """
+ config = self.config
+ # dt is proportional to resolution
+ dt_per_km = config.getfloat('parabolic_bowl', 'dt_per_km')
+
+ dt = dt_per_km * self.resolution
+ # https://stackoverflow.com/a/1384565/7728169
+ dt = time.strftime('%H:%M:%S', time.gmtime(dt))
+
+ return dt
+
+ def _get_resources(self):
+ """ get the these properties from the config options """
+ config = self.config
+ self.ntasks = config.getint('parabolic_bowl',
+ f'{self.resolution}km_ntasks')
+ self.min_tasks = config.getint('parabolic_bowl',
+ f'{self.resolution}km_min_tasks')
+ self.openmp_threads = 1
diff --git a/compass/ocean/tests/parabolic_bowl/initial_state.py b/compass/ocean/tests/parabolic_bowl/initial_state.py
new file mode 100644
index 0000000000..a182a824ec
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/initial_state.py
@@ -0,0 +1,92 @@
+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 compass.model import run_model
+from compass.step import Step
+
+
+class InitialState(Step):
+ """
+ A step for creating a mesh and initial condition for drying slope test
+ cases
+ """
+ def __init__(self, test_case, name, resolution,
+ coord_type='single_layer', wetdry='standard'):
+ """
+ Create the step
+
+ Parameters
+ ----------
+ test_case : compass.ocean.tests.parabolic_bowl.default.Default
+ The test case this step belongs to
+ """
+ self.coord_type = coord_type
+ self.resolution = resolution
+
+ super().__init__(test_case=test_case, name=name, ntasks=1,
+ min_tasks=1, openmp_threads=1)
+
+ self.add_namelist_file('compass.ocean.tests.parabolic_bowl',
+ 'namelist.init', mode='init')
+
+ self.add_streams_file('compass.ocean.tests.parabolic_bowl',
+ 'streams.init', mode='init')
+
+ for file in ['base_mesh.nc', 'culled_mesh.nc', 'culled_graph.info',
+ 'ocean.nc']:
+ self.add_output_file(file)
+
+ self.add_model_as_input()
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+ config = self.config
+ logger = self.logger
+
+ # Set vertical levels
+ coord_type = self.coord_type
+ if coord_type == 'single_layer':
+ options = {'config_parabolic_bowl_vert_levels': '1'}
+ else:
+ vert_levels = config.getint('vertical_grid', 'vert_levels')
+ options = {'config_parabolic_bowl_vert_levels': f'{vert_levels}'}
+ self.update_namelist_at_runtime(options)
+
+ # Set test case parameters
+ eta_max = config.get('parabolic_bowl', 'eta_max')
+ depth_max = config.get('parabolic_bowl', 'depth_max')
+ coriolis = config.get('parabolic_bowl', 'coriolis_parameter')
+ omega = config.get('parabolic_bowl', 'omega')
+ gravity = config.get('parabolic_bowl', 'gravity')
+ options = {'config_parabolic_bowl_Coriolis_parameter': coriolis,
+ 'config_parabolic_bowl_eta0': eta_max,
+ 'config_parabolic_bowl_b0': depth_max,
+ 'config_parabolic_bowl_omega': omega,
+ 'config_parabolic_bowl_gravity': gravity}
+ self.update_namelist_at_runtime(options)
+
+ # Determine mesh parameters
+ Lx = config.getint('parabolic_bowl', 'Lx')
+ Ly = config.getint('parabolic_bowl', 'Ly')
+ nx = round(Lx / self.resolution)
+ ny = round(Ly / self.resolution)
+ dc = 1e3 * self.resolution
+
+ logger.info(' * Make planar hex mesh')
+ dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, nonperiodic_x=True,
+ nonperiodic_y=True)
+ logger.info(' * Completed Make planar hex mesh')
+ write_netcdf(dsMesh, 'base_mesh.nc')
+
+ logger.info(' * Cull mesh')
+ dsMesh = cull(dsMesh, logger=logger)
+ logger.info(' * Convert mesh')
+ dsMesh = convert(dsMesh, graphInfoFileName='culled_graph.info',
+ logger=logger)
+ logger.info(' * Completed Convert mesh')
+ write_netcdf(dsMesh, 'culled_mesh.nc')
+ run_model(self, namelist='namelist.ocean',
+ streams='streams.ocean')
diff --git a/compass/ocean/tests/parabolic_bowl/namelist.10km.forward b/compass/ocean/tests/parabolic_bowl/namelist.10km.forward
new file mode 100644
index 0000000000..05af815f7e
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/namelist.10km.forward
@@ -0,0 +1 @@
+config_dt='0000_00:00:05'
diff --git a/compass/ocean/tests/parabolic_bowl/namelist.20km.forward b/compass/ocean/tests/parabolic_bowl/namelist.20km.forward
new file mode 100644
index 0000000000..6594d67a91
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/namelist.20km.forward
@@ -0,0 +1 @@
+config_dt='0000_00:00:10'
diff --git a/compass/ocean/tests/parabolic_bowl/namelist.5km.forward b/compass/ocean/tests/parabolic_bowl/namelist.5km.forward
new file mode 100644
index 0000000000..b6cb4195f9
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/namelist.5km.forward
@@ -0,0 +1 @@
+config_dt='0000_00:00:02.5'
diff --git a/compass/ocean/tests/parabolic_bowl/namelist.forward b/compass/ocean/tests/parabolic_bowl/namelist.forward
new file mode 100644
index 0000000000..a8ec1ff4c4
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/namelist.forward
@@ -0,0 +1,13 @@
+config_run_duration='0003_00:00:00'
+config_time_integrator='RK4'
+config_vert_coord_movement='impermeable_interfaces'
+config_ALE_thickness_proportionality='weights_only'
+config_use_wetting_drying=.true.
+config_prevent_drying=.true.
+config_drying_min_cell_height=2.0e-2
+config_zero_drying_velocity=.true.
+config_verify_not_dry=.true.
+config_thickness_flux_type='upwind'
+config_use_cvmix=.false.
+config_check_ssh_consistency=.true.
+config_implicit_constant_bottom_drag_coeff = 0.0
diff --git a/compass/ocean/tests/parabolic_bowl/namelist.init b/compass/ocean/tests/parabolic_bowl/namelist.init
new file mode 100644
index 0000000000..d755143e04
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/namelist.init
@@ -0,0 +1,6 @@
+config_init_configuration = 'parabolic_bowl'
+config_ocean_run_mode = 'init'
+config_use_wetting_drying = .true.
+config_drying_min_cell_height = 2e-2
+config_write_cull_cell_mask = .false.
+config_parabolic_bowl_vert_levels = 1
diff --git a/compass/ocean/tests/parabolic_bowl/namelist.ramp.forward b/compass/ocean/tests/parabolic_bowl/namelist.ramp.forward
new file mode 100644
index 0000000000..83673b1213
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/namelist.ramp.forward
@@ -0,0 +1,3 @@
+config_zero_drying_velocity_ramp = .true.
+config_zero_drying_velocity_ramp_hmin = 2e-2
+config_zero_drying_velocity_ramp_hmax = 4e-2
diff --git a/compass/ocean/tests/parabolic_bowl/namelist.single_layer.forward b/compass/ocean/tests/parabolic_bowl/namelist.single_layer.forward
new file mode 100644
index 0000000000..7378b7eee7
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/namelist.single_layer.forward
@@ -0,0 +1,7 @@
+config_disable_thick_vadv = .true.
+config_disable_thick_sflux = .true.
+config_disable_vel_vmix = .true.
+config_disable_vel_vadv = .true.
+config_disable_tr_all_tend = .true.
+config_disable_vel_hmix = .true.
+config_pressure_gradient_type = 'ssh_gradient'
diff --git a/compass/ocean/tests/parabolic_bowl/parabolic_bowl.cfg b/compass/ocean/tests/parabolic_bowl/parabolic_bowl.cfg
new file mode 100644
index 0000000000..89e7ae2ca2
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/parabolic_bowl.cfg
@@ -0,0 +1,48 @@
+[job]
+
+wall_time = 2:00:00
+
+
+# config options for drying slope test cases
+[parabolic_bowl]
+
+# dimensions of domain in x and y directions (km)
+Lx = 1440
+Ly = 1560
+
+# Coriolis parameter
+coriolis_parameter = 1.031e-4
+
+# maximum initial ssh magnitude
+eta_max = 2.0
+
+# maximum water depth
+depth_max = 50.0
+
+# angular fequency of oscillation
+omega = 1.4544e-4
+
+# gravitational acceleration
+gravity = 9.81
+
+# a list of resolutions (km) to test
+resolutions = 20, 10, 5
+
+# time step per resolution (s/km), since dt is proportional to resolution
+dt_per_km = 0.5
+
+# 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
+
+# config options for visualizing drying slope ouptut
+[parabolic_bowl_viz]
+
+# coordinates (in km) for timeseries plot
+points = [0,0], [300,0], [610,0]
+
+# generate contour plots at a specified interval between output timesnaps
+plot_interval = 10
diff --git a/compass/ocean/tests/parabolic_bowl/streams.forward b/compass/ocean/tests/parabolic_bowl/streams.forward
new file mode 100644
index 0000000000..a412ee4c9f
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/streams.forward
@@ -0,0 +1,34 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/compass/ocean/tests/parabolic_bowl/streams.init b/compass/ocean/tests/parabolic_bowl/streams.init
new file mode 100644
index 0000000000..1efd479999
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/streams.init
@@ -0,0 +1,32 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/compass/ocean/tests/parabolic_bowl/viz/__init__.py b/compass/ocean/tests/parabolic_bowl/viz/__init__.py
new file mode 100644
index 0000000000..92d8bcad30
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/viz/__init__.py
@@ -0,0 +1,325 @@
+import datetime as dt
+import os
+import subprocess
+
+import matplotlib.pyplot as plt
+import numpy as np
+import xarray as xr
+from scipy.interpolate import LinearNDInterpolator
+
+from compass.step import Step
+
+
+class Viz(Step):
+ """
+ A step for visualizing parabolic bowl results and
+ comparing with analytical solution
+
+ Attributes
+ ----------
+ """
+ def __init__(self, test_case, resolutions):
+ """
+ Create the step
+
+ Parameters
+ ----------
+ test_case : compass.TestCase
+ The test case this step belongs to
+ """
+ super().__init__(test_case=test_case, name='viz')
+
+ self.resolutions = resolutions
+
+ for res in resolutions:
+ self.add_input_file(filename=f'output_{res}km.nc',
+ target=f'../forward_{res}km/output.nc')
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+
+ points = self.get_points()
+ self.timeseries_plots(points)
+ self.inject_exact_solution()
+ self.contour_plots(points)
+ self.rmse_plots()
+
+ def get_points(self):
+ """
+ Get the point coordinates for plotting solution timeseries
+ """
+
+ points = self.config.get('parabolic_bowl_viz', 'points')
+ points = points.replace('[', '').replace(']', '').split(',')
+ points = np.asarray(points, dtype=float).reshape(-1, 2)
+ points = points * 1000
+
+ return points
+
+ def timeseries_plots(self, points):
+ """
+ Plot solution timeseries at a given number of points
+ for each resolution
+ """
+
+ fig, ax = plt.subplots(nrows=len(points), ncols=1)
+
+ for res in self.resolutions:
+ ds = xr.open_dataset(f'output_{res}km.nc')
+
+ time = [dt.datetime.strptime(x.decode(), '%Y-%m-%d_%H:%M:%S')
+ for x in ds.xtime.values]
+ t = np.asarray([(x - time[0]).total_seconds() for x in time])
+
+ xy = np.vstack((ds.xCell.values, ds.yCell.values)).T
+ interp = LinearNDInterpolator(xy, ds.ssh.values.T)
+
+ for i, pt in enumerate(points):
+
+ ssh = interp(pt).T
+ ax[i].plot(t / 86400, ssh, label=f'{res}km')
+
+ for i, pt in enumerate(points):
+ ssh_exact = self.exact_solution('zeta', pt[0], pt[1], t)
+ ax[i].plot(t / 86400, ssh_exact, label='exact')
+
+ for i, pt in enumerate(points):
+ ax[i].set_xlabel('t (days)')
+ ax[i].set_ylabel('ssh (m)')
+ ax[i].set_title(f'Point ({pt[0]/1000}, {pt[1]/1000})')
+ if i == len(points) - 1:
+ lines, labels = ax[i].get_legend_handles_labels()
+
+ fig.tight_layout()
+ fig.subplots_adjust(bottom=0.2)
+ fig.legend(lines, labels,
+ loc='lower center', ncol=4)
+ fig.savefig('points.png')
+
+ def inject_exact_solution(self):
+ """
+ Save exact solution to output nc file
+ """
+
+ for res in self.resolutions:
+ ds = xr.open_dataset(f'output_{res}km.nc')
+
+ if 'ssh_exact' and 'layerThickness_exact' not in ds:
+ time = [dt.datetime.strptime(x.decode(), '%Y-%m-%d_%H:%M:%S')
+ for x in ds.xtime.values]
+ ssh_exact = ds.ssh.copy(deep=True)
+ layerThickness_exact = ds.layerThickness.copy(deep=True)
+ for i, tstep in enumerate(time):
+ t = (time[i] - time[0]).total_seconds()
+
+ ssh_exact[i, :] = self.exact_solution(
+ 'zeta', ds.xCell.values, ds.yCell.values, t)
+ layerThickness_exact[i, :, 0] = self.exact_solution(
+ 'h', ds.xCell.values, ds.yCell.values, t)
+ ds['ssh_exact'] = ssh_exact
+ ds['layerThickness_exact'] = layerThickness_exact
+ ds.ssh_exact.encoding['_FillValue'] = None
+ ds.layerThickness_exact.encoding['_FillValue'] = None
+ ds.to_netcdf(f'output_{res}km.nc',
+ format="NETCDF3_64BIT_OFFSET", mode='a')
+ ds.close()
+
+ def contour_plots(self, points):
+ """
+ Plot contour plots at a specified output interval for each resolution
+ and show where the points used in `points.png` are located.
+ """
+
+ sol_min = -2
+ sol_max = 2
+ clevels = np.linspace(sol_min, sol_max, 50)
+ cmap = plt.get_cmap('RdBu')
+
+ ds = xr.open_dataset(f'output_{self.resolutions[0]}km.nc')
+ time = [dt.datetime.strptime(x.decode(), '%Y-%m-%d_%H:%M:%S')
+ for x in ds.xtime.values]
+ ds.close()
+
+ plot_interval = self.config.getint('parabolic_bowl_viz',
+ 'plot_interval')
+ for i, tstep in enumerate(time):
+
+ if i % plot_interval != 0:
+ continue
+
+ ncols = len(self.resolutions) + 1
+ fig, ax = plt.subplots(nrows=1, ncols=ncols,
+ figsize=(5 * ncols, 5),
+ constrained_layout=True)
+
+ for j, res in enumerate(self.resolutions):
+ ds = xr.open_dataset(f'output_{res}km.nc')
+ ax[j].tricontourf(ds.xCell / 1000, ds.yCell / 1000,
+ ds['ssh'][i, :],
+ levels=clevels, cmap=cmap,
+ vmin=sol_min, vmax=sol_max, extend='both')
+ ax[j].set_aspect('equal', 'box')
+ ax[j].set_title(f'{res}km resolution')
+ ax[j].set_xlabel('x (km)')
+ ax[j].set_ylabel('y (km)')
+ ds.close()
+
+ ds = xr.open_dataset(f'output_{min(self.resolutions)}km.nc')
+ cm = ax[ncols - 1].tricontourf(ds.xCell / 1000, ds.yCell / 1000,
+ ds['ssh_exact'][i, :],
+ levels=clevels, cmap=cmap,
+ vmin=sol_min, vmax=sol_max,
+ extend='both')
+ ax[ncols - 1].set_aspect('equal', 'box')
+ ax[ncols - 1].scatter(points[:, 0] / 1000,
+ points[:, 1] / 1000, 15, 'k')
+
+ ax[ncols - 1].set_title('Analytical solution')
+ ax[ncols - 1].set_xlabel('x (km)')
+ ax[ncols - 1].set_ylabel('y (km)')
+ ds.close()
+
+ cb = fig.colorbar(cm, ax=ax[-1], shrink=0.6)
+ cb.set_label('ssh (m)')
+ t = round((time[i] - time[0]).total_seconds() / 86400., 2)
+ fig.suptitle(f'ssh solution at t={t} days')
+ fig.savefig(f'solution_{i:03d}.png')
+ plt.close()
+
+ def rmse_plots(self):
+ """
+ Plot convergence curves
+ """
+
+ comparisons = []
+ cases = {'standard_ramp': '../../../standard/ramp/viz',
+ 'standard_noramp': '../../../standard/noramp/viz'}
+ for case in cases:
+ include = True
+ for res in self.resolutions:
+ if not os.path.exists(f'{cases[case]}/output_{res}km.nc'):
+ include = False
+ if include:
+ comparisons.append(case)
+
+ fig, ax = plt.subplots(nrows=1, ncols=1)
+
+ max_rmse = 0
+ for j, comp in enumerate(comparisons):
+ rmse = np.zeros(len(self.resolutions))
+ for i, res in enumerate(self.resolutions):
+
+ rmse[i] = self.compute_rmse(
+ 'h',
+ f'{cases[comp]}/output_{res}km.nc')
+
+ if rmse[i] > max_rmse:
+ max_rmse = rmse[i]
+
+ ax.loglog(self.resolutions, rmse,
+ linestyle='-', marker='o', label=comp)
+
+ rmse_1st_order = np.zeros(len(self.resolutions))
+ rmse_1st_order[0] = max_rmse
+ for i in range(len(self.resolutions) - 1):
+ rmse_1st_order[i + 1] = rmse_1st_order[i] / 2.0
+
+ ax.loglog(self.resolutions, np.flip(rmse_1st_order),
+ linestyle='-', color='k', alpha=.25, label='1st order')
+
+ ax.set_xlabel('Cell size (km)')
+ ax.set_ylabel('RMSE (m)')
+
+ ax.legend(loc='lower right')
+ ax.set_title('Layer thickness convergence')
+ fig.tight_layout()
+ fig.savefig('error.png')
+
+ def compute_rmse(self, varname, filename):
+ """
+ Compute the rmse between the modeled and exact solutions
+ """
+
+ ds = xr.open_dataset(filename)
+
+ time = [dt.datetime.strptime(x.decode(), '%Y-%m-%d_%H:%M:%S')
+ for x in ds.xtime.values]
+ ind = time.index(dt.datetime.strptime('0001-01-03_18:00:00',
+ '%Y-%m-%d_%H:%M:%S'))
+ if varname == 'zeta':
+ var = ds['ssh'].values[ind, :]
+ elif varname == 'h':
+ var = ds['layerThickness'].values[ind, :, 0]
+
+ t = (time[ind] - time[0]).total_seconds()
+ var_exact = self.exact_solution(varname, ds.xCell.values,
+ ds.yCell.values, t)
+ rmse = np.sqrt(np.mean(np.square(var - var_exact)))
+
+ return rmse
+
+ def exact_solution(self, var, x, y, t):
+ """
+ Evaluate the exact solution
+ """
+
+ config = self.config
+
+ f = config.getfloat('parabolic_bowl', 'coriolis_parameter')
+ eta0 = config.getfloat('parabolic_bowl', 'eta_max')
+ b0 = config.getfloat('parabolic_bowl', 'depth_max')
+ omega = config.getfloat('parabolic_bowl', 'omega')
+ g = config.getfloat('parabolic_bowl', 'gravity')
+
+ x = np.array(x)
+ y = np.array(y)
+ t = np.array(t)
+
+ x = np.atleast_1d(x)
+ y = np.atleast_1d(y)
+ t = np.atleast_1d(t)
+
+ if t.size > 1:
+ x = np.resize(x, t.shape)
+ y = np.resize(y, t.shape)
+
+ eps = 1.0e-12
+ r = np.sqrt(np.square(x) + np.square(y))
+ L = np.sqrt(8.0 * g * b0 / (omega**2 - f**2))
+ C = ((b0 + eta0)**2 - b0**2) / ((b0 + eta0)**2 + b0**2)
+ b = b0 * (1.0 - r**2 / L**2)
+ num = 1.0 - C**2
+ den = 1.0 / (1.0 - C * np.cos(omega * t))
+ h = b0 * (den * np.sqrt(num) - den**2 * (r**2 / L**2) * num)
+ h[h < eps] = 0.0
+
+ if var == 'h':
+ soln = h
+
+ elif var == 'zeta':
+ soln = b0 * (den * np.sqrt(num) - 1.0 -
+ (r**2 / L**2) * (den**2 * num - 1.0))
+ soln[h < eps] = -b[h < eps]
+
+ elif var == 'u':
+ soln = 0.5 * den * (omega * x * C * np.sin(omega * t) -
+ f * y * (np.sqrt(num) +
+ C * np.cos(omega * t) - 1.0))
+ soln[h < eps] = 0
+
+ elif var == 'v':
+ soln = 0.5 * den * (omega * y * C * np.sin(omega * t) +
+ f * x * (np.sqrt(num) +
+ C * np.cos(omega * t) - 1.0))
+ soln[h < eps] = 0
+
+ elif var == 'r':
+ soln = L * np.sqrt((1.0 - C * np.cos(omega * t)) /
+ np.sqrt(1.0 - C**2))
+
+ else:
+ print('Variable name not supported')
+
+ return soln
diff --git a/docs/developers_guide/ocean/api.rst b/docs/developers_guide/ocean/api.rst
index d9e52dfe88..baf545f4d5 100644
--- a/docs/developers_guide/ocean/api.rst
+++ b/docs/developers_guide/ocean/api.rst
@@ -604,6 +604,38 @@ overflow
hydro_vs_nonhydro.visualize.Visualize
hydro_vs_nonhydro.visualize.Visualize.run
+parabolic_bowl
+~~~~~~~~~~~~~~
+
+.. currentmodule:: compass.ocean.tests.parabolic_bowl
+
+.. autosummary::
+ :toctree: generated/
+
+ ParabolicBowl
+
+ default.Default
+ default.Default.configure
+ default.Default.update_cores
+ default.Default.validate
+
+ forward.Forward
+ forward.Forward.run
+ forward.Forward.setup
+ forward.Forward.get_dt
+
+ initial_state.InitialState
+ initial_state.InitialState.run
+
+ viz.Viz
+ viz.Viz.run
+ viz.Viz.get_points
+ viz.Viz.timeseries_plots
+ viz.Viz.inject_exact_solution
+ viz.Viz.contour_plots
+ viz.Viz.rmse_plots
+ viz.Viz.compute_rmse
+ viz.Viz.exact_solution
planar_convergence
~~~~~~~~~~~~~~~~~~
diff --git a/docs/developers_guide/ocean/test_groups/index.rst b/docs/developers_guide/ocean/test_groups/index.rst
index 144df323c2..11302227cf 100644
--- a/docs/developers_guide/ocean/test_groups/index.rst
+++ b/docs/developers_guide/ocean/test_groups/index.rst
@@ -21,6 +21,7 @@ Test groups
merry_go_round
nonhydro
overflow
+ parabolic_bowl
planar_convergence
soma
sphere_transport
diff --git a/docs/developers_guide/ocean/test_groups/parabolic_bowl.rst b/docs/developers_guide/ocean/test_groups/parabolic_bowl.rst
new file mode 100644
index 0000000000..ae8ea55ee9
--- /dev/null
+++ b/docs/developers_guide/ocean/test_groups/parabolic_bowl.rst
@@ -0,0 +1,45 @@
+.. _dev_ocean_parabolic_bowl:
+
+parabolic_bowl
+==================
+
+The ``parabolic_bowl`` test group
+(:py:class:`compass.ocean.tests.parabolic_bowl.ParabolicBowl`)
+implements convergence study for wetting and drying in a parabolic bowl.
+
+.. _dev_ocean_parabolic_bowl_default:
+
+default
+-------
+
+The :py:class:`compass.ocean.tests.parabolic_bowl.default.Default`
+test performs a series of 3 day-long runs where an initial mound of water
+oscillates in a parabolic bowl. The resolution of the problem (by default, between
+5 and 20 km). Modeled results are compared with a known exact solution to
+demonstrate convergence. See :ref:`ocean_parabolic_bowl_default`.
+for config options and more details on the test case.
+
+init
+~~~~
+
+The class :py:class:`compass.ocean.tests.parabolic_bowl.init.Init`
+defines a step for setting up for generating planar hex meshes and creating
+the initial state for each test case.
+
+forward
+~~~~~~~
+
+The class :py:class:`compass.ocean.tests.parabolic_bowl.forward.Forward`
+defines a step for running MPAS-Ocean from the initial condition produced in
+the ``initial_state`` step. The time step is determined from the resolution
+based on the ``dt_per_km`` config option. Other namelist options are taken
+from the test case's ``namelist.forward``.
+
+viz
+~~~
+
+The class :py:class:`compass.ocean.tests.parabolic_bowl.viz.Viz`
+defines a step for creating contour plots of the solution at different
+times (``solution_***.png``), time series plots at different points
+(``points.png``). It also computes the root mean squared error for
+the results at each resolution and plots them in ``error.png``.
diff --git a/docs/users_guide/ocean/test_groups/images/parabolic_bowl_error.png b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_error.png
new file mode 100644
index 0000000000..e106ee9ed7
Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_error.png differ
diff --git a/docs/users_guide/ocean/test_groups/images/parabolic_bowl_points.png b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_points.png
new file mode 100644
index 0000000000..7af6ccffbf
Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_points.png differ
diff --git a/docs/users_guide/ocean/test_groups/images/parabolic_bowl_solution_000.png b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_solution_000.png
new file mode 100644
index 0000000000..04b44ef898
Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_solution_000.png differ
diff --git a/docs/users_guide/ocean/test_groups/images/parabolic_bowl_solution_360.png b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_solution_360.png
new file mode 100644
index 0000000000..f7eae1e2b0
Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_solution_360.png differ
diff --git a/docs/users_guide/ocean/test_groups/index.rst b/docs/users_guide/ocean/test_groups/index.rst
index 7a2bbaec68..4b3ce264b3 100644
--- a/docs/users_guide/ocean/test_groups/index.rst
+++ b/docs/users_guide/ocean/test_groups/index.rst
@@ -25,6 +25,7 @@ coming months.
merry_go_round
nonhydro
overflow
+ parabolic_bowl
planar_convergence
soma
sphere_transport
diff --git a/docs/users_guide/ocean/test_groups/parabolic_bowl.rst b/docs/users_guide/ocean/test_groups/parabolic_bowl.rst
new file mode 100644
index 0000000000..8030a395ef
--- /dev/null
+++ b/docs/users_guide/ocean/test_groups/parabolic_bowl.rst
@@ -0,0 +1,171 @@
+.. _ocean_parabolic_bowl:
+
+parabolic_bowl
+==============
+
+The ``parabolic_bowl`` test group implements convergence study for
+wetting and drying. Currently, the only test case is the default case
+
+.. _ocean_parabolic_bowl_default:
+
+default
+-------
+
+The ``default`` test case implements the parabolic bowl test case found in
+`Thacker 1981 `_. The problem
+consists an initial mound of water which propagates outward and reflects off
+a wet/dry boundary in a parabolic-shaped basin. The presence of a Coriolis
+factor causes the wave to rotate around the bowl as it oscillates. The
+bathymetry is given as:
+
+.. math::
+
+ b(x,y) = b_0\left(1 - \frac{R^2}{L^2}\right),
+
+with
+
+.. math::
+
+ R &= \sqrt{x^2 + y^2}, \\
+ L &= \sqrt{\frac{8gb_0}{\omega^2 - f^2}}.
+
+An exact solution for this problem exists for the frictionless, nonlinear shallow water equations:
+
+.. math::
+
+ \eta(x,y,t) &= b_0\left[ \frac{\sqrt{1-C^2}}{1-C\cos(\omega t)} -1 - \left( \frac{R^2}{L^2} \right)\left( \frac{1-C^2}{(1-C\cos(\omega t))^2}-1\right)\right], \\
+ u(x,y,t) &= \frac{1}{1-C\cos(\omega t)}\left( \frac{1}{2}\omega x C \sin(\omega t) - \frac{1}{2}fy\left( \sqrt{1-C^2} + C\cos(\omega t) - 1\right)\right), \\
+ v(x,y,t) &= \frac{1}{1-C\cos(\omega t)}\left(\frac{1}{2}fx\left(\sqrt{1-C^2} + A\cos(\omega t) -1 \right) + \frac{1}{2}\omega y C \sin(\omega t)\right),
+
+where, :math:`C` is defined as:
+
+.. math::
+
+ C = \frac{(b_0+\eta_0)^2 - b_0^2}{(b_0+\eta_0)^2 + b_0^2}.
+
+Since this is a single layer case, the solution for the total depth, :math:`h = \eta + b`, is
+
+.. math::
+
+ h(x,y,t) = b_0\left[ \frac{\sqrt{1-C^2}}{1-C\cos(\omega t)} - \left( \frac{R^2}{L^2} \right)\left( \frac{1-C^2}{(1-C\cos(\omega t))^2}\right)\right].
+
+
+By default, the resolution is varied from 20 km to 5 km by doubling the resolution,
+with the time step proportional to resolution.
+The result of the ``viz`` step of the test case is are plots of the solution at
+different times, a time series at various points, and a convergence plot.
+
+
+.. image:: images/parabolic_bowl_solution_000.png
+ :width: 500 px
+ :align: center
+
+.. image:: images/parabolic_bowl_solution_360.png
+ :width: 500 px
+ :align: center
+
+.. image:: images/parabolic_bowl_points.png
+ :width: 500 px
+ :align: center
+
+.. image:: images/parabolic_bowl_error.png
+ :width: 500 px
+ :align: center
+
+config options
+~~~~~~~~~~~~~~
+
+The ``parabolic_bowl`` config options include:
+
+.. code-block:: cfg
+
+ # config options for drying slope test cases
+ [parabolic_bowl]
+
+ # dimensions of domain in x and y directions (km)
+ Lx = 1440
+ Ly = 1560
+
+ # Coriolis parameter
+ coriolis_parameter = 1.031e-4
+
+ # Maximum initial ssh magnitude
+ eta_max = 2.0
+
+ # Maximum water depth
+ depth_max = 50.0
+
+ # Angular fequency of oscillation
+ omega = 1.4544e-4
+
+ # Gravitational acceleration
+ gravity = 9.81
+
+ # a list of resolutions (km) to test
+ resolutions = 20, 10, 5
+
+ # time step per resolution (s/km), since dt is proportional to resolution
+ dt_per_km = 0.5
+
+ # 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
+
+ # config options for visualizing drying slope ouptut
+ [parabolic_bowl_viz]
+
+ # coordinates (in km) for timeseries plot
+ points = [0,0], [300,0], [610,0]
+
+ # generate contour plots at a specified interval between output timesnaps
+ plot_interval = 10
+
+
+The last 7 options are used to control properties of the cosine bell and the
+background properties. The first 4 options are discussed below.
+
+resolutions
+~~~~~~~~~~~
+
+The default resolutions (in km) used in the test case are:
+
+.. code-block:: cfg
+
+ resolutions = 20, 10, 5
+
+To alter the resolutions used in this test, you will need to create your own
+config file (or add a ``parabolic_bowl`` section to a config file if you're
+already using one). The resolutions are a comma-separated list of the
+resolution of the mesh in km. If you specify a different list
+before setting up ``parabolic_bowl``, steps will be generated with the requested
+resolutions. (If you alter ``resolutions`` in the test case's config file in
+the work directory, nothing will happen.)
+
+time step
+~~~~~~~~~
+
+The time step for forward integration is determined by multiplying the
+resolution by ``dt_per_km``, so that coarser meshes have longer time steps.
+You can alter this before setup (in a user config file) or before running the
+test case (in the config file in the work directory).
+
+cores
+~~~~~
+
+The number of cores (and the minimum) is proportional to the number of cells,
+so that the number of cells per core is roughly constant. You can alter how
+many cells are allocated to each core with ``goal_cells_per_core``. You can
+control the maximum number of cells that are allowed to be placed on a single
+core (before the test case will fail) with ``max_cells_per_core``. If there
+aren't enough processors to handle the finest resolution, you will see that
+the step (and therefore the test case) has failed.
+
+viz
+~~~
+
+The visualization step can be configured to plot the timeseries for an
+arbitrary set of coordinates by setting ``points``. Also, the interval
+between contour plot time snaps can be controlled with ``plot_interval``.