Merge pull request #249 from cbegeman/add-internal-wave-test

Add internal wave test

.gitignore

@@ -8,9 +8,6 @@ local.config
 *.o
 *.mod
 
-# Ignore all PNGs
-*.png
-
 # Each tool should also have it's own .gitignore file that ignores the build files for that tool.
 
 # Byte-compiled / optimized / DLL files

compass/ocean/__init__.py

@@ -3,6 +3,7 @@
 from compass.ocean.tests.global_convergence import GlobalConvergence
 from compass.ocean.tests.global_ocean import GlobalOcean
 from compass.ocean.tests.gotm import Gotm
+from compass.ocean.tests.internal_wave import InternalWave
 from compass.ocean.tests.ice_shelf_2d import IceShelf2d
 from compass.ocean.tests.isomip_plus import IsomipPlus
 from compass.ocean.tests.ziso import Ziso
@@ -23,6 +24,7 @@ def __init__(self):
         self.add_test_group(GlobalConvergence(mpas_core=self))
         self.add_test_group(GlobalOcean(mpas_core=self))
         self.add_test_group(Gotm(mpas_core=self))
+        self.add_test_group(InternalWave(mpas_core=self))
         self.add_test_group(IceShelf2d(mpas_core=self))
         self.add_test_group(IsomipPlus(mpas_core=self))
         self.add_test_group(Ziso(mpas_core=self))

compass/ocean/tests/internal_wave/__init__.py

@@ -0,0 +1,21 @@
+from compass.testgroup import TestGroup
+from compass.ocean.tests.internal_wave.default import Default
+from compass.ocean.tests.internal_wave.rpe_test import RpeTest
+from compass.ocean.tests.internal_wave.ten_day_test import TenDayTest
+
+
+class InternalWave(TestGroup):
+    """
+    A test group for General Ocean Turbulence Model (GOTM) test cases
+    """
+
+    def __init__(self, mpas_core):
+        """
+        mpas_core : compass.MpasCore
+            the MPAS core that this test group belongs to
+        """
+        super().__init__(mpas_core=mpas_core, name='internal_wave')
+
+        self.add_test_case(Default(test_group=self))
+        self.add_test_case(RpeTest(test_group=self))
+        self.add_test_case(TenDayTest(test_group=self))

compass/ocean/tests/internal_wave/default/__init__.py

@@ -0,0 +1,33 @@
+from compass.testcase import TestCase
+from compass.ocean.tests.internal_wave.initial_state import InitialState
+from compass.ocean.tests.internal_wave.forward import Forward
+from compass.ocean.tests.internal_wave.viz import Viz
+from compass.validate import compare_variables
+
+
+class Default(TestCase):
+    """
+    The default test case for the internal wave test
+    """
+
+    def __init__(self, test_group):
+        """
+        Create the test case
+
+        Parameters
+        ----------
+        test_group : compass.ocean.tests.internal_wave.InternalWave
+            The test group that this test case belongs to
+        """
+        super().__init__(test_group=test_group, name='default')
+        self.add_step(InitialState(test_case=self))
+        self.add_step(Forward(test_case=self, cores=4, threads=1))
+        self.add_step(Viz(test_case=self), run_by_default=False)
+
+    def validate(self):
+        """
+        Validate variables against a baseline
+        """
+        compare_variables(test_case=self,
+                          variables=['layerThickness', 'normalVelocity'],
+                          filename1='forward/output.nc')

compass/ocean/tests/internal_wave/forward.py

@@ -0,0 +1,72 @@
+from compass.model import run_model
+from compass.step import Step
+
+
+class Forward(Step):
+    """
+    A step for performing forward MPAS-Ocean runs as part of internal wave
+    test cases.
+    """
+    def __init__(self, test_case, name='forward', subdir=None, cores=1,
+                 min_cores=None, threads=1, nu=None):
+        """
+        Create a new test case
+
+        Parameters
+        ----------
+        test_case : compass.TestCase
+            The test case this step belongs to
+
+        name : str
+            the name of the test case
+
+        subdir : str, optional
+            the subdirectory for the step.  The default is ``name``
+
+        cores : int, optional
+            the number of cores the step would ideally use.  If fewer cores
+            are available on the system, the step will run on all available
+            cores as long as this is not below ``min_cores``
+
+        min_cores : int, optional
+            the number of cores the step requires.  If the system has fewer
+            than this number of cores, the step will fail
+
+        threads : int, optional
+            the number of threads the step will use
+
+        nu : float, optional
+            the viscosity (if different from the default for the test group)
+        """
+        if min_cores is None:
+            min_cores = cores
+        super().__init__(test_case=test_case, name=name, subdir=subdir,
+                         cores=cores, min_cores=min_cores, threads=threads)
+        self.add_namelist_file('compass.ocean.tests.internal_wave',
+                               'namelist.forward')
+        if nu is not None:
+            # update the viscosity to the requested value
+            options = {'config_mom_del2': '{}'.format(nu)}
+            self.add_namelist_options(options)
+
+        self.add_streams_file('compass.ocean.tests.internal_wave',
+                              'streams.forward')
+
+        self.add_input_file(filename='init.nc',
+                            target='../initial_state/ocean.nc')
+        self.add_input_file(filename='mesh.nc',
+                            target='../initial_state/culled_mesh.nc')
+        self.add_input_file(filename='graph.info',
+                            target='../initial_state/culled_graph.info')
+
+        self.add_model_as_input()
+
+        self.add_output_file(filename='output.nc')
+
+    # no setup() is needed
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        run_model(self)

compass/ocean/tests/internal_wave/initial_state.py

@@ -0,0 +1,132 @@
+import xarray
+import numpy
+
+from mpas_tools.planar_hex import make_planar_hex_mesh
+from mpas_tools.io import write_netcdf
+from mpas_tools.mesh.conversion import convert, cull
+
+from compass.ocean.vertical import init_vertical_coord
+from compass.step import Step
+
+
+class InitialState(Step):
+    """
+    A step for creating a mesh and initial condition for internal wave test
+    cases
+    """
+    def __init__(self, test_case):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        test_case : compass.ocean.tests.internal_wave.default.Default
+            The test case this step belongs to
+        """
+        super().__init__(test_case=test_case, name='initial_state', cores=1,
+                         min_cores=1, threads=1)
+
+        self.add_namelist_file('compass.ocean.tests.internal_wave',
+                               'namelist.init')
+
+        self.add_streams_file('compass.ocean.tests.internal_wave',
+                              'streams.init')
+
+        for file in ['base_mesh.nc', 'culled_mesh.nc', 'culled_graph.info',
+                     'ocean.nc']:
+            self.add_output_file(file)
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        config = self.config
+        logger = self.logger
+
+        replacements = dict()
+        replacements['config_periodic_planar_vert_levels'] = \
+            config.getfloat('vertical_grid', 'vert_levels')
+        replacements['config_periodic_planar_bottom_depth'] = \
+            config.getfloat('vertical_grid', 'bottom_depth')
+        self.update_namelist_at_runtime(options=replacements)
+
+        section = config['vertical_grid']
+        vert_levels = section.getint('vert_levels')
+        bottom_depth = section.getfloat('bottom_depth')
+
+        section = config['internal_wave']
+        nx = section.getint('nx')
+        ny = section.getint('ny')
+        dc = section.getfloat('dc')
+        use_distances = section.getboolean('use_distances')
+        amplitude_width_dist = section.getfloat('amplitude_width_dist')
+        amplitude_width_frac = section.getfloat('amplitude_width_frac')
+        bottom_temperature = section.getfloat('bottom_temperature')
+        surface_temperature = section.getfloat('surface_temperature')
+        temperature_difference = section.getfloat('temperature_difference')
+        salinity = section.getfloat('salinity')
+
+        logger.info(' * Make planar hex mesh')
+        dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, nonperiodic_x=False,
+                                      nonperiodic_y=True)
+        logger.info(' * Completed Make planar hex mesh')
+        write_netcdf(dsMesh, 'base_mesh.nc')
+
+        logger.info(' * Cull mesh')
+        dsMesh = cull(dsMesh, logger=logger)
+        logger.info(' * Convert mesh')
+        dsMesh = convert(dsMesh, graphInfoFileName='culled_graph.info',
+                         logger=logger)
+        logger.info(' * Completed Convert mesh')
+        write_netcdf(dsMesh, 'culled_mesh.nc')
+
+        ds = dsMesh.copy()
+        yCell = ds.yCell
+
+        ds['bottomDepth'] = bottom_depth * xarray.ones_like(yCell)
+        ds['ssh'] = xarray.zeros_like(yCell)
+
+        init_vertical_coord(config, ds)
+
+        yMin = yCell.min().values
+        yMax = yCell.max().values
+
+        yMid = 0.5*(yMin + yMax)
+
+        if use_distances:
+            perturbation_width = amplitude_width_dist
+        else:
+            perturbation_width = (yMax - yMin) * amplitude_width_frac
+
+        # Set stratified temperature
+        temp_vert = (bottom_temperature
+                     + (surface_temperature - bottom_temperature) *
+                       ((ds.refZMid + bottom_depth) / bottom_depth))
+
+        depth_frac = xarray.zeros_like(temp_vert)
+        refBottomDepth = ds['refBottomDepth']
+        for k in range(1, vert_levels):
+            depth_frac[k] = refBottomDepth[k-1] / refBottomDepth[vert_levels-1]
+
+        # If cell is in the southern half, outside the sin width, subtract
+        # temperature difference
+        frac = xarray.where(numpy.abs(yCell - yMid) < perturbation_width,
+                            numpy.cos(0.5 * numpy.pi * (yCell - yMid) /
+                                      perturbation_width) *
+                            numpy.sin(numpy.pi * depth_frac),
+                            0.)
+
+        temperature = temp_vert - temperature_difference * frac
+        temperature = temperature.transpose('nCells', 'nVertLevels')
+        temperature = temperature.expand_dims(dim='Time', axis=0)
+
+        normalVelocity = xarray.zeros_like(ds.xEdge)
+        normalVelocity, _ = xarray.broadcast(normalVelocity, ds.refBottomDepth)
+        normalVelocity = normalVelocity.transpose('nEdges', 'nVertLevels')
+        normalVelocity = normalVelocity.expand_dims(dim='Time', axis=0)
+
+        ds['temperature'] = temperature
+        ds['salinity'] = salinity * xarray.ones_like(temperature)
+        ds['normalVelocity'] = normalVelocity
+
+        write_netcdf(ds, 'ocean.nc')

compass/ocean/tests/internal_wave/internal_wave.cfg

@@ -0,0 +1,56 @@
+# Options related to the vertical grid
+[vertical_grid]
+
+# the type of vertical grid
+grid_type = uniform
+
+# Number of vertical levels
+vert_levels = 20
+
+# Depth of the bottom of the ocean
+bottom_depth = 500.0
+
+# The type of vertical coordinate (e.g. z-level, z-star)
+coord_type = z-level
+
+# 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
+
+# config options for internal wave test cases
+[internal_wave]
+
+# the number of grid cells in x and y
+nx = 4
+ny = 50
+
+# the size of grid cells (m)
+dc = 5000.0
+
+# Logical flag that determines if locations of features are defined by distance
+# or fractions. False means fractions.
+use_distances = False
+
+# Temperature of the surface in the northern half of the domain.
+surface_temperature = 20.1
+
+# Temperature of the bottom in the northern half of the domain.
+bottom_temperature = 10.1
+
+# Difference in the temperature field between top and bottom
+temperature_difference = 2.0
+
+# Fraction of domain in Y direction the temperature gradient should be linear
+# over.
+amplitude_width_frac = 0.33
+
+# Width of the temperature gradient around the center sin wave. Default value
+# is relative to a 500km domain in Y.
+amplitude_width_dist = 50e3
+
+# Salinity of the water in the entire domain.
+salinity = 35.0
+
+# isopycnal displacement = 125.0 Isopycnal configuration not implemented

compass/ocean/tests/internal_wave/namelist.forward

@@ -0,0 +1,10 @@
+config_dt = '00:05:00'
+config_btr_dt = '00:00:15'
+config_time_integrator = 'split_explicit'
+config_run_duration = '0000_00:15:00'
+config_use_cvmix = .true.
+config_pio_stride = 4
+config_use_mom_del2 = .true.
+config_mom_del2 = 10.0
+config_implicit_bottom_drag_coeff = 1.0e-2
+config_use_cvmix_convection = .true.

compass/ocean/tests/internal_wave/namelist.init

@@ -0,0 +1,3 @@
+config_init_configuration = 'internal_waves'
+config_vert_levels = 1
+config_write_cull_cell_mask = .true.