diff --git a/docs/users_guide/ocean/tasks/manufactured_solution.md b/docs/users_guide/ocean/tasks/manufactured_solution.md index 866f2d85d..3d58a1ae8 100644 --- a/docs/users_guide/ocean/tasks/manufactured_solution.md +++ b/docs/users_guide/ocean/tasks/manufactured_solution.md @@ -16,9 +16,9 @@ Currently, the there is only one task, the convergence test from These tasks support both MPAS-Ocean and Omega. -(ocean-manufactured-solution-convergence)= +(ocean-manufactured-solution-default)= -## convergence +## default ### description @@ -127,3 +127,22 @@ conv_thresh = 1.8 The number of cores is determined according to the config options ``max_cells_per_core`` and ``goal_cells_per_core``. + +(ocean-manufactured-solution-del2)= + +## del2 + +### description + +All settings are the same as the {ref}`ocean-manufactured-solution-default` case +except laplacian viscosity is turned on. + +(ocean-manufactured-solution-del4)= + +## del4 + +### description + +All settings are the same as the {ref}`ocean-manufactured-solution-default` case +except hyperviscosity is turned on. The expected convergence should not be +significantly different from the {ref}`ocean-manufactured-solution-default` case. diff --git a/e3sm_submodules/E3SM-Project b/e3sm_submodules/E3SM-Project index 67ae4d5d7..3ad78d70b 160000 --- a/e3sm_submodules/E3SM-Project +++ b/e3sm_submodules/E3SM-Project @@ -1 +1 @@ -Subproject commit 67ae4d5d770e97ecf689aae239500de428d7153b +Subproject commit 3ad78d70bb60f4fc42d90b60b31b1e195bd7be75 diff --git a/polaris/ocean/convergence/analysis.py b/polaris/ocean/convergence/analysis.py index c644c02d0..f3008091a 100644 --- a/polaris/ocean/convergence/analysis.py +++ b/polaris/ocean/convergence/analysis.py @@ -218,7 +218,7 @@ def plot_convergence(self, variable_name, title, zidx): convergence_failed = False poly = np.polyfit(np.log10(refinement_array), np.log10(error_array), 1) convergence = poly[0] - conv_round = np.round(convergence, 3) + conv_round = convergence fit = refinement_array ** poly[0] * 10 ** poly[1] @@ -239,7 +239,7 @@ def plot_convergence(self, variable_name, title, zidx): ax.loglog(resolutions, order2, 'k', label='second order', alpha=0.3) ax.loglog(refinement_array, fit, 'k', - label=f'linear fit (order={conv_round})') + label=f'linear fit (order={convergence:1.3f})') ax.loglog(refinement_array, error_array, 'o', label='numerical') if self.baseline_dir is not None: @@ -255,14 +255,14 @@ def plot_convergence(self, variable_name, title, zidx): poly = np.polyfit( np.log10(refinement_array), np.log10(error_array), 1) base_convergence = poly[0] - conv_round = np.round(base_convergence, 3) + conv_round = base_convergence fit = refinement_array ** poly[0] * 10 ** poly[1] ax.loglog(refinement_array, error_array, 'o', color='#ff7f0e', label='baseline') ax.loglog(refinement_array, fit, color='#ff7f0e', label=f'linear fit, baseline ' - f'(order={conv_round})') + f'(order={base_convergence:1.3f})') ax.set_xlabel(x_label) ax.set_ylabel(f'{error_title}') ax.set_title(f'Error Convergence of {title}') @@ -273,11 +273,11 @@ def plot_convergence(self, variable_name, title, zidx): plt.close() logger.info(f'Order of convergence for {title}: ' - f'{conv_round}') + f'{conv_round:1.3f}') if conv_round < conv_thresh: logger.error(f'Error: order of convergence for {title}\n' - f' {conv_round} < min tolerance ' + f' {conv_round:1.3f} < min tolerance ' f'{conv_thresh}') convergence_failed = True @@ -410,7 +410,7 @@ def get_output_field(self, refinement_factor, field_name, time, zidx=None): ds_out = ds_out.isel(Time=tidx) field_mpas = ds_out[field_name] - if zidx is not None: + if zidx is not None and 'nVertLevels' in field_mpas.dims: field_mpas = field_mpas.isel(nVertLevels=zidx) return field_mpas diff --git a/polaris/ocean/suites/convergence.txt b/polaris/ocean/suites/convergence.txt index 9a92ef47e..af9940a77 100644 --- a/polaris/ocean/suites/convergence.txt +++ b/polaris/ocean/suites/convergence.txt @@ -1,5 +1,7 @@ ocean/planar/inertial_gravity_wave/convergence_both -ocean/planar/manufactured_solution/convergence_both +ocean/planar/manufactured_solution/convergence_both/default +ocean/planar/manufactured_solution/convergence_both/del2 +ocean/planar/manufactured_solution/convergence_both/del4 ocean/spherical/icos/correlated_tracers_2d cached: icos_base_mesh_60km icos_base_mesh_120km icos_base_mesh_240km icos_base_mesh_480km ocean/spherical/qu/correlated_tracers_2d diff --git a/polaris/ocean/suites/manufactured_solution.txt b/polaris/ocean/suites/manufactured_solution.txt new file mode 100644 index 000000000..ff95dcc96 --- /dev/null +++ b/polaris/ocean/suites/manufactured_solution.txt @@ -0,0 +1,5 @@ +ocean/planar/manufactured_solution/convergence_space/default +# ocean/planar/manufactured_solution/convergence_time/default +ocean/planar/manufactured_solution/convergence_both/default +ocean/planar/manufactured_solution/convergence_both/del2 +ocean/planar/manufactured_solution/convergence_both/del4 diff --git a/polaris/ocean/tasks/manufactured_solution/__init__.py b/polaris/ocean/tasks/manufactured_solution/__init__.py index 965fc0a9e..b6cac98c2 100644 --- a/polaris/ocean/tasks/manufactured_solution/__init__.py +++ b/polaris/ocean/tasks/manufactured_solution/__init__.py @@ -35,6 +35,14 @@ def add_manufactured_solution_tasks(component): component.add_task(ManufacturedSolution(component=component, config=config, refinement=refinement)) + component.add_task(ManufacturedSolution(component=component, + config=config, + refinement='both', + del2=True)) + component.add_task(ManufacturedSolution(component=component, + config=config, + refinement='both', + del4=True)) class ManufacturedSolution(Task): @@ -42,7 +50,8 @@ class ManufacturedSolution(Task): The convergence test case using the method of manufactured solutions """ - def __init__(self, component, config, refinement='both'): + def __init__(self, component, config, refinement='both', del2=False, + del4=False): """ Create the test case @@ -56,11 +65,33 @@ def __init__(self, component, config, refinement='both'): refinement : str, optional Whether to refine in space, time or both space and time + + del2 : bool + Whether to evaluate the momentum del2 operator + + del4 : bool + Whether to evaluate the momentum del4 operator """ - name = f'manufactured_solution_convergence_{refinement}' basedir = 'planar/manufactured_solution' subdir = f'{basedir}/convergence_{refinement}' + name = f'manufactured_solution_convergence_{refinement}' + + if del2: + test_name = 'del2' + if del4: + del4 = False + print('Manufactured solution test does not currently support' + 'both del2 and del4 convergence; testing del2 alone.') + elif del4: + test_name = 'del4' + else: + test_name = 'default' + + name = f'{name}_{test_name}' + subdir = f'{subdir}/{test_name}' + config_filename = 'manufactured_solution.cfg' + super().__init__(component=component, name=name, subdir=subdir) self.set_shared_config(config, link=config_filename) @@ -86,7 +117,7 @@ def __init__(self, component, config, refinement='both'): init_step = component.steps[subdir] else: init_step = Init(component=component, resolution=resolution, - subdir=subdir) + name=f'init_{mesh_name}', subdir=subdir) init_step.set_shared_config(self.config, link=config_filename) if resolution not in resolutions: self.add_step(init_step, symlink=symlink) @@ -97,7 +128,7 @@ def __init__(self, component, config, refinement='both'): timestep = ceil(timestep) timesteps.append(timestep) - subdir = f'{basedir}/forward/{mesh_name}_{timestep}s' + subdir = f'{basedir}/{test_name}/forward/{mesh_name}_{timestep}s' symlink = f'forward/{mesh_name}_{timestep}s' if subdir in component.steps: forward_step = component.steps[subdir] @@ -108,7 +139,8 @@ def __init__(self, component, config, refinement='both'): refinement_factor=refinement_factor, name=f'forward_{mesh_name}_{timestep}s', subdir=subdir, - init=init_step) + init=init_step, + del2=del2, del4=del4) forward_step.set_shared_config( config, link=config_filename) self.add_step(forward_step, symlink=symlink) diff --git a/polaris/ocean/tasks/manufactured_solution/forward.py b/polaris/ocean/tasks/manufactured_solution/forward.py index 2b7b9882b..e73c145b4 100644 --- a/polaris/ocean/tasks/manufactured_solution/forward.py +++ b/polaris/ocean/tasks/manufactured_solution/forward.py @@ -21,9 +21,18 @@ class Forward(ConvergenceForward): refinement : str Refinement type. One of 'space', 'time' or 'both' indicating both space and time + + resolution : float + The resolution of the test case in km + + del2 : bool + Whether to evaluate the momentum del2 operator + + del4 : bool + Whether to evaluate the momentum del4 operator """ def __init__(self, component, name, refinement_factor, subdir, - init, refinement='both'): + init, refinement='both', del2=False, del4=False): """ Create a new test case @@ -47,6 +56,12 @@ def __init__(self, component, name, refinement_factor, subdir, refinement : str, optional Refinement type. One of 'space', 'time' or 'both' indicating both space and time + + del2 : bool + Whether to evaluate the momentum del2 operator + + del4 : bool + Whether to evaluate the momentum del4 operator """ super().__init__(component=component, name=name, subdir=subdir, @@ -57,6 +72,8 @@ def __init__(self, component, name, refinement_factor, subdir, graph_target=f'{init.path}/culled_graph.info', output_filename='output.nc', validate_vars=['layerThickness', 'normalVelocity']) + self.del2 = del2 + self.del4 = del4 def setup(self): """ @@ -102,10 +119,25 @@ def dynamic_model_config(self, at_setup): super().dynamic_model_config(at_setup=at_setup) exact_solution = ExactSolution(self.config) - options = {'config_manufactured_solution_amplitude': - float(exact_solution.eta0), - 'config_manufactured_solution_wavelength_x': - float(exact_solution.lambda_x), - 'config_manufactured_solution_wavelength_y': - float(exact_solution.lambda_y)} - self.add_model_config_options(options, config_model='mpas-ocean') + mpas_options = {'config_manufactured_solution_amplitude': + float(exact_solution.eta0), + 'config_manufactured_solution_wavelength_x': + float(exact_solution.lambda_x), + 'config_manufactured_solution_wavelength_y': + float(exact_solution.lambda_y)} + shared_options = {} + if self.del2: + mpas_options['config_disable_vel_hmix'] = False + shared_options['config_use_mom_del2'] = True + shared_options['config_use_mom_del4'] = False + elif self.del4: + mpas_options['config_disable_vel_hmix'] = False + shared_options['config_use_mom_del2'] = False + shared_options['config_use_mom_del4'] = True + else: + mpas_options['config_disable_vel_hmix'] = True + shared_options['config_use_mom_del2'] = False + shared_options['config_use_mom_del4'] = False + + self.add_model_config_options(mpas_options, config_model='mpas-ocean') + self.add_model_config_options(shared_options, config_model='ocean') diff --git a/polaris/ocean/tasks/manufactured_solution/forward.yaml b/polaris/ocean/tasks/manufactured_solution/forward.yaml index 98e2ce8cb..f8d5141eb 100644 --- a/polaris/ocean/tasks/manufactured_solution/forward.yaml +++ b/polaris/ocean/tasks/manufactured_solution/forward.yaml @@ -11,10 +11,16 @@ mpas-ocean: config_bottom_drag_mode: implicit config_implicit_bottom_drag_type: constant config_implicit_constant_bottom_drag_coeff: 0.0 + hmix_del2: + config_mom_del2: 1.5e6 + hmix_del4: + config_mom_del4: 5.e13 manufactured_solution: config_use_manufactured_solution: true debug: - config_disable_vel_hmix: true + config_compute_active_tracer_budgets: false + config_disable_vel_vmix: true + config_disable_tr_all_tend: true streams: mesh: filename_template: init.nc diff --git a/polaris/ocean/tasks/manufactured_solution/init.py b/polaris/ocean/tasks/manufactured_solution/init.py index d08cf15fc..fa820b727 100644 --- a/polaris/ocean/tasks/manufactured_solution/init.py +++ b/polaris/ocean/tasks/manufactured_solution/init.py @@ -5,7 +5,6 @@ from polaris.mesh.planar import compute_planar_hex_nx_ny from polaris.ocean.model import OceanIOStep -from polaris.ocean.resolution import resolution_to_subdir from polaris.ocean.tasks.manufactured_solution.exact_solution import ( ExactSolution, ) @@ -22,7 +21,7 @@ class Init(OceanIOStep): resolution : float The resolution of the test case in km """ - def __init__(self, component, resolution, subdir): + def __init__(self, component, resolution, subdir, name): """ Create the step @@ -37,9 +36,8 @@ def __init__(self, component, resolution, subdir): taskdir : str The subdirectory that the task belongs to """ - mesh_name = resolution_to_subdir(resolution) super().__init__(component=component, - name=f'init_{mesh_name}', + name=name, subdir=subdir) self.resolution = resolution diff --git a/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg b/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg index c6b270211..479555530 100644 --- a/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg +++ b/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg @@ -25,7 +25,7 @@ n_wavelengths_x = 2 n_wavelengths_y = 2 # Time step per resolution (s/km), since dt is proportional to resolution -dt_per_km = 3.0 +dt_per_km = 1.5 # Convergence threshold below which the test fails conv_thresh = 1.8 diff --git a/polaris/ocean/tasks/manufactured_solution/viz.py b/polaris/ocean/tasks/manufactured_solution/viz.py index 1df29a89f..a87306da9 100644 --- a/polaris/ocean/tasks/manufactured_solution/viz.py +++ b/polaris/ocean/tasks/manufactured_solution/viz.py @@ -1,9 +1,8 @@ -import datetime - import cmocean # noqa: F401 import matplotlib.pyplot as plt import numpy as np +from polaris.mpas import time_since_start from polaris.ocean.convergence import ( get_resolution_for_task, get_timestep_for_task, @@ -109,7 +108,7 @@ def setup(self): forward = dependencies['forward'][refinement_factor] self.add_input_file( filename=f'mesh_r{refinement_factor:02g}.nc', - work_dir_target=f'{base_mesh.path}/base_mesh.nc') + work_dir_target=f'{base_mesh.path}/culled_mesh.nc') self.add_input_file( filename=f'init_r{refinement_factor:02g}.nc', work_dir_target=f'{init.path}/initial_state.nc') @@ -139,9 +138,13 @@ def run(self): use_mplstyle() fig, axes = plt.subplots(nrows=nres, ncols=3, figsize=(12, 2 * nres)) - rmse = [] error_range = None + section = config['convergence'] + eval_time = section.getfloat('convergence_eval_time') + s_per_hour = 3600.0 + time = eval_time * s_per_hour + for i, refinement_factor in enumerate(refinement_factors): ds_mesh = self.open_model_dataset( f'mesh_r{refinement_factor:02g}.nc') @@ -149,29 +152,25 @@ def run(self): f'init_r{refinement_factor:02g}.nc') ds = self.open_model_dataset( f'output_r{refinement_factor:02g}.nc', decode_times=False) + exact = ExactSolution(config, ds_init) if model == 'mpas-o': - t0 = datetime.datetime.strptime(ds.xtime.values[0].decode(), - '%Y-%m-%d_%H:%M:%S') - tf = datetime.datetime.strptime(ds.xtime.values[-1].decode(), - '%Y-%m-%d_%H:%M:%S') - t = (tf - t0).total_seconds() - + dt = time_since_start(ds.xtime.values) else: # time is seconds since the start of the simulation in Omega - t = ds.Time[-1].values + dt = ds.Time.values + tidx = np.argmin(np.abs(dt - time)) - ssh_model = ds.ssh.isel(Time=-1) + ssh_model = ds.ssh.isel(Time=tidx) if 'nVertLevels' in ssh_model.dims: # Omega v0 uses stacked shallow water where ssh has nVertLevels ssh_model = ssh_model.isel(nVertLevels=0) - rmse.append(np.sqrt(np.mean( - (ssh_model.values - exact.ssh(t).values)**2))) # Comparison plots - ds['ssh_exact'] = exact.ssh(t) - ds['ssh_error'] = ssh_model - exact.ssh(t) + ssh_exact = exact.ssh(time) + ds['ssh_exact'] = ssh_exact + ds['ssh_error'] = ssh_model - ssh_exact if error_range is None: error_range = np.max(np.abs(ds.ssh_error.values))