Merge pull request #241 from cbegeman/add-barotropic-gyre-task

Add barotropic gyre task

The barotropic gyre test case is added in order to evaluate linearized barotropic ocean dynamics alongside laplacian viscosity and surface wind forcing. This is the Munk model (1950). The variant of this test case with a balance between wind stress and bottom stress due to drag (the Stommel model) is not included.

docs/developers_guide/ocean/api.md

@@ -53,6 +53,34 @@
    rpe.analysis.Analysis.run
 ```
 
+### barotropic_gyre
+
+```{eval-rst}
+.. currentmodule:: polaris.ocean.tasks.barotropic_gyre
+
+.. autosummary::
+   :toctree: generated/
+
+   add_barotropic_gyre_tasks
+
+   BarotropicGyre
+
+   analysis.Analysis
+   analysis.Analysis.compute_error
+   analysis.Analysis.exact_solution
+   analysis.Analysis.run
+
+   forward.compute_max_time_step
+   forward.Forward
+   forward.Forward.compute_cell_count
+   forward.Forward.dynamic_model_config
+
+   init.Init
+   init.Init.setup
+   init.Init.run
+
+```
+
 ### cosine_bell
 
 ```{eval-rst}

docs/developers_guide/ocean/tasks/barotropic_gyre.md

@@ -0,0 +1,59 @@
+(dev-ocean-barotropic-gyre)=
+
+# barotropic_gyre
+
+The barotropic_gyre task group is currently comprised of one `default` task.
+
+## framework
+
+The shared config options for `barotropic_gyre` tests  are described in
+{ref}`ocean-barotropic-gyre` in the User's Guide.
+
+Additionally, the tests share a `forward.yaml` file with a few common model
+config options related to time management, time integration, and Laplacian
+viscosity, as well as defining `mesh`, `input`, `restart`, and `output`
+streams.
+
+### init
+
+The class {py:class}`polaris.ocean.tasks.barotropic_gyre.init.Init`
+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 x and y directions.  A vertical grid is
+generated, with 1 layer by default.  Next, the wind stress forcing field is
+generated.
+
+### forward
+
+The class {py:class}`polaris.ocean.tasks.barotropic_gyre.forward.Forward`
+defines a step for running the ocean from the initial condition produced in
+the `init` step. Namelist and streams files are updated in
+{py:meth}`polaris.ocean.tasks.barotropic_gyre.forward.Forward.dynamic_model_config()`
+with time steps determined algorithmically based on config options.  The
+number of cells is approximated from config options in
+{py:meth}`polaris.ocean.tasks.barotropic_gyre.forward.Forward.compute_cell_count()`
+so that this can be used to constrain the number of MPI tasks that Polaris 
+tasks 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()`.  Next, the ocean 
+model is run. If `run_time_steps` is provided then this determines the run
+duration, otherwise the duration is 3 years. Finally, validation of
+`layerThickness` and `normalVelocity` in the `output.nc` file are performed 
+against a baseline if one is provided when calling {ref}`dev-polaris-setup`.
+
+### analysis
+
+The {py:class}`polaris.ocean.tasks.barotropic_gyre.analysis.Analysis`
+computes the L2 error norm at the final time step of the simulation against
+the analytical solution for the linearized dynamics. This step also produces a
+figure with the model solution, the analytical solution, and the difference
+between the two.
+
+(dev-ocean-baroclinic-channel-default)=
+
+## default
+
+The {py:class}`polaris.ocean.tasks.baroclinic_channel.default.Default`
+test performs a test of the linearized dynamics. 

docs/developers_guide/ocean/tasks/index.md

@@ -6,6 +6,7 @@
 :titlesonly: true
 
 baroclinic_channel
+barotropic_gyre
 correlated_tracers_2d
 cosine_bell
 geostrophic

docs/users_guide/ocean/tasks/barotropic_gyre.md

@@ -0,0 +1,156 @@
+(ocean-barotropic-gyre)=
+
+# barotropic gyre
+
+The barotropic gyre test case implements the Munk model for a simplified,
+wind-driven gyre
+([Munk 1950](https://doi.org/10.1175/1520-0469%281950%29007%3C0080:OTWDOC%3E2.0.CO;2)).
+It is a variant of the Stommel model except it uses Laplacian horizontal
+viscosity rather than drag to balance the wind stress. As such, this test case
+may be used to verify the Laplacian viscosity term and the wind stress forcing
+term.
+
+The test case is barotropic.
+
+Here, we linearize the momentum equations to permit comparison against an
+analytical solution. The nonlinear solution would be a logical extension of
+this case and should be qualitatively similar.
+
+## description
+
+A successful test shows the development of a western boundary current and
+clockwise circulation. The comparison with the analytical solution is performed
+via a streamfunction. The L2 norm is printed and the following figure is
+produced:
+
+```{image} images/barotropic_gyre_solution.png
+:align: center
+:width: 800 px
+```
+
+## mesh
+
+The mesh is planar with dimensions given by the config options `lx` and `ly`
+and the resolution given by the config option `resolution`. The defaults are
+1200 km x 1200 km and a resolution of 20 km but different mesh dimensions and
+resolutions may be tested. The resolution will be compared with the estimated
+lateral boundary layer width to ensure that the western boundary current will
+be adequately resolved.
+
+The boundary conditions are non-periodic in both x- and y-dimensions.
+
+## vertical grid
+
+As this test case is barotropic, the default number of vertical levels is 1.
+A different number of vertical levels may be specified by the user. The
+bottom topography is flat so this is not a test case for which different
+vertical coordinates would be informative.
+
+```cfg
+# Options related to the vertical grid
+[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 = 5000.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 condition is at rest. Since the test case is barotropic, only a
+reference density is provided with the config option `rho_0`. A beta-plane
+is indicated by the coriolis parameter and determined by the config options
+`f_0` and `beta`.
+
+## forcing
+
+This test has a constant wind stress field following:
+
+$$
+\tau_x &= \tau_0 * \cos(\pi \frac{y}{l_y}) \\
+\tau_y &= 0
+$$
+
+where `tau_0` is given by a config option.
+
+```{image} images/barotropic_gyre_forcing.png
+:align: center
+:width: 500 px
+```
+
+### config options
+
+```cfg
+[barotropic_gyre]
+
+# distance in kilometers between cell centers
+resolution = 20
+
+# Longituidinal domain length in kilometers
+lx = 1200
+
+# Latitudinal domain length in kilometers
+ly = 1200
+
+# Maximum amplitude of the zonal wind stress [N m-2]
+tau_0 = 0.1
+
+# Horizontal visocity [m2 s-1]
+nu_2 = 4e2
+
+# [s-1 m-1]
+beta = 10e-11
+
+# [s-1]
+f_0 = 10e-4
+
+# homogenous fluid density [kg m-3]
+rho_0 = 1000
+```
+
+The config option `nu_2` specifies the del2 horizontal viscosity. This value
+will be compared against the resolution to check for stability at set-up. All
+other config options are explained further in previous sections.
+
+### cores
+
+The number of cores is determined by `goal_cells_per_core` and
+`max_cells_per_core` in the `ocean` section of the config file.
+
+(ocean-barotropic-gyre-default)=
+
+## supported models
+
+These tasks support only MPAS-Ocean.
+
+## default
+
+The default case is designed such that only a short forward run is performed.
+However, when users set up the default case, a long forward run is also set
+up and may be run along with the analysis step to compare the numerical
+solution with an analytical solution. The user should change the
+`barotropic_gyre: steps_to_run` config option to include desired steps.
+
+### time step and run duration
+
+The time step for forward integration is set as a function of the CFL
+condition assuming a maximum velocity of 1 m/s and a function of the coriolis
+parameter `f_0`. The stability parameter is set to 0.25 based on the stability
+of MPAS-Ocean for this test case.
+
+The run duration is 3 time steps for the `short_forward` step and 3 years for
+the `long_forward` step.

docs/users_guide/ocean/tasks/index.md

@@ -6,6 +6,7 @@
 :titlesonly: true
 
 baroclinic_channel
+barotropic_gyre
 correlated_tracers_2d
 cosine_bell
 geostrophic

polaris/ocean/__init__.py

@@ -1,5 +1,6 @@
 from polaris import Component
 from polaris.ocean.tasks.baroclinic_channel import add_baroclinic_channel_tasks
+from polaris.ocean.tasks.barotropic_gyre import add_barotropic_gyre_tasks
 from polaris.ocean.tasks.cosine_bell import add_cosine_bell_tasks
 from polaris.ocean.tasks.geostrophic import add_geostrophic_tasks
 from polaris.ocean.tasks.ice_shelf_2d import add_ice_shelf_2d_tasks
@@ -28,6 +29,7 @@ def __init__(self):
 
         # planar: please keep these in alphabetical order
         add_baroclinic_channel_tasks(component=self)
+        add_barotropic_gyre_tasks(component=self)
         add_ice_shelf_2d_tasks(component=self)
         add_inertial_gravity_wave_tasks(component=self)
         add_internal_wave_tasks(component=self)

polaris/ocean/tasks/barotropic_gyre/__init__.py

@@ -0,0 +1,67 @@
+from polaris import Step, Task
+from polaris.config import PolarisConfigParser
+from polaris.ocean.tasks.barotropic_gyre.analysis import Analysis
+from polaris.ocean.tasks.barotropic_gyre.forward import Forward
+from polaris.ocean.tasks.barotropic_gyre.init import Init
+
+
+def add_barotropic_gyre_tasks(component):
+    """
+    Add a task that defines a convergence test for inertial gravity waves
+
+    component : polaris.ocean.Ocean
+        the ocean component that the task will be added to
+    """
+    test_name = 'default'
+    component.add_task(BarotropicGyre(component=component,
+                                      test_name=test_name))
+
+
+class BarotropicGyre(Task):
+    """
+    The convergence test case for inertial gravity waves
+    """
+
+    def __init__(self, component, test_name):
+        """
+        Create the test case
+
+        Parameters
+        ----------
+        component : polaris.ocean.Ocean
+            The ocean component that this task belongs to
+        """
+        group_name = 'barotropic_gyre'
+        name = f'{group_name}_{test_name}'
+        subdir = f'planar/{group_name}/{test_name}'
+        super().__init__(component=component, name=name, subdir=subdir)
+
+        config = self.config
+        config_filename = f'{group_name}.cfg'
+        config.filepath = f'{subdir}/{config_filename}'
+        config.add_from_package(f'polaris.ocean.tasks.{group_name}',
+                                config_filename)
+        self.set_shared_config(config, link=config_filename)
+
+        init = Init(component=component, subdir=subdir)
+        init.set_shared_config(config, link=config_filename)
+        self.add_step(init)
+
+        forward = Forward(component=component, indir=self.subdir, ntasks=None,
+                          min_tasks=None, openmp_threads=1,
+                          name='short_forward', run_time_steps=3,
+                          graph_target=f'{init.path}/culled_graph.info')
+        forward.set_shared_config(config, link=config_filename)
+        self.add_step(forward)
+
+        forward = Forward(component=component, indir=self.subdir, ntasks=None,
+                          min_tasks=None, openmp_threads=1,
+                          name='long_forward',
+                          graph_target=f'{init.path}/culled_graph.info')
+        forward.set_shared_config(config, link=config_filename)
+        self.add_step(forward, run_by_default=False)
+
+        analysis = Analysis(component=component, indir=self.subdir,
+                            boundary_condition='free slip')
+        analysis.set_shared_config(config, link=config_filename)
+        self.add_step(analysis, run_by_default=False)