diff --git a/compass/ocean/__init__.py b/compass/ocean/__init__.py index 929720e50f..7bf9296c18 100644 --- a/compass/ocean/__init__.py +++ b/compass/ocean/__init__.py @@ -1,5 +1,6 @@ from compass.mpas_core import MpasCore from compass.ocean.tests.baroclinic_channel import BaroclinicChannel +from compass.ocean.tests.buttermilk_bay import ButtermilkBay from compass.ocean.tests.dam_break import DamBreak from compass.ocean.tests.drying_slope import DryingSlope from compass.ocean.tests.global_convergence import GlobalConvergence @@ -37,6 +38,7 @@ def __init__(self): super().__init__(name='ocean') self.add_test_group(BaroclinicChannel(mpas_core=self)) + self.add_test_group(ButtermilkBay(mpas_core=self)) self.add_test_group(DamBreak(mpas_core=self)) self.add_test_group(DryingSlope(mpas_core=self)) self.add_test_group(GlobalConvergence(mpas_core=self)) diff --git a/compass/ocean/tests/buttermilk_bay/__init__.py b/compass/ocean/tests/buttermilk_bay/__init__.py new file mode 100644 index 0000000000..74e80d4160 --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/__init__.py @@ -0,0 +1,18 @@ +from compass.ocean.tests.buttermilk_bay.default import Default +from compass.testgroup import TestGroup + + +class ButtermilkBay(TestGroup): + """ + A test group for Buttermilk Bay (subgrid 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='buttermilk_bay') + for wetdry in ['standard', 'subgrid']: + self.add_test_case( + Default(test_group=self, wetdry=wetdry)) diff --git a/compass/ocean/tests/buttermilk_bay/buttermilk_bay.cfg b/compass/ocean/tests/buttermilk_bay/buttermilk_bay.cfg new file mode 100644 index 0000000000..852bc42e33 --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/buttermilk_bay.cfg @@ -0,0 +1,33 @@ +[job] + +wall_time = 2:00:00 + + +# config options for buttermilk bay +[buttermilk_bay] + +# dimensions of domain in x and y directions (m) +Lx = 4608 +Ly = 4608 + +# a list of resolutions (m) to test +resolutions = 256, 128, 64 + +# time step per resolution (s/m), since dt is proportional to resolution +dt_per_m = 0.02 + +# 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 +[buttermilk_bay_viz] + +# coordinates (in km) for timeseries plot +points = [2.8, 0.53], [1.9, 1.66], [2.4, 3.029], [2.51, 3.027], [1.26, 1.56] + +# generate contour plots at a specified interval between output timesnaps +plot_interval = 1 diff --git a/compass/ocean/tests/buttermilk_bay/default/__init__.py b/compass/ocean/tests/buttermilk_bay/default/__init__.py new file mode 100644 index 0000000000..549ce1cff9 --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/default/__init__.py @@ -0,0 +1,141 @@ +import numpy as np + +from compass.config import CompassConfigParser +from compass.ocean.tests.buttermilk_bay.forward import Forward +from compass.ocean.tests.buttermilk_bay.initial_state import InitialState +from compass.ocean.tests.buttermilk_bay.viz import Viz +from compass.testcase import TestCase +from compass.validate import compare_variables + + +class Default(TestCase): + """ + The default buttermilk_bay test case + + Attributes + ---------- + wetdry : str + The type of wetting and drying (``standard``, ``subgrid``) + """ + + def __init__(self, test_group, wetdry): + """ + Create the test case + + Parameters + ---------- + test_group : compass.ocean.tests.buttermilk_bay.ButtermilkBay + The test group that this test case belongs to + + wetdry : str + The type of wetting and drying used (``standard``, ``subgrid``) + """ + name = wetdry + subdir = wetdry + super().__init__(test_group=test_group, name=name, + subdir=subdir) + + self.resolutions = None + self.wetdry = wetdry + # add the steps with default resolutions so they can be listed + config = CompassConfigParser() + config.add_from_package('compass.ocean.tests.buttermilk_bay', + 'buttermilk_bay.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('buttermilk_bay', + 'goal_cells_per_core') + max_cells_per_core = config.getfloat('buttermilk_bay', + 'max_cells_per_core') + lx = config.getfloat('buttermilk_bay', 'Lx') + ly = config.getfloat('buttermilk_bay', 'Ly') + + for resolution in self.resolutions: + + 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}m' + step = self.steps[f'forward_{res_name}'] + step.ntasks = ntasks + step.min_tasks = min_tasks + + config.set('buttermilk_bay', f'{res_name}_ntasks', str(ntasks), + comment=f'Target core count for {res_name} mesh') + config.set('buttermilk_bay', 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 = '256, 128, 64' + + # set the default values that a user may change before setup + config.set('buttermilk_bay', 'resolutions', default_resolutions, + comment='a list of resolutions (m) to test') + + # get the resolutions back, perhaps with values set in the user's + # config file + resolutions = config.getlist('buttermilk_bay', + '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}m' + + init_step = InitialState(test_case=self, + name=f'initial_state_{res_name}', + resolution=resolution, + wetdry=self.wetdry) + self.add_step(init_step) + self.add_step(Forward(test_case=self, + name=f'forward_{res_name}', + resolution=resolution, + wetdry=self.wetdry)) + self.add_step(Viz(test_case=self, + wetdry=self.wetdry, + 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}m/output.nc') diff --git a/compass/ocean/tests/buttermilk_bay/forward.py b/compass/ocean/tests/buttermilk_bay/forward.py new file mode 100644 index 0000000000..92d38d14c4 --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/forward.py @@ -0,0 +1,122 @@ +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 buttermilk bay + test cases. + """ + def __init__(self, test_case, resolution, + name, + 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 + """ + + self.resolution = resolution + + super().__init__(test_case=test_case, name=name) + + self.add_namelist_file('compass.ocean.tests.buttermilk_bay', + 'namelist.forward') + + res_name = f'{resolution}m' + self.add_namelist_file('compass.ocean.tests.buttermilk_bay', + f'namelist.{coord_type}.forward') + + if wetdry == 'subgrid': + self.add_namelist_file('compass.ocean.tests.buttermilk_bay', + 'namelist.subgrid.forward') + self.add_streams_file('compass.ocean.tests.buttermilk_bay', + 'streams.subgrid.forward') + else: + self.add_streams_file('compass.ocean.tests.buttermilk_bay', + '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_input_file(filename='forcing.nc', + target=f'{input_path}/init_mode_forcing_data.nc') + + 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_m + 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_m = config.getfloat('buttermilk_bay', 'dt_per_m') + + dt = dt_per_m * 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('buttermilk_bay', + f'{self.resolution}m_ntasks') + self.min_tasks = config.getint('buttermilk_bay', + f'{self.resolution}m_min_tasks') + self.openmp_threads = 1 diff --git a/compass/ocean/tests/buttermilk_bay/initial_state.py b/compass/ocean/tests/buttermilk_bay/initial_state.py new file mode 100644 index 0000000000..e9a8b12e93 --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/initial_state.py @@ -0,0 +1,89 @@ +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 buttermilk bay test + cases + """ + def __init__(self, test_case, name, resolution, + coord_type='single_layer', wetdry='standard'): + """ + Create the step + + Parameters + ---------- + test_case : compass.ocean.tests.buttermilk_bay.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.buttermilk_bay', + 'namelist.init', mode='init') + + if wetdry == 'subgrid': + self.add_namelist_file('compass.ocean.tests.buttermilk_bay', + 'namelist.subgrid.init', mode='init') + + self.add_streams_file('compass.ocean.tests.buttermilk_bay', + 'streams.init', mode='init') + + self.add_input_file( + filename='buttermilk_bathy.nc', + target='buttermilk_bathy.nc', + database='bathymetry_database') + + for file in ['base_mesh.nc', 'culled_mesh.nc', 'culled_graph.info', + 'ocean.nc', 'init_mode_forcing_data.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_buttermilk_bay_vert_levels': '1'} + else: + vert_levels = config.getint('vertical_grid', 'vert_levels') + options = {'config_buttermilk_bay_vert_levels': f'{vert_levels}'} + self.update_namelist_at_runtime(options) + + # Determine mesh parameters + Lx = config.getint('buttermilk_bay', 'Lx') + Ly = config.getint('buttermilk_bay', 'Ly') + nx = round(Lx / self.resolution) + ny = round(Ly / self.resolution) + dc = 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/buttermilk_bay/namelist.forward b/compass/ocean/tests/buttermilk_bay/namelist.forward new file mode 100644 index 0000000000..1e278c4f80 --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/namelist.forward @@ -0,0 +1,23 @@ +config_run_duration='0002_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=0.05 +config_zero_drying_velocity=.true. +config_verify_not_dry=.false. +config_thickness_flux_type='upwind' +config_use_cvmix=.false. +config_check_ssh_consistency=.true. +config_implicit_constant_bottom_drag_coeff = 0.003 +config_use_ssh_gradient_wetting_drying = .true. +config_use_mom_del2 = .true. +config_mom_del2 = 2.0e1 +config_use_tidal_forcing = .true. +config_tidal_forcing_type = 'direct' +config_tidal_forcing_model = 'monochromatic' +config_zero_drying_velocity_ramp = .true. +config_zero_drying_velocity_ramp_hmin = 0.05 +config_zero_drying_velocity_ramp_hmax = 0.10 +config_drying_safety_height = 2.0e-10 diff --git a/compass/ocean/tests/buttermilk_bay/namelist.init b/compass/ocean/tests/buttermilk_bay/namelist.init new file mode 100644 index 0000000000..6600f40fc1 --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/namelist.init @@ -0,0 +1,9 @@ +config_init_configuration = 'buttermilk_bay' +config_ocean_run_mode = 'init' +config_use_wetting_drying = .true. +config_drying_min_cell_height = 0.05 +config_write_cull_cell_mask = .false. +config_Buttermilk_bay_vert_levels = 1 +config_Buttermilk_bay_topography_file = 'buttermilk_bathy.nc' +config_write_cull_cell_mask = .true. +config_use_tidal_forcing = .true. diff --git a/compass/ocean/tests/buttermilk_bay/namelist.ramp.forward b/compass/ocean/tests/buttermilk_bay/namelist.ramp.forward new file mode 100644 index 0000000000..83673b1213 --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/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/buttermilk_bay/namelist.single_layer.forward b/compass/ocean/tests/buttermilk_bay/namelist.single_layer.forward new file mode 100644 index 0000000000..fb70c887ae --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/namelist.single_layer.forward @@ -0,0 +1,6 @@ +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_pressure_gradient_type = 'ssh_gradient' diff --git a/compass/ocean/tests/buttermilk_bay/namelist.subgrid.forward b/compass/ocean/tests/buttermilk_bay/namelist.subgrid.forward new file mode 100644 index 0000000000..476e68bc77 --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/namelist.subgrid.forward @@ -0,0 +1 @@ +config_use_subgrid_wetting_drying = .true. diff --git a/compass/ocean/tests/buttermilk_bay/namelist.subgrid.init b/compass/ocean/tests/buttermilk_bay/namelist.subgrid.init new file mode 100644 index 0000000000..85e8a4cca6 --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/namelist.subgrid.init @@ -0,0 +1,6 @@ +config_use_subgrid_wetting_drying = .true. +config_subgrid_table_levels = -1 +config_Buttermilk_bay_subgrid_table_levels = 100 +config_Buttermilk_bay_subgrid_refinement_level = 10 +config_Buttermilk_bay_subgrid_use_thin_layer = .false. +config_Buttermilk_bay_subgrid_edge_bathymetry_max_pixel = .true. diff --git a/compass/ocean/tests/buttermilk_bay/streams.forward b/compass/ocean/tests/buttermilk_bay/streams.forward new file mode 100644 index 0000000000..7458065f52 --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/streams.forward @@ -0,0 +1,51 @@ + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/compass/ocean/tests/buttermilk_bay/streams.init b/compass/ocean/tests/buttermilk_bay/streams.init new file mode 100644 index 0000000000..2f3d87dd20 --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/streams.init @@ -0,0 +1,55 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/compass/ocean/tests/buttermilk_bay/streams.subgrid.forward b/compass/ocean/tests/buttermilk_bay/streams.subgrid.forward new file mode 100644 index 0000000000..26d7675b5f --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/streams.subgrid.forward @@ -0,0 +1,74 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/compass/ocean/tests/buttermilk_bay/viz/__init__.py b/compass/ocean/tests/buttermilk_bay/viz/__init__.py new file mode 100644 index 0000000000..300c9f61a7 --- /dev/null +++ b/compass/ocean/tests/buttermilk_bay/viz/__init__.py @@ -0,0 +1,159 @@ +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 buttermilk bay results + + Attributes + ---------- + wetdry : str + The wetting and drying approach used + + resolutions : list + The grid resolutions run for this case + """ + def __init__(self, test_case, wetdry, resolutions): + """ + Create the step + + Parameters + ---------- + test_case : compass.TestCase + The test case this step belongs to + + wetdry : str + The wetting and drying approach used + + resolutions : list + The grid resolutions run for this case + """ + super().__init__(test_case=test_case, name='viz') + + self.resolutions = resolutions + self.wetdry = wetdry + + for res in resolutions: + self.add_input_file(filename=f'output_{res}m.nc', + target=f'../forward_{res}m/output.nc') + + def run(self): + """ + Run this step of the test case + """ + + points = self.get_points() + self.timeseries_plots(points) + self.contour_plots(points) + + def get_points(self): + """ + Get the point coordinates for plotting solution timeseries + """ + + points = self.config.get('buttermilk_bay_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, + figsize=(6, 2 * len(points))) + + for res in self.resolutions: + ds = xr.open_dataset(f'output_{res}m.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}m') + + 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.suptitle(f'Buttermilk Bay ({self.wetdry})') + fig.tight_layout() + fig.subplots_adjust(bottom=0.2) + fig.legend(lines, labels, + loc='lower center', ncol=4) + fig.savefig('points.png') + + 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]}m.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('buttermilk_bay_viz', + 'plot_interval') + for i, tstep in enumerate(time): + + if i % plot_interval != 0: + continue + + ncols = len(self.resolutions) + 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}m.nc') + cm = 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}m resolution') + ax[j].set_xlabel('x (km)') + ax[j].set_ylabel('y (km)') + ds.close() + + ax[j].set_aspect('equal', 'box') + ax[j].scatter(points[:, 0] / 1000, + points[:, 1] / 1000, 15, 'k') + + 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'Buttermilk Bay ({self.wetdry}) ' + f'ssh solution at t={t} days')) + fig.savefig(f'solution_{i:03d}.png') + plt.close() diff --git a/compass/ocean/tests/parabolic_bowl/__init__.py b/compass/ocean/tests/parabolic_bowl/__init__.py index 23b24fed93..854a702e15 100644 --- a/compass/ocean/tests/parabolic_bowl/__init__.py +++ b/compass/ocean/tests/parabolic_bowl/__init__.py @@ -13,12 +13,15 @@ def __init__(self, mpas_core): the MPAS core that this test group belongs to """ super().__init__(mpas_core=mpas_core, name='parabolic_bowl') - for wetdry in ['standard']: # , 'subgrid']: + for wetdry in ['standard', 'subgrid']: for ramp_type in ['ramp', 'noramp']: # note: LTS has only standard W/D for use_lts in [True, False]: - self.add_test_case( - Default(test_group=self, - ramp_type=ramp_type, - wetdry=wetdry, - use_lts=use_lts)) + if use_lts and wetdry == 'subgrid': + continue + else: + self.add_test_case( + Default(test_group=self, + ramp_type=ramp_type, + wetdry=wetdry, + use_lts=use_lts)) diff --git a/compass/ocean/tests/parabolic_bowl/default/__init__.py b/compass/ocean/tests/parabolic_bowl/default/__init__.py index c1b48cf8fc..4e9d7d7e4d 100644 --- a/compass/ocean/tests/parabolic_bowl/default/__init__.py +++ b/compass/ocean/tests/parabolic_bowl/default/__init__.py @@ -17,6 +17,15 @@ class Default(TestCase): ---------- ramp_type : str The type of vertical coordinate (``ramp``, ``noramp``, etc.) + + wetdry : str + The type of wetting and drying used (``standard``, ``subgrid``) + + use_lts : bool + Whether local time-stepping is used + + resolutions : list + List of resolutions to run """ def __init__(self, test_group, ramp_type, wetdry, use_lts): @@ -152,7 +161,10 @@ def _setup_steps(self, config, use_lts): resolution=resolution, ramp_type=self.ramp_type, wetdry=self.wetdry)) - self.add_step(Viz(test_case=self, resolutions=resolutions, + self.add_step(Viz(test_case=self, + ramp_type=self.ramp_type, + wetdry=self.wetdry, + resolutions=resolutions, use_lts=use_lts)) def validate(self): diff --git a/compass/ocean/tests/parabolic_bowl/forward.py b/compass/ocean/tests/parabolic_bowl/forward.py index 5863bd9de2..b7228f5e76 100644 --- a/compass/ocean/tests/parabolic_bowl/forward.py +++ b/compass/ocean/tests/parabolic_bowl/forward.py @@ -38,6 +38,9 @@ def __init__(self, test_case, resolution, ramp_type : str, optional Vertical coordinate configuration + + wetdry : str, optional + The wetting and drying approach uesd """ self.resolution = resolution @@ -54,7 +57,9 @@ def __init__(self, test_case, resolution, if ramp_type == 'ramp': self.add_namelist_file('compass.ocean.tests.parabolic_bowl', 'namelist.ramp.forward') - + if wetdry == 'subgrid': + self.add_namelist_file('compass.ocean.tests.parabolic_bowl', + 'namelist.subgrid.forward') if use_lts: self.add_namelist_options( {'config_time_integrator': "'LTS'", @@ -72,8 +77,12 @@ def __init__(self, test_case, resolution, target=f'{input_path}/lts_ocean.nc') else: - self.add_streams_file('compass.ocean.tests.parabolic_bowl', - 'streams.forward') + if wetdry == 'subgrid': + self.add_streams_file('compass.ocean.tests.parabolic_bowl', + 'streams.subgrid.forward') + else: + 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') diff --git a/compass/ocean/tests/parabolic_bowl/initial_state.py b/compass/ocean/tests/parabolic_bowl/initial_state.py index a182a824ec..87860eec15 100644 --- a/compass/ocean/tests/parabolic_bowl/initial_state.py +++ b/compass/ocean/tests/parabolic_bowl/initial_state.py @@ -20,6 +20,18 @@ def __init__(self, test_case, name, resolution, ---------- test_case : compass.ocean.tests.parabolic_bowl.default.Default The test case this step belongs to + + name : str + The name of the test case + + resolution : float + The grid resolution of the test case + + coord_type : str + The type vertical coordinate used + + wetdry : str + The wetting and drying approach used """ self.coord_type = coord_type self.resolution = resolution @@ -30,6 +42,10 @@ def __init__(self, test_case, name, resolution, self.add_namelist_file('compass.ocean.tests.parabolic_bowl', 'namelist.init', mode='init') + if wetdry == 'subgrid': + self.add_namelist_file('compass.ocean.tests.parabolic_bowl', + 'namelist.subgrid.init', mode='init') + self.add_streams_file('compass.ocean.tests.parabolic_bowl', 'streams.init', mode='init') diff --git a/compass/ocean/tests/parabolic_bowl/namelist.forward b/compass/ocean/tests/parabolic_bowl/namelist.forward index a8ec1ff4c4..d6567604cf 100644 --- a/compass/ocean/tests/parabolic_bowl/namelist.forward +++ b/compass/ocean/tests/parabolic_bowl/namelist.forward @@ -9,5 +9,6 @@ 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_check_ssh_consistency=.false. config_implicit_constant_bottom_drag_coeff = 0.0 +config_use_ssh_gradient_wetting_drying = .true. diff --git a/compass/ocean/tests/parabolic_bowl/namelist.subgrid.forward b/compass/ocean/tests/parabolic_bowl/namelist.subgrid.forward new file mode 100644 index 0000000000..476e68bc77 --- /dev/null +++ b/compass/ocean/tests/parabolic_bowl/namelist.subgrid.forward @@ -0,0 +1 @@ +config_use_subgrid_wetting_drying = .true. diff --git a/compass/ocean/tests/parabolic_bowl/namelist.subgrid.init b/compass/ocean/tests/parabolic_bowl/namelist.subgrid.init new file mode 100644 index 0000000000..c56e921b51 --- /dev/null +++ b/compass/ocean/tests/parabolic_bowl/namelist.subgrid.init @@ -0,0 +1,6 @@ +config_use_subgrid_wetting_drying = .true. +config_subgrid_table_levels = -1 +config_parabolic_bowl_subgrid_table_levels = 100 +config_parabolic_bowl_subgrid_refinement_level = 10 +config_parabolic_bowl_subgrid_use_thin_layer = .true. +config_parabolic_bowl_subgrid_edge_bathymetry_max_pixel = .true. diff --git a/compass/ocean/tests/parabolic_bowl/streams.init b/compass/ocean/tests/parabolic_bowl/streams.init index 1efd479999..f7d2ce0bd1 100644 --- a/compass/ocean/tests/parabolic_bowl/streams.init +++ b/compass/ocean/tests/parabolic_bowl/streams.init @@ -26,6 +26,21 @@ + + + + + + + + + + + + + + + diff --git a/compass/ocean/tests/parabolic_bowl/streams.subgrid.forward b/compass/ocean/tests/parabolic_bowl/streams.subgrid.forward new file mode 100644 index 0000000000..8a8afb12e7 --- /dev/null +++ b/compass/ocean/tests/parabolic_bowl/streams.subgrid.forward @@ -0,0 +1,61 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/compass/ocean/tests/parabolic_bowl/viz/__init__.py b/compass/ocean/tests/parabolic_bowl/viz/__init__.py index 26bc2b9de4..2c16e77d61 100644 --- a/compass/ocean/tests/parabolic_bowl/viz/__init__.py +++ b/compass/ocean/tests/parabolic_bowl/viz/__init__.py @@ -18,7 +18,7 @@ class Viz(Step): Attributes ---------- """ - def __init__(self, test_case, resolutions, use_lts): + def __init__(self, test_case, ramp_type, wetdry, resolutions, use_lts): """ Create the step @@ -31,6 +31,8 @@ def __init__(self, test_case, resolutions, use_lts): self.resolutions = resolutions self.use_lts = use_lts + self.wetdry = wetdry + self.ramp_type = ramp_type for res in resolutions: self.add_input_file(filename=f'output_{res}km.nc', @@ -93,6 +95,7 @@ def timeseries_plots(self, points): if i == len(points) - 1: lines, labels = ax[i].get_legend_handles_labels() + fig.suptitle(f'{self.wetdry} ({self.ramp_type})') fig.tight_layout() fig.subplots_adjust(bottom=0.2) fig.legend(lines, labels, @@ -185,7 +188,8 @@ def contour_plots(self, points): 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.suptitle((f'{self.wetdry} ({self.ramp_type}) ' + f'ssh solution at t={t} days')) fig.savefig(f'solution_{i:03d}.png') plt.close() @@ -202,7 +206,9 @@ def rmse_plots(self): comparisons = [] cases = {'standard_ramp': f'../../../standard/{ramp_name}/viz', - 'standard_noramp': f'../../../standard/{noramp_name}/viz'} + 'standard_noramp': f'../../../standard/{noramp_name}/viz', + 'subgrid_ramp': f'../../../subgrid/{ramp_name}/viz', + 'subgrid_noramp': f'../../../subgrid/{noramp_name}/viz'} for case in cases: include = True for res in self.resolutions: @@ -233,7 +239,7 @@ def rmse_plots(self): 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), + ax.loglog(self.resolutions, rmse_1st_order, linestyle='-', color='k', alpha=.25, label='1st order') ax.set_xlabel('Cell size (km)') diff --git a/docs/developers_guide/ocean/api.rst b/docs/developers_guide/ocean/api.rst index 733aa0f21f..66afd72a84 100644 --- a/docs/developers_guide/ocean/api.rst +++ b/docs/developers_guide/ocean/api.rst @@ -53,6 +53,35 @@ baroclinic_channel initial_state.InitialState.setup initial_state.InitialState.run +buttermilk_bay +~~~~~~~~~~~~~~ + +.. currentmodule:: compass.ocean.tests.buttermilk_bay + +.. autosummary:: + :toctree: generated/ + + ButtermilkBay + + 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.contour_plots + dam_break ~~~~~~~~~ diff --git a/docs/developers_guide/ocean/test_groups/buttermilk_bay.rst b/docs/developers_guide/ocean/test_groups/buttermilk_bay.rst new file mode 100644 index 0000000000..eb5fe0a826 --- /dev/null +++ b/docs/developers_guide/ocean/test_groups/buttermilk_bay.rst @@ -0,0 +1,46 @@ +.. _dev_ocean_buttermilk_bay: + +buttermilk_bay +================== + +The ``buttermilk_bay`` test group +(:py:class:`compass.ocean.tests.buttermilk_bay.ButtermilkBay`) +implements a realistic regional test for wetting and drying in +a multi-bay system. + +.. _dev_ocean_buttermilk_bay_default: + +default +------- + +The :py:class:`compass.ocean.tests.buttermilk_bay.default.Default` +test performs a series of 2 day-long runs where an initially still bay +is forced at the boundary with a diurnal tidal signal. +The resolution of the problem (by default, between 64 and 256 m) is +varied to demonstrate the accuracy improvement of a given wetting +and drying approach. See :ref:`ocean_buttermilk_bay_default`. +for config options and more details on the test case. + +init +~~~~ + +The class :py:class:`compass.ocean.tests.buttermilk_bay.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.buttermilk_bay.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_m`` config option. Other namelist options are taken +from the test case's ``namelist.forward``. + +viz +~~~ + +The class :py:class:`compass.ocean.tests.buttermilk_bay.viz.Viz` +defines a step for creating contour plots of the solution at different +times (``solution_***.png``), and time series plots at different points +(``points.png``) diff --git a/docs/developers_guide/ocean/test_groups/index.rst b/docs/developers_guide/ocean/test_groups/index.rst index 11302227cf..258bda2911 100644 --- a/docs/developers_guide/ocean/test_groups/index.rst +++ b/docs/developers_guide/ocean/test_groups/index.rst @@ -8,6 +8,7 @@ Test groups :titlesonly: baroclinic_channel + buttermilk_bay dam_break drying_slope global_convergence diff --git a/docs/users_guide/ocean/test_groups/buttermilk_bay.rst b/docs/users_guide/ocean/test_groups/buttermilk_bay.rst new file mode 100644 index 0000000000..48bddec58c --- /dev/null +++ b/docs/users_guide/ocean/test_groups/buttermilk_bay.rst @@ -0,0 +1,140 @@ +.. _ocean_buttermilk_bay: + +buttermilk_bay +============== + +The ``buttermilk_bay`` test group implements a realistic regional +simulation of tidally forced circulation in a multi-bay system +connected by narrow channels. The computational domain is +`Buttermilk Bay, Massachusetts `_. + +.. _ocean_buttermilk_bay_default: + +default +------- + +The ``default`` test case implements the Buttermilk Bay test case found in +`Kennedy et al. (2019) `_. +The tidal forcing is applied to a portion of the bottom boundary with an +amplitude of 2 m and a period of 0.5 days. + +By default, the resolution is varied from 256 m to 64 m by doubling the resolution, +with the time step proportional to resolution. +The wetting/drying ramp feature (``config_zero_drying_velocity_ramp = .true.``) +is used by default for both the standard and subgrid versions of this test case. +The result of the ``viz`` step of the test case is plots of the solution at +different times, a time series at various points, and a convergence plot. + +standard +~~~~~~~~ +The standard version of the test case uses the standard wetting and drying +approach. At coarser resolutions, the bathymetric representation lacks +sufficient detail to resolve the connectivity between the upper and lower +bays. As a consequence, the results of the simulation are quite poor: + +.. image:: images/buttermilk_bay_standard_solution_169.png + :width: 500 px + :align: center + +.. image:: images/buttermilk_bay_standard_points.png + :width: 500 px + :align: center + +subgrid +~~~~~~~ +The subgrid version of the test case includes a subgrid correction scheme +to incorporate the effects of subgrid scale bathymetric variations on +the flow. This allows the small-scale connectivity between bays to be +represented, leading to significant improvements over the standard +wetting and drying approach: + +.. image:: images/buttermilk_bay_subgrid_solution_169.png + :width: 500 px + :align: center + +.. image:: images/buttermilk_bay_subgrid_points.png + :width: 500 px + :align: center + +config options +~~~~~~~~~~~~~~ + +The ``buttermilk_bay`` config options include: + +.. code-block:: cfg + + # config options for buttermilk bay + [buttermilk_bay] + + # dimensions of domain in x and y directions (m) + Lx = 4608 + Ly = 4608 + + # a list of resolutions (m) to test + resolutions = 256, 128, 64 + + # time step per resolution (s/m), since dt is proportional to resolution + dt_per_m = 0.02 + + # 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 + [buttermilk_bay_viz] + + # coordinates (in km) for timeseries plot + points = [2.8, 0.53], [1.9, 1.66], [2.4, 3.029], [2.51, 3.027], [1.26, 1.56] + + # generate contour plots at a specified interval between output timesnaps + plot_interval = 1 + + +resolutions +~~~~~~~~~~~ + +The default resolutions (in m) used in the test case are: + +.. code-block:: cfg + + resolutions = 256, 128, 64 + +To alter the resolutions used in this test, you will need to create your own +config file (or add a ``buttermilk_bay`` 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 meters. If you specify a different list +before setting up ``buttermilk_bay``, 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.) The resolution value along the values +of ``Lx`` and ``Ly`` are used to determine the number of cells in the x and y +directions used to generate the mesh. + +time step +~~~~~~~~~ + +The time step for forward integration is determined by multiplying the +resolution by ``dt_per_m``, 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``. By default, these +are set to the locations used in Kennedy et al. 2019. Also, the interval +between contour plot time snaps can be controlled with ``plot_interval``. diff --git a/docs/users_guide/ocean/test_groups/images/buttermilk_bay_standard_points.png b/docs/users_guide/ocean/test_groups/images/buttermilk_bay_standard_points.png new file mode 100644 index 0000000000..5a1012100c Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/buttermilk_bay_standard_points.png differ diff --git a/docs/users_guide/ocean/test_groups/images/buttermilk_bay_standard_solution_169.png b/docs/users_guide/ocean/test_groups/images/buttermilk_bay_standard_solution_169.png new file mode 100644 index 0000000000..170f7ef48b Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/buttermilk_bay_standard_solution_169.png differ diff --git a/docs/users_guide/ocean/test_groups/images/buttermilk_bay_subgrid_points.png b/docs/users_guide/ocean/test_groups/images/buttermilk_bay_subgrid_points.png new file mode 100644 index 0000000000..edf9aeb397 Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/buttermilk_bay_subgrid_points.png differ diff --git a/docs/users_guide/ocean/test_groups/images/buttermilk_bay_subgrid_solution_169.png b/docs/users_guide/ocean/test_groups/images/buttermilk_bay_subgrid_solution_169.png new file mode 100644 index 0000000000..9e94de7056 Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/buttermilk_bay_subgrid_solution_169.png differ diff --git a/docs/users_guide/ocean/test_groups/images/parabolic_bowl_subgrid_error.png b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_subgrid_error.png new file mode 100644 index 0000000000..2a8f51e152 Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_subgrid_error.png differ diff --git a/docs/users_guide/ocean/test_groups/images/parabolic_bowl_subgrid_points.png b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_subgrid_points.png new file mode 100644 index 0000000000..f9a4dc7054 Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_subgrid_points.png differ diff --git a/docs/users_guide/ocean/test_groups/images/parabolic_bowl_subgrid_solution_360.png b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_subgrid_solution_360.png new file mode 100644 index 0000000000..4047e9f5b3 Binary files /dev/null and b/docs/users_guide/ocean/test_groups/images/parabolic_bowl_subgrid_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 4b3ce264b3..56e55dbbcc 100644 --- a/docs/users_guide/ocean/test_groups/index.rst +++ b/docs/users_guide/ocean/test_groups/index.rst @@ -12,6 +12,7 @@ coming months. :titlesonly: baroclinic_channel + buttermilk_bay dam_break drying_slope global_convergence diff --git a/docs/users_guide/ocean/test_groups/parabolic_bowl.rst b/docs/users_guide/ocean/test_groups/parabolic_bowl.rst index 2093f4d677..8ab4ff0309 100644 --- a/docs/users_guide/ocean/test_groups/parabolic_bowl.rst +++ b/docs/users_guide/ocean/test_groups/parabolic_bowl.rst @@ -72,6 +72,42 @@ different times, a time series at various points, and a convergence plot. :width: 500 px :align: center +lts +~~~ + +Both the ``ramp`` and ``noramp`` test cases can be run with the ``lts`` variant +which uses local time-stepping (LTS) as time integrator. Note that the tests +verify the ability of the LTS scheme to run correctly with wetting and drying +and are not designed to leverage the LTS capability of producing faster runs. + +subgrid +~~~~~~~ + +Both the ``ramp`` and ``noramp`` test cases can be run with a subgrid scale +correction scheme that accounts for the effects of subgrid scale flow in +partially wet cells due to fine-scale bathymetric variation. This approach +is useful becuase it allows connectivity due to unresolved features +to be represented. In many coastal applications, this scheme can enable +coarse resolution models to capture flooding with accuracy comprable to +what higher-resolution simulations achieve without the subgrid corrections. +Using the subgrid corrections for the parabolic bowl should result in +reduced errors vs. using only the standard wetting and drying. For more +details on the subgrid scale correction scheme see: +`Kennedy et al. (2019) `_. +Results for the subgrid test cases are shown below: + +.. image:: images/parabolic_bowl_subgrid_solution_360.png + :width: 500 px + :align: center + +.. image:: images/parabolic_bowl_subgrid_points.png + :width: 500 px + :align: center + +.. image:: images/parabolic_bowl_subgrid_error.png + :width: 500 px + :align: center + config options ~~~~~~~~~~~~~~ @@ -124,8 +160,8 @@ The ``parabolic_bowl`` config options include: 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. +The first 6 options are used to control properties of the initial/analytical solution. +The remaining options are discussed below. resolutions ~~~~~~~~~~~ @@ -169,11 +205,7 @@ 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``. +An error convergence plot is also generated. Errors for the ``ramp`` +and ``noramp`` cases for both the ``standard`` and ``subgrid`` cases, +if the output exists at the time it is run. -lts -~~~ - -Both the ``ramp`` and ``noramp`` test cases can be run with the ``lts`` variant -which uses local time-stepping (LTS) as time integrator. Note that the tests -verify the ability of the LTS scheme to run correctly with wetting and drying -and are not designed to leverage the LTS capability of producing faster runs.