Merge pull request #262 from cbegeman/add-omega-sphere-transport

Add support for omega sphere transport tests

Add omega support for the sphere_transport test suite

polaris/ocean/config/single_layer.yaml

@@ -9,6 +9,5 @@ mpas-ocean:
     config_disable_thick_vadv: true
     config_disable_vel_vadv: true
     config_disable_vel_vmix: true
-    config_disable_tr_all_tend: true
   pressure_gradient:
     config_pressure_gradient_type: ssh_gradient

polaris/ocean/convergence/analysis.py

@@ -150,11 +150,17 @@ def run(self):
         """
         plt.switch_backend('Agg')
         convergence_vars = self.convergence_vars
+        variables_failed = []
         for var in convergence_vars:
-            self.plot_convergence(
+            convergence_failed = self.plot_convergence(
                 variable_name=var["name"],
                 title=var["title"],
                 zidx=var["zidx"])
+            if convergence_failed:
+                variables_failed.append(var["name"])
+        if len(variables_failed) >= 1:
+            raise ValueError('Convergence rate below minimum tolerance for '
+                             f'variables {", ".join(variables_failed)}.')
 
     def plot_convergence(self, variable_name, title, zidx):
         """
@@ -281,8 +287,7 @@ def plot_convergence(self, variable_name, title, zidx):
                          f'{conv_thresh}')
             convergence_failed = True
 
-        if convergence_failed:
-            raise ValueError('Convergence rate below minimum tolerance.')
+        return convergence_failed
 
     def compute_error(self, refinement_factor, variable_name, zidx=None,
                       error_type='l2'):

polaris/ocean/suites/omega_pr.txt

@@ -1,3 +1,4 @@
 ocean/planar/manufactured_solution/convergence_both
+ocean/spherical/icos/rotation_2d/convergence_both
 ocean/spherical/icos/cosine_bell/decomp
 ocean/spherical/icos/cosine_bell/restart

polaris/ocean/tasks/sphere_transport/analysis.py

@@ -39,13 +39,13 @@ def __init__(self, component, subdir, case_name,
         self.case_name = case_name
         convergence_vars = [{'name': 'tracer1',
                              'title': 'tracer1',
-                             'zidx': 1},
+                             'zidx': 0},
                             {'name': 'tracer2',
                              'title': 'tracer2',
-                             'zidx': 1},
+                             'zidx': 0},
                             {'name': 'tracer3',
                              'title': 'tracer3',
-                             'zidx': 1}]
+                             'zidx': 0}]
         super().__init__(component=component, subdir=subdir,
                          dependencies=dependencies,
                          convergence_vars=convergence_vars,
@@ -71,8 +71,14 @@ def convergence_parameters(self, field_name=None):
             The maximum convergence rate
         """
         config = self.config
+        model = config.get('ocean', 'model')
+        if model == 'omega' and self.case_name == 'rotation_2d':
+            order = 2
+        else:
+            order = 3
         section = config[self.case_name]
-        conv_thresh = section.getfloat(f'convergence_thresh_{field_name}')
+        conv_thresh = section.getfloat(
+            f'convergence_thresh_{field_name}_order{order}')
 
         section = config['convergence']
         error_type = section.get('error_type')

polaris/ocean/tasks/sphere_transport/correlated_tracers_2d.cfg

@@ -9,9 +9,9 @@ vel_amp = 10.
 correlation_coefficients = -0.8, 0.0, 0.9
 
 # convergence threshold below which the test fails
-convergence_thresh_tracer1 = 1.5
-convergence_thresh_tracer2 = 1.3
-convergence_thresh_tracer3 = 1.0
+convergence_thresh_tracer1_order3 = 1.5
+convergence_thresh_tracer2_order3 = 1.3
+convergence_thresh_tracer3_order3 = 1.0
 
 # time in days at which to evaluate mixing
 mixing_evaluation_time = 6.0

polaris/ocean/tasks/sphere_transport/divergent_2d.cfg

@@ -5,6 +5,6 @@
 vel_amp = 5.
 
 # convergence threshold below which the test fails
-convergence_thresh_tracer1 = 1.8
-convergence_thresh_tracer2 = 1.5
-convergence_thresh_tracer3 = 0.3
+convergence_thresh_tracer1_order3 = 1.8
+convergence_thresh_tracer2_order3 = 1.5
+convergence_thresh_tracer3_order3 = 0.3

polaris/ocean/tasks/sphere_transport/forward.py

@@ -60,3 +60,13 @@ def __init__(self, component, name, subdir, base_mesh, init,
                          graph_target=f'{base_mesh.path}/graph.info',
                          refinement_factor=refinement_factor,
                          refinement=refinement)
+
+    def setup(self):
+        """
+        TEMP: symlink initial condition to name hard-coded in Omega
+        """
+        super().setup()
+        model = self.config.get('ocean', 'model')
+        # TODO: remove as soon as Omega no longer hard-codes this file
+        if model == 'omega':
+            self.add_input_file(filename='OmegaMesh.nc', target='init.nc')

polaris/ocean/tasks/sphere_transport/forward.yaml

@@ -1,16 +1,18 @@
-mpas-ocean:
-  run_modes:
-    config_ocean_run_mode: forward
+ocean:
   time_management:
     config_stop_time: none
     config_run_duration: {{ run_duration }}
+  time_integration:
+    config_dt: {{ dt }}
+    config_time_integrator: {{ time_integrator }}
+
+mpas-ocean:
+  run_modes:
+    config_ocean_run_mode: forward
   decomposition:
     config_block_decomp_file_prefix: graph.info.part.
   advection:
     config_vert_coord_movement: impermeable_interfaces
-  time_integration:
-    config_dt: {{ dt }}
-    config_time_integrator: {{ time_integrator }}
   debug:
     config_disable_thick_sflux: true
     config_disable_vel_all_tend: true
@@ -57,3 +59,31 @@ mpas-ocean:
       - refLayerThickness
       - kineticEnergyCell
       - relativeVorticityCell
+
+Omega:
+  Tendencies:
+    ThicknessFluxTendencyEnable : false
+    PVTendencyEnable : false
+    KETendencyEnable : false
+    SSHTendencyEnable : false
+    VelDiffTendencyEnable: false
+    VelHyperDiffTendencyEnable: false
+    TracerHorzAdvTendencyEnable : true
+    TracerDiffTendencyEnable : false
+    TracerHyperDiffTendencyEnable : false
+  Dimension:
+    NVertLevels: 1
+  IOStreams:
+    InitialState:
+      Filename: init.nc
+      Contents:
+        - State
+        - Tracers
+    History:
+      Filename: output.nc
+      Freq: {{ output_freq }}
+      FreqUnits: seconds
+      Contents:
+        - Tracers
+        - LayerThickness
+        - NormalVelocity

polaris/ocean/tasks/sphere_transport/init.py

@@ -1,8 +1,7 @@
 import numpy as np
 import xarray as xr
-from mpas_tools.io import write_netcdf
 
-from polaris import Step
+from polaris.ocean.model import OceanIOStep
 from polaris.ocean.tasks.sphere_transport.resources.flow_types import (
     flow_divergent,
     flow_nondivergent,
@@ -17,7 +16,7 @@
 from polaris.ocean.vertical import init_vertical_coord
 
 
-class Init(Step):
+class Init(OceanIOStep):
     """
     A step for an initial condition for for the cosine bell test case
     """
@@ -156,4 +155,4 @@ def run(self):
         ds['fEdge'] = xr.zeros_like(ds_mesh.xEdge)
         ds['fVertex'] = xr.zeros_like(ds_mesh.xVertex)
 
-        write_netcdf(ds, 'initial_state.nc')
+        self.write_model_dataset(ds, 'initial_state.nc')

polaris/ocean/tasks/sphere_transport/nondivergent_2d.cfg

@@ -5,9 +5,9 @@
 vel_amp = 10.
 
 # convergence threshold below which the test fails
-convergence_thresh_tracer1 = 1.5
-convergence_thresh_tracer2 = 1.1
-convergence_thresh_tracer3 = 0.3
+convergence_thresh_tracer1_order3 = 1.5
+convergence_thresh_tracer2_order3 = 1.1
+convergence_thresh_tracer3_order3 = 0.3
 
 # time in days at which to evaluate filament preservation
 filament_evaluation_time = 6.0

polaris/ocean/tasks/sphere_transport/rotation_2d.cfg

@@ -4,7 +4,12 @@
 # rotation vector in cartesian coordinates
 rotation_vector = 0.2, 0.7, 1.0
 
-# convergence threshold below which the test fails
-convergence_thresh_tracer1 = 1.8
-convergence_thresh_tracer2 = 2.0
-convergence_thresh_tracer3 = 0.4
+# convergence threshold below which the test fails for order 3
+convergence_thresh_tracer1_order3 = 1.8
+convergence_thresh_tracer2_order3 = 2.0
+convergence_thresh_tracer3_order3 = 0.4
+
+# convergence threshold below which the test fails for order 2
+convergence_thresh_tracer1_order2 = 1.0
+convergence_thresh_tracer2_order2 = 1.6
+convergence_thresh_tracer3_order2 = 0.3

polaris/ocean/tasks/sphere_transport/sphere_transport.cfg

@@ -5,7 +5,7 @@
 grid_type = uniform
 
 # Number of vertical levels
-vert_levels = 3
+vert_levels = 1
 
 # Depth of the bottom of the ocean
 bottom_depth = 300.0
@@ -33,7 +33,10 @@ error_type = l2
 # config options for spherical convergence tests
 [convergence_forward]
 
-# time integrator: {'split_explicit', 'RK4'}
+# time integrator
+#  either: {'RK4'}
+#  mpas-ocean: {'split_explicit'}
+#  omega: {'Forward-Backward', 'RungeKutta2'}
 time_integrator = RK4
 
 # RK4 time step per resolution (s/km), since dt is proportional to resolution

polaris/ocean/tasks/sphere_transport/viz.py

@@ -1,12 +1,12 @@
 import cmocean  # noqa: F401
-import xarray as xr
+import numpy as np
 
-from polaris import Step
-from polaris.mpas import time_index_from_xtime
+from polaris.mpas import time_since_start
+from polaris.ocean.model import OceanIOStep
 from polaris.viz import plot_global_mpas_field
 
 
-class Viz(Step):
+class Viz(OceanIOStep):
     """
     A step for plotting fields from the cosine bell output
 
@@ -69,27 +69,32 @@ def run(self):
         Run this step of the test case
         """
         config = self.config
+        model = config.get('ocean', 'model')
         mesh_name = self.mesh_name
         run_duration = config.getfloat('convergence_forward',
                                        'run_duration')
 
         variables_to_plot = self.variables_to_plot
-        ds_init = xr.open_dataset('initial_state.nc')
+        ds_init = self.open_model_dataset('initial_state.nc')
         ds_init = ds_init[variables_to_plot.keys()].isel(Time=0, nVertLevels=0)
 
-        ds_out = xr.open_dataset('output.nc')
+        ds_out = self.open_model_dataset('output.nc', decode_times=False)
         s_per_hour = 3600.0
 
         # Visualization at halfway around the globe (provided run duration is
         # set to the time needed to circumnavigate the globe)
-        tidx = time_index_from_xtime(ds_out.xtime.values,
-                                     run_duration * s_per_hour / 2.)
+        time = run_duration * s_per_hour
+        if model == 'mpas-o':
+            dt = time_since_start(ds_out.xtime.values)
+        else:
+            # time is seconds since the start of the simulation in Omega
+            dt = ds_out.Time.values
+        tidx = np.argmin(np.abs(dt - time / 2.))
         ds_mid = ds_out[variables_to_plot.keys()].isel(Time=tidx,
                                                        nVertLevels=0)
 
         # Visualization at all the way around the globe
-        tidx = time_index_from_xtime(ds_out.xtime.values,
-                                     run_duration * s_per_hour)
+        tidx = np.argmin(np.abs(dt - time))
         ds_final = ds_out[variables_to_plot.keys()].isel(Time=tidx,
                                                          nVertLevels=0)