Extend manufactured solution to del2 and del4

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.

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
 

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

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

polaris/ocean/tasks/manufactured_solution/__init__.py

@@ -35,14 +35,23 @@ 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):
     """
     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)

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')

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

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
 

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

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,39 +138,39 @@ 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')
             ds_init = self.open_model_dataset(
                 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))