Merge pull request E3SM-Project#72 from cbegeman/add-manufactured-sol…

…ution-test

Add manufactured solution test group

This PR adds the manufactured solution test group and the convergence test case as described in [Bishnu et al.(2023)](https://doi.org/10.22541/essoar.167100170.03833124/v1).

docs/developers_guide/ocean/api.md

@@ -130,6 +130,37 @@
    viz.Viz.run
 ```
 
+### manufactured_solution
+
+```{eval-rst}
+.. currentmodule:: polaris.ocean.tests.manufactured_solution
+
+.. autosummary::
+   :toctree: generated/
+
+   ManufacturedSolution
+
+   analysis.Analysis
+   analysis.Analysis.run
+
+   exact_solution.ExactSolution
+   exact_solution.ExactSolution.ssh
+   exact_solution.ExactSolution.normalVelocity
+
+   forward.Forward
+   forward.Forward.compute_cell_count
+   forward.Forward.dynamic_model_config
+
+   initial_state.InitialState
+   initial_state.InitialState.run
+
+   viz.Viz
+   viz.Viz.run
+
+   convergence.Convergence
+   convergence.Convergence.validate
+```
+
 ### single_column
 
 ```{eval-rst}

docs/developers_guide/ocean/test_groups/index.md

@@ -8,5 +8,6 @@
 baroclinic_channel
 global_convergence
 inertial_gravity_wave
+manufactured_solution
 single_column
 ```

docs/developers_guide/ocean/test_groups/manufactured_solution.md

@@ -0,0 +1,82 @@
+(dev-ocean-manufactured_solution)=
+
+# manufactured_solution
+
+The `manufactured_solution` test group
+({py:class}`polaris.ocean.tests.manufactured_solution.ManufacturedSolution`)
+implements a test case according to the Method of Manufactured Solutions
+(see {ref}`ocean-manufactured-solution`) at 4 resolutions (200, 100, 50, and 25 km).  Here,
+we describe the shared framework for this test group and the 1 test case.
+
+(dev-ocean-manufactured-solution)=
+
+## framework
+
+The shared config options for the `manufactured_solution` test group
+are described in {ref}`ocean-manufactured-solution` in the User's Guide.
+
+Additionally, the test group has a shared `forward.yaml` file with
+a few common model config options related to run duration and default 
+horizontal  and vertical momentum and tracer diffusion, as well as defining 
+`mesh`, `input`, `restart`, and `output` streams.
+
+### exact_solution
+
+The class {py:class}`polaris.ocean.tests.manufactured_solution.exact_solution.ExactSolution`
+defines a class for storing attributes and methods relevant to computing the
+exact manufactured solution.  The constructor obtains the parameters from the
+config file. The
+{py:meth}`polaris.ocean.tests.manufactured_solution.exact_solution.ExactSolution.ssh(t)`
+method computes the SSH field at time `t`.  The
+{py:meth}`polaris.ocean.tests.manufactured_solution.exact_solution.ExactSolution.normalVelocity(t)`
+method computes the `normalVelocity` field at time `t`.
+
+### initial_state
+
+The class {py:class}`polaris.ocean.tests.manufactured_solution.initial_state.InitialState`
+defines a step for setting up the initial state for each test case.
+
+First, a mesh appropriate for the resolution is generated using
+{py:func}`mpas_tools.planar_hex.make_planar_hex_mesh()`.  Then, the mesh is
+culled to remove periodicity in the y direction.  A vertical grid is generated
+with 1 layer.  Finally, the initial layerThickness field is computed from the
+exact solution for the SSH field and the initial velocity is also assigned to
+the exact solution. The tracer and coriolis fields are uniform in space.
+
+### forward
+
+The class {py:class}`polaris.ocean.tests.manufactured_solution.forward.Forward`
+defines a step for running MPAS-Ocean from the initial condition produced in
+the `initial_state` step.  Namelist and streams files are updated in
+{py:meth}`polaris.ocean.tests.manufactured_solution.forward.Forward.dynamic_model_config()`
+with time steps determined algorithmically based on config options.  The
+number of cells is computed from config options in
+{py:meth}`polaris.ocean.tests.manufactured_solution.forward.Forward.compute_cell_count()`
+so that this can be used to constrain the number of MPI tasks that tests
+have as their target and minimum (if the resources are not explicitly
+prescribed).  For MPAS-Ocean, PIO namelist options are modified and a
+graph partition is generated as part of `runtime_setup()`.  Finally, the ocean 
+model is run.
+
+### analysis
+
+The class {py:class}`polaris.ocean.tests.manufactured_solution.analysis.Analysis`
+defines a step for computing the root mean-square-error from the final
+simulated field and the exact solution. It uses the config options to determine
+whether the convergence rate falls within acceptable bounds.
+
+### viz
+
+The class {py:class}`polaris.ocean.tests.manufactured_solution.viz.Viz`
+defines a step for visualization. It produces two plots: the convergence of the
+RMSE with resolution and a plan-view of the simulated, exact, and (simulated -
+exact) SSH fields.
+
+(dev-ocean-manufactured-solution-convergence)=
+
+### convergence
+
+The {py:class}`polaris.ocean.tests.manufactured_solution.convergence.Convergence`
+test performs a 10-hour run.  Then, validation of `temperature`, 
+`layerThickness` and `normalVelocity` are performed against a
+baseline if one is provided when calling {ref}`dev-polaris-setup`.

docs/users_guide/ocean/test_groups/index.md

@@ -8,5 +8,6 @@
 baroclinic_channel
 global_convergence
 inertial_gravity_wave
+manufactured_solution
 single_column
 ```

docs/users_guide/ocean/test_groups/manufactured_solution.md

@@ -0,0 +1,123 @@
+(ocean-manufactured-solution)=
+
+# manufactured_solution
+
+The manufactured solution test group implements configurations for surface wave
+propagation with the rotating, nonlinear shallow water equations on a doubly
+periodic domain. This test group is intended to utilize tendency terms embedded
+in the forward ocean model in order to produce the manufactured solution. This
+solution can be then used to assess the numerical accuracy and convergence of
+the discretized nonlinear momentum equation. 
+
+Currently, the test group contains one test case, which is a convergence test
+from
+[Bishnu et al.(2023)](https://doi.org/10.22541/essoar.167100170.03833124/v1)
+
+(ocean-manufactured-solution-convergence)=
+
+## convergence
+
+### description
+
+The `convergence` test case runs the manufactured solution simulation for 4
+different resolutions: 200, 100, 50, and 25 km.
+
+The forward step for each resolution runs the simulation for 10 hours. The
+model is configured without vertical advection and mixing. No tracers are enabled
+and the pressure gradient used is the gradient of the sea surface height.
+Horizontal mixing and bottom friction are also neglected.
+
+The analysis step computes the root mean-square-error of the difference between
+the simulated SSH field and the exact solution at the end of the simulation. It
+also computes the convergence rate with resolution
+
+The visualization step produces two plots: the convergence of the RMSE with
+resolution and a plan-view of the simulated, exact, and (simulated-exact)
+SSH fields.
+
+### mesh
+
+For each resolution, the `initial_state` step generates and planar hexagonal
+mesh that is periodic in both the x and y directions.
+
+### vertical grid
+
+Since this test case is a shallow water case, the vertical grid is set to a
+single layer configuration.
+
+```cfg
+[vertical_grid]
+
+# The type of vertical grid
+grid_type = uniform
+
+# Number of vertical levels
+vert_levels = 1 
+
+# Depth of the bottom of the ocean
+bottom_depth = 1000.0
+
+# The type of vertical coordinate (e.g. z-level, z-star)
+coord_type = z-star
+
+# Whether to use "partial" or "full", or "None" to not alter the topography
+partial_cell_type = None
+
+# The minimum fraction of a layer for partial cells
+min_pc_fraction = 0.1
+```
+
+### initial conditions
+
+The initial conditions are set to the following:
+$$
+\eta = \eta_0 \sin(k_x x + k_y y - \omega t)\\
+u = \eta_0 \cos(k_x x + k_y y - \omega t)\\
+v = u
+$$
+
+### forcing
+
+N/A
+
+### time step and run duration
+
+The time step is determined by the config option ``dt_per_km`` according to the
+mesh resolution. The run duration is 10 hours.
+
+### config options
+
+The following config options are availiable for this case:
+
+```cfg
+[manufactured_solution]
+
+# the size of the domain in km in the x and y directions
+lx = 10000.0
+
+# the coriolis parameter
+coriolis_parameter = 1.0e-4
+
+# the amplitude of the sea surface height perturbation
+ssh_amplitude = 1.0
+
+# Number of wavelengths in x direction
+n_wavelengths_x = 2
+
+# Number of wavelengths in y direction
+n_wavelengths_y = 2
+
+# Time step per resolution (s/km), since dt is proportional to resolution
+dt_per_km = 3.0
+
+# Convergence threshold below which the test fails
+conv_thresh = 1.8
+
+# Convergence rate above which a warning is issued
+conv_max = 2.2
+```
+
+### cores
+
+The number of cores is determined according to the config options
+``max_cells_per_core`` and ``goal_cells_per_core``.

polaris/ocean/__init__.py

@@ -2,6 +2,7 @@
 from polaris.ocean.tests.baroclinic_channel import BaroclinicChannel
 from polaris.ocean.tests.global_convergence import GlobalConvergence
 from polaris.ocean.tests.inertial_gravity_wave import InertialGravityWave
+from polaris.ocean.tests.manufactured_solution import ManufacturedSolution
 from polaris.ocean.tests.single_column import SingleColumn
 
 
@@ -20,6 +21,7 @@ def __init__(self):
         self.add_test_group(BaroclinicChannel(component=self))
         self.add_test_group(GlobalConvergence(component=self))
         self.add_test_group(InertialGravityWave(component=self))
+        self.add_test_group(ManufacturedSolution(component=self))
         self.add_test_group(SingleColumn(component=self))
 
     def configure(self, config):

polaris/ocean/suites/nightly.txt

@@ -2,3 +2,4 @@ ocean/baroclinic_channel/10km/threads_test
 ocean/baroclinic_channel/10km/decomp_test
 ocean/baroclinic_channel/10km/restart_test
 ocean/inertial_gravity_wave/convergence
+ocean/manufactured_solution/convergence

polaris/ocean/suites/pr.txt

@@ -4,3 +4,4 @@ ocean/baroclinic_channel/10km/restart_test
 ocean/inertial_gravity_wave/convergence
 ocean/single_column/960km/cvmix
 ocean/single_column/960km/ideal_age
+ocean/manufactured_solution/convergence

polaris/ocean/tests/manufactured_solution/__init__.py

@@ -0,0 +1,16 @@
+from polaris import TestGroup
+from polaris.ocean.tests.manufactured_solution.convergence import Convergence
+
+
+class ManufacturedSolution(TestGroup):
+    """
+    A test group for manufactured solution test cases
+    """
+    def __init__(self, component):
+        """
+        component : polaris.ocean.Ocean
+            the ocean component that this test group belongs to
+        """
+        super().__init__(component=component, name='manufactured_solution')
+
+        self.add_test_case(Convergence(test_group=self))

polaris/ocean/tests/manufactured_solution/analysis.py

@@ -0,0 +1,81 @@
+import datetime
+import warnings
+
+import cmocean  # noqa: F401
+import numpy as np
+import xarray as xr
+
+from polaris import Step
+from polaris.ocean.tests.manufactured_solution.exact_solution import (
+    ExactSolution,
+)
+
+
+class Analysis(Step):
+    """
+    A step for analysing the output from the manufactured solution
+    test case
+
+    Attributes
+    ----------
+    resolutions : list of int
+        The resolutions of the meshes that have been run
+    """
+    def __init__(self, test_case, resolutions):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        test_case : polaris.TestCase
+            The test case this step belongs to
+
+        resolutions : list of int
+            The resolutions of the meshes that have been run
+        """
+        super().__init__(test_case=test_case, name='analysis')
+        self.resolutions = resolutions
+
+        for resolution in resolutions:
+            self.add_input_file(
+                filename=f'init_{resolution}km.nc',
+                target=f'../{resolution}km/initial_state/initial_state.nc')
+            self.add_input_file(
+                filename=f'output_{resolution}km.nc',
+                target=f'../{resolution}km/forward/output.nc')
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        config = self.config
+        resolutions = self.resolutions
+
+        section = config['manufactured_solution']
+        conv_thresh = section.getfloat('conv_thresh')
+        conv_max = section.getfloat('conv_max')
+
+        rmse = []
+        for i, res in enumerate(resolutions):
+            init = xr.open_dataset(f'init_{res}km.nc')
+            ds = xr.open_dataset(f'output_{res}km.nc')
+            exact = ExactSolution(config, init)
+
+            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()
+            ssh_model = ds.ssh.values[-1, :]
+            rmse.append(np.sqrt(np.mean((ssh_model - exact.ssh(t).values)**2)))
+
+        p = np.polyfit(np.log10(resolutions), np.log10(rmse), 1)
+        conv = p[0]
+
+        if conv < conv_thresh:
+            raise ValueError(f'order of convergence '
+                             f' {conv} < min tolerence {conv_thresh}')
+
+        if conv > conv_max:
+            warnings.warn(f'order of convergence '
+                          f'{conv} > max tolerence {conv_max}')