From dadeb44c2512a2a71ff6c6c97866d4d5ffad1302 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 22 Jun 2021 17:36:11 -0500 Subject: [PATCH 01/16] Add inactive top cell option to init step --- compass/ocean/tests/global_ocean/__init__.py | 15 +++++++-- .../ocean/tests/global_ocean/init/__init__.py | 5 +-- .../tests/global_ocean/init/initial_state.py | 33 ++++++++++++++++++- .../tests/global_ocean/init/streams.init | 1 + 4 files changed, 48 insertions(+), 6 deletions(-) diff --git a/compass/ocean/tests/global_ocean/__init__.py b/compass/ocean/tests/global_ocean/__init__.py index b68a3f6d58..7b2c83a539 100644 --- a/compass/ocean/tests/global_ocean/__init__.py +++ b/compass/ocean/tests/global_ocean/__init__.py @@ -57,6 +57,13 @@ def __init__(self, mpas_core): include_regression=True, include_phc=True, include_en4_1900=True) + self._add_tests(mesh_names=['QU240', 'Icos240', 'QUwISC240'], + DynamicAdjustment=QU240DynamicAdjustment, + high_res_topography=False, + include_rk4=True, + include_regression=True, + include_phc=True, + with_inactive_top_cells=True) # for other meshes, we do fewer tests self._add_tests(mesh_names=['QU', 'Icos', 'QUwISC', 'IcoswISC'], @@ -88,7 +95,7 @@ def __init__(self, mpas_core): def _add_tests(self, mesh_names, DynamicAdjustment, high_res_topography=True, include_rk4=False, include_regression=False, include_phc=False, - include_en4_1900=False): + include_en4_1900=False, with_inactive_top_cells=False): """ Add test cases for the given mesh(es) """ default_ic = 'WOA23' @@ -100,7 +107,8 @@ def _add_tests(self, mesh_names, DynamicAdjustment, self.add_test_case(mesh_test) init_test = Init(test_group=self, mesh=mesh_test, - initial_condition=default_ic) + initial_condition=default_ic, + with_inactive_top_cells=with_inactive_top_cells) self.add_test_case(init_test) time_integrator = default_time_int @@ -173,7 +181,8 @@ def _add_tests(self, mesh_names, DynamicAdjustment, # additional initial conditions (if any) time_integrator = default_time_int init_test = Init(test_group=self, mesh=mesh_test, - initial_condition=initial_condition) + initial_condition=initial_condition, + with_inactive_top_cells=with_inactive_top_cells) self.add_test_case(init_test) self.add_test_case( diff --git a/compass/ocean/tests/global_ocean/init/__init__.py b/compass/ocean/tests/global_ocean/init/__init__.py index c59841abeb..779e5c6f4e 100644 --- a/compass/ocean/tests/global_ocean/init/__init__.py +++ b/compass/ocean/tests/global_ocean/init/__init__.py @@ -25,7 +25,7 @@ class Init(TestCase): The subdirectory within the test group for all test cases with this initial condition """ - def __init__(self, test_group, mesh, initial_condition): + def __init__(self, test_group, mesh, initial_condition, with_inactive_top_cells): """ Create the test case @@ -53,7 +53,8 @@ def __init__(self, test_group, mesh, initial_condition): self.add_step( InitialState( test_case=self, mesh=mesh, - initial_condition=initial_condition)) + initial_condition=initial_condition, + with_inactive_top_cells=with_inactive_top_cells)) if mesh.with_ice_shelf_cavities: self.add_step(RemapIceShelfMelt(test_case=self, mesh=mesh)) diff --git a/compass/ocean/tests/global_ocean/init/initial_state.py b/compass/ocean/tests/global_ocean/init/initial_state.py index d9ea96eedf..efdb212cf7 100644 --- a/compass/ocean/tests/global_ocean/init/initial_state.py +++ b/compass/ocean/tests/global_ocean/init/initial_state.py @@ -1,5 +1,6 @@ import os from importlib.resources import contents +import xarray from compass.model import run_model from compass.ocean.plot import plot_initial_state, plot_vertical_grid @@ -9,6 +10,8 @@ from compass.ocean.vertical.grid_1d import generate_1d_grid, write_1d_grid from compass.step import Step +from mpas_tools.io import write_netcdf + class InitialState(Step): """ @@ -23,7 +26,8 @@ class InitialState(Step): initial_condition : {'WOA23', 'PHC', 'EN4_1900'} The initial condition dataset to use """ - def __init__(self, test_case, mesh, initial_condition): + def __init__(self, test_case, mesh, initial_condition, + with_inactive_top_cells): """ Create the step @@ -44,6 +48,7 @@ def __init__(self, test_case, mesh, initial_condition): super().__init__(test_case=test_case, name='initial_state') self.mesh = mesh self.initial_condition = initial_condition + self.with_inactive_top_cells = with_inactive_top_cells package = 'compass.ocean.tests.global_ocean.init' @@ -170,6 +175,32 @@ def run(self): update_pio = config.getboolean('global_ocean', 'init_update_pio') run_model(self, update_pio=update_pio) + + if self.with_inactive_top_cells: + + logger = self.logger + logger.info(" * Updating min,maxLevelCell for inactive top cells") + + in_filename = 'initial_state.nc' + out_filename = in_filename + + with xarray.open_dataset(in_filename) as ds: + + # keep the data set with Time for output + ds_out = ds + + ds = ds.isel(Time=0) + + if ('maxLevelCell' in ds) and ('minLevelCell' in ds): + minLevelCell = ds.minLevelCell+1 + maxLevelCell = ds.maxLevelCell+1 + ds_out['minLevelCell'] = minLevelCell + ds_out['maxLevelCell'] = maxLevelCell + else: + logger.info(" - Streams missing for inactive top cells") + + write_netcdf(ds_out, out_filename) + logger.info(" - Complete") add_mesh_and_init_metadata(self.outputs, config, init_filename='initial_state.nc') diff --git a/compass/ocean/tests/global_ocean/init/streams.init b/compass/ocean/tests/global_ocean/init/streams.init index c105b6041e..aff8c8bc29 100644 --- a/compass/ocean/tests/global_ocean/init/streams.init +++ b/compass/ocean/tests/global_ocean/init/streams.init @@ -23,6 +23,7 @@ + From 1da84c90e3cb179d0cfd774982b75eb6a162ff98 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 23 Jun 2021 09:15:52 -0500 Subject: [PATCH 02/16] Make inactive top cells false by default --- compass/ocean/tests/global_ocean/init/__init__.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/compass/ocean/tests/global_ocean/init/__init__.py b/compass/ocean/tests/global_ocean/init/__init__.py index 779e5c6f4e..f8fe13c106 100644 --- a/compass/ocean/tests/global_ocean/init/__init__.py +++ b/compass/ocean/tests/global_ocean/init/__init__.py @@ -25,7 +25,8 @@ class Init(TestCase): The subdirectory within the test group for all test cases with this initial condition """ - def __init__(self, test_group, mesh, initial_condition, with_inactive_top_cells): + def __init__(self, test_group, mesh, initial_condition, + with_inactive_top_cells=False): """ Create the test case @@ -43,12 +44,16 @@ def __init__(self, test_group, mesh, initial_condition, with_inactive_top_cells) name = 'init' mesh_name = mesh.mesh_name ic_dir = initial_condition - self.init_subdir = os.path.join(mesh_name, ic_dir) + init_subdir = os.path.join(mesh_name, ic_dir) + if with_inactive_top_cells: + init_subdir = os.path.join(init_subdir, inactive_top) + self.init_subdir = init_subdir subdir = os.path.join(self.init_subdir, name) super().__init__(test_group=test_group, name=name, subdir=subdir) self.mesh = mesh self.initial_condition = initial_condition + self.with_inactive_top_cells = with_inactive_top_cells self.add_step( InitialState( From 83233decb91174363e89503c8cb29d5e65a0acac Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 23 Jun 2021 09:28:51 -0500 Subject: [PATCH 03/16] Increase nVertLevels for inactive top layers --- compass/ocean/tests/global_ocean/init/initial_state.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/compass/ocean/tests/global_ocean/init/initial_state.py b/compass/ocean/tests/global_ocean/init/initial_state.py index efdb212cf7..d27d3191bb 100644 --- a/compass/ocean/tests/global_ocean/init/initial_state.py +++ b/compass/ocean/tests/global_ocean/init/initial_state.py @@ -167,6 +167,13 @@ def run(self): Run this step of the testcase """ config = self.config + if self.with_inactive_top_cells: + # Since we start at minLevelCell = 2, we need to increase the + # number of vertical levels in the cfg file to end up with the + # intended number in the initial state + vert_levels = config.getint('vertical_grid', 'vert_levels') + config.set('vertical_grid', 'vert_levels', f'{vert_levels + 1}', + comment='the number of vertical levels + 1') interfaces = generate_1d_grid(config=config) write_1d_grid(interfaces=interfaces, out_filename='vertical_grid.nc') From c059ee5db6fe137f200c1b70c076d97254946d06 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 23 Jun 2021 13:20:14 -0500 Subject: [PATCH 04/16] Update initial_state with minLevelCell --- compass/ocean/tests/global_ocean/init/initial_state.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/ocean/tests/global_ocean/init/initial_state.py b/compass/ocean/tests/global_ocean/init/initial_state.py index d27d3191bb..ce842f6e2f 100644 --- a/compass/ocean/tests/global_ocean/init/initial_state.py +++ b/compass/ocean/tests/global_ocean/init/initial_state.py @@ -190,7 +190,7 @@ def run(self): in_filename = 'initial_state.nc' out_filename = in_filename - + with xarray.open_dataset(in_filename) as ds: # keep the data set with Time for output From 5d7e757f7114529904ed65134a79dff91ca97328 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 23 Jun 2021 16:21:50 -0500 Subject: [PATCH 05/16] Add config option for number of inactive layers --- .../tests/global_ocean/init/initial_state.py | 3 +- compass/ocean/vertical/grid_1d.py | 290 ++++++++++++++++++ 2 files changed, 292 insertions(+), 1 deletion(-) create mode 100644 compass/ocean/vertical/grid_1d.py diff --git a/compass/ocean/tests/global_ocean/init/initial_state.py b/compass/ocean/tests/global_ocean/init/initial_state.py index ce842f6e2f..4662c03b6b 100644 --- a/compass/ocean/tests/global_ocean/init/initial_state.py +++ b/compass/ocean/tests/global_ocean/init/initial_state.py @@ -174,6 +174,7 @@ def run(self): vert_levels = config.getint('vertical_grid', 'vert_levels') config.set('vertical_grid', 'vert_levels', f'{vert_levels + 1}', comment='the number of vertical levels + 1') + config.set('vertical_grid', 'inactive_top_cells', '1') interfaces = generate_1d_grid(config=config) write_1d_grid(interfaces=interfaces, out_filename='vertical_grid.nc') @@ -190,7 +191,7 @@ def run(self): in_filename = 'initial_state.nc' out_filename = in_filename - + with xarray.open_dataset(in_filename) as ds: # keep the data set with Time for output diff --git a/compass/ocean/vertical/grid_1d.py b/compass/ocean/vertical/grid_1d.py new file mode 100644 index 0000000000..e71d35f2a7 --- /dev/null +++ b/compass/ocean/vertical/grid_1d.py @@ -0,0 +1,290 @@ +import numpy +from importlib import resources +import json +from netCDF4 import Dataset +import numpy as np +from scipy.optimize import root_scalar + + +def generate_1d_grid(config): + """ + Generate a vertical grid for a test case, using the config options in the + ``vertical_grid`` section + + Parameters + ---------- + config : configparser.ConfigParser + Configuration options with parameters used to construct the vertical + grid + + Returns + ------- + interfaces : numpy.ndarray + A 1D array of positive depths for layer interfaces in meters + """ + offset = 0 + if config.has_option('vertical_grid', 'inactive_top_cells'): + offset = config.getint('vertical_grid', 'inactive_top_cells') + + section = config['vertical_grid'] + grid_type = section.get('grid_type') + if grid_type == 'uniform': + vert_levels = section.getint('vert_levels') + interfaces = _generate_uniform(vert_levels - offset) + elif grid_type == 'tanh_dz': + vert_levels = section.getint('vert_levels') + min_layer_thickness = section.getfloat('min_layer_thickness') + max_layer_thickness = section.getfloat('max_layer_thickness') + bottom_depth = section.getfloat('bottom_depth') + interfaces = _create_tanh_dz_grid(vert_levels - offset, + bottom_depth, + min_layer_thickness, + max_layer_thickness) + + elif grid_type in ['60layerPHC', '100layerE3SMv1']: + interfaces = _read_json(grid_type) + else: + raise ValueError('Unexpected grid type: {}'.format(grid_type)) + + if config.has_option('vertical_grid', 'bottom_depth') and \ + grid_type != 'tanh_dz': + bottom_depth = section.getfloat('bottom_depth') + # renormalize to the requested range + interfaces = (bottom_depth/interfaces[-1]) * interfaces + + if config.has_option('vertical_grid', 'inactive_top_cells'): + interfaces = np.append(np.zeros((offset)), interfaces) + + return interfaces + + +def write_1d_grid(interfaces, out_filename): + """ + write the vertical grid to a file + + Parameters + ---------- + interfaces : numpy.ndarray + A 1D array of positive depths for layer interfaces in meters + + out_filename : str + MPAS file name for output of vertical grid + """ + + nz = len(interfaces) - 1 + + # open a new netCDF file for writing. + ncfile = Dataset(out_filename, 'w') + # create the depth_t dimension. + ncfile.createDimension('nVertLevels', nz) + + refBottomDepth = ncfile.createVariable( + 'refBottomDepth', np.dtype('float64').char, ('nVertLevels',)) + refMidDepth = ncfile.createVariable( + 'refMidDepth', np.dtype('float64').char, ('nVertLevels',)) + refLayerThickness = ncfile.createVariable( + 'refLayerThickness', np.dtype('float64').char, ('nVertLevels',)) + + botDepth = interfaces[1:] + midDepth = 0.5 * (interfaces[0:-1] + interfaces[1:]) + + refBottomDepth[:] = botDepth + refMidDepth[:] = midDepth + refLayerThickness[:] = interfaces[1:] - interfaces[0:-1] + ncfile.close() + + +def add_1d_grid(config, ds): + """ + Add a 1D vertical grid based on the config options in the ``vertical_grid`` + section to a mesh data set + + The following variables are added to the mesh: + * ``refTopDepth`` - the positive-down depth of the top of each ref. level + * ``refZMid`` - the positive-down depth of the middle of each ref. level + * ``refBottomDepth`` - the positive-down depth of the bottom of each ref. + level + * ``refInterfaces`` - the positive-down depth of the interfaces between + ref. levels (with ``nVertLevels`` + 1 elements). + There is considerable redundancy between these variables but each is + sometimes convenient. + + Parameters + ---------- + config : configparser.ConfigParser + Configuration options with parameters used to construct the vertical + grid + + ds : xarray.Dataset + A data set to add the grid variables to + """ + + interfaces = generate_1d_grid(config=config) + + ds['refTopDepth'] = ('nVertLevels', interfaces[0:-1]) + ds['refZMid'] = ('nVertLevels', -0.5 * (interfaces[1:] + interfaces[0:-1])) + ds['refBottomDepth'] = ('nVertLevels', interfaces[1:]) + ds['refInterfaces'] = ('nVertLevelsP1', interfaces) + + +def _generate_uniform(vert_levels): + """ Generate uniform layer interfaces between 0 and 1 """ + interfaces = numpy.linspace(0., 1., vert_levels+1) + return interfaces + + +def _read_json(grid_type): + """ Read the grid interfaces from a json file """ + + filename = '{}.json'.format(grid_type) + with resources.open_text("compass.ocean.vertical", filename) as data_file: + data = json.load(data_file) + interfaces = numpy.array(data) + + return interfaces + + +def _create_tanh_dz_grid(num_vert_levels, bottom_depth, min_layer_thickness, + max_layer_thickness): + """ + Creates the vertical grid for MPAS-Ocean and writes it to a NetCDF file + + Parameters + ---------- + num_vert_levels : int + Number of vertical levels for the grid + + bottom_depth : float + bottom depth for the chosen vertical coordinate [m] + + min_layer_thickness : float + Target thickness of the first layer [m] + + max_layer_thickness : float + Target maximum thickness in column [m] + + Returns + ------- + interfaces : numpy.ndarray + A 1D array of positive depths for layer interfaces in meters + """ + + nz = num_vert_levels + dz1 = min_layer_thickness + dz2 = max_layer_thickness + + # the bracket here is large enough that it should hopefully encompass any + # reasonable value of delta, the characteristic length scale over which + # dz varies. The args are passed on to the match_bottom function below, + # and the root finder will determine a value of delta (sol.root) such that + # match_bottom is within a tolerance of zero, meaning the bottom of the + # coordinate computed by cumsum_z hits bottom_depth almost exactly + sol = root_scalar(_match_bottom, method='brentq', + bracket=[dz1, 10 * bottom_depth], + args=(nz, dz1, dz2, bottom_depth)) + + delta = sol.root + layerThickness, z = _cumsum_z(delta, nz, dz1, dz2) + interfaces = -z + + return interfaces + + +def _match_bottom(delta, nz, dz1, dz2, bottom_depth): + """ + Compute the difference between the bottom depth computed with the given + parameters and the target ``bottom_depth``, used in the root finding + algorithm to determine which value of ``delta`` to use. + + Parameters + ---------- + delta : float + The characteristic length scale over which dz varies (this parameter + will be optimized to hit a target depth in a target number of layers) + + nz : int + The number of layers + + dz1 : float + The layer thickness at the top of the ocean (z = 0) + + dz2 : float + The layer thickness at z --> -infinity + + bottom_depth: float + depth of the bottom of the ocean that should match the bottom layer + interface. Note: the bottom_depth is positive, whereas the layer + interfaces are negative. + + Returns + ------- + diff : float + The computed bottom depth minus the target ``bottom_depth``. ``diff`` + should be zero when we have found the desired ``delta``. + """ + _, z = _cumsum_z(delta, nz, dz1, dz2) + diff = -bottom_depth - z[-1] + return diff + + +def _cumsum_z(delta, nz, dz1, dz2): + """ + Compute layer interface depths and layer thicknesses over ``nz`` layers + + Parameters + ---------- + delta : float + The characteristic length scale over which dz varies (this parameter + will be optimized to hit a target depth in a target number of layers) + + nz : int + The number of layers + + dz1 : float + The layer thickness at the top of the ocean (z = 0) + + dz2 : float + The layer thickness at z --> -infinity + + Returns + ------- + dz : numpy.ndarray + The layer thicknesses for each layer + + z : numpy.ndarray + The depth (positive up) of each layer interface (``nz + 1`` total + elements) + """ + dz = np.zeros(nz) + z = np.zeros(nz + 1) + for zindex in range(nz): + dz[zindex] = _dz_z(z[zindex], dz1, dz2, delta) + z[zindex + 1] = z[zindex] - dz[zindex] + return dz, z + + +def _dz_z(z, dz1, dz2, delta): + """ + layer thickness as a function of depth + + Parameters + ---------- + z : float + Depth coordinate (positive up) at which to find the layer thickness + + dz1 : float + The layer thickness at the top of the ocean (z = 0) + + dz2 : float + The layer thickness at z --> -infinity + + delta : float + The characteristic length scale over which dz varies (this parameter + will be optimized to hit a target depth in a target numer of layers) + + Returns + ------- + dz : float + The layer thickness + """ + return (dz2 - dz1) * np.tanh(-z * np.pi / delta) + dz1 From c7a8094bfb4854907dc8de1cdb8098986078941b Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 24 Jun 2021 15:40:45 -0500 Subject: [PATCH 06/16] Remove maxLevelCell assignment --- compass/ocean/tests/global_ocean/init/initial_state.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/compass/ocean/tests/global_ocean/init/initial_state.py b/compass/ocean/tests/global_ocean/init/initial_state.py index 4662c03b6b..515f35bbf5 100644 --- a/compass/ocean/tests/global_ocean/init/initial_state.py +++ b/compass/ocean/tests/global_ocean/init/initial_state.py @@ -187,7 +187,7 @@ def run(self): if self.with_inactive_top_cells: logger = self.logger - logger.info(" * Updating min,maxLevelCell for inactive top cells") + logger.info(" * Updating minLevelCell for inactive top cells") in_filename = 'initial_state.nc' out_filename = in_filename @@ -199,11 +199,9 @@ def run(self): ds = ds.isel(Time=0) - if ('maxLevelCell' in ds) and ('minLevelCell' in ds): + if ('minLevelCell' in ds): minLevelCell = ds.minLevelCell+1 - maxLevelCell = ds.maxLevelCell+1 ds_out['minLevelCell'] = minLevelCell - ds_out['maxLevelCell'] = maxLevelCell else: logger.info(" - Streams missing for inactive top cells") From 95e0ca2a7eda2ba38347a67d2a6cb9d6133abf28 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 9 Jul 2021 15:11:37 -0500 Subject: [PATCH 07/16] Update documentation for inactive top cell test --- docs/developers_guide/ocean/test_groups/global_ocean.rst | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/developers_guide/ocean/test_groups/global_ocean.rst b/docs/developers_guide/ocean/test_groups/global_ocean.rst index 91e28e8af4..5cc2646880 100644 --- a/docs/developers_guide/ocean/test_groups/global_ocean.rst +++ b/docs/developers_guide/ocean/test_groups/global_ocean.rst @@ -1013,6 +1013,11 @@ If a baseline is provided, prognostic variables and ice-shelf melt fluxes (if ice-shelf cavities are included in the mesh) are compared with a baseline, and the ``time integration`` timer is compared with that of the baseline. +This test case can also serve as a verification of the inactive top cell +implementation. If ``inactive_top_cells`` is given in the config file with +values greater than zero, the vertical domain is shifted downward by the given +number of layers. + .. _dev_ocean_global_ocean_restart_test: restart_test test case From c75a01b57ec3124d4dd881107e5973f5ba4d16cb Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 19 Jul 2021 13:49:54 -0600 Subject: [PATCH 08/16] Add validation with cropped output --- compass/ocean/inactive_top_cells.py | 28 +++++++++++++++++++ .../ocean/tests/global_ocean/init/__init__.py | 8 ++++++ .../tests/global_ocean/init/initial_state.py | 15 ++++++++-- .../global_ocean/performance_test/__init__.py | 14 ++++++++++ 4 files changed, 63 insertions(+), 2 deletions(-) create mode 100644 compass/ocean/inactive_top_cells.py diff --git a/compass/ocean/inactive_top_cells.py b/compass/ocean/inactive_top_cells.py new file mode 100644 index 0000000000..fff8a311e1 --- /dev/null +++ b/compass/ocean/inactive_top_cells.py @@ -0,0 +1,28 @@ +import os +import xarray + +from mpas_tools.io import write_netcdf + + +def remove_inactive_top_cells_output(in_filename, inactive_top_cells=1): + """ + Remove inactive top cells from the output netCDF file + + Parameters + ---------- + in_filename : str + Filename for the netCDF file to be cropped + + inactive_top_cells : str, optional + The number of inactive top cell layers. It should be equal to + ``config.inactive_top_cells``. + """ + if not os.path.exists(in_filename): + raise OSError('File {} does not exist.'.format(in_filename)) + + out_filename = in_filename.split('.')[0] + '_crop.nc' + + with xarray.open_dataset(in_filename) as ds_in: + ds_out = ds_in.isel(nVertLevels=slice(1, None)) + + write_netcdf(ds_out, out_filename) diff --git a/compass/ocean/tests/global_ocean/init/__init__.py b/compass/ocean/tests/global_ocean/init/__init__.py index f8fe13c106..f26063ca49 100644 --- a/compass/ocean/tests/global_ocean/init/__init__.py +++ b/compass/ocean/tests/global_ocean/init/__init__.py @@ -99,6 +99,14 @@ def validate(self): compare_variables(test_case=self, variables=variables, filename1='initial_state/initial_state.nc') + if self.with_inactive_top_cells: + variables = [ + 'temperature', 'salinity', 'layerThickness', 'normalVelocity'] + compare_variables(test_case=self, variables=variables, + filename1='initial_state/initial_state_crop.nc', + filename2='initial_state/initial_state_comp.nc', + quiet=False, check_outputs=False) + if self.mesh.with_ice_shelf_cavities: variables = ['ssh', 'landIcePressure'] compare_variables(test_case=self, variables=variables, diff --git a/compass/ocean/tests/global_ocean/init/initial_state.py b/compass/ocean/tests/global_ocean/init/initial_state.py index 515f35bbf5..94854d3e0d 100644 --- a/compass/ocean/tests/global_ocean/init/initial_state.py +++ b/compass/ocean/tests/global_ocean/init/initial_state.py @@ -7,6 +7,7 @@ from compass.ocean.tests.global_ocean.metadata import ( add_mesh_and_init_metadata, ) +from compass.ocean.inactive_top_cells import remove_inactive_top_cells_output from compass.ocean.vertical.grid_1d import generate_1d_grid, write_1d_grid from compass.step import Step @@ -139,6 +140,9 @@ def __init__(self, test_case, mesh, initial_condition, 'graph.info']: self.add_output_file(filename=file) + if with_inactive_top_cells: + self.add_output_file(filename='initial_state_crop.nc') + def setup(self): """ Get resources at setup from config options @@ -175,6 +179,7 @@ def run(self): config.set('vertical_grid', 'vert_levels', f'{vert_levels + 1}', comment='the number of vertical levels + 1') config.set('vertical_grid', 'inactive_top_cells', '1') + logger = self.logger interfaces = generate_1d_grid(config=config) write_1d_grid(interfaces=interfaces, out_filename='vertical_grid.nc') @@ -186,7 +191,6 @@ def run(self): if self.with_inactive_top_cells: - logger = self.logger logger.info(" * Updating minLevelCell for inactive top cells") in_filename = 'initial_state.nc' @@ -200,11 +204,18 @@ def run(self): ds = ds.isel(Time=0) if ('minLevelCell' in ds): - minLevelCell = ds.minLevelCell+1 + if config.has_option('vertical_grid', + 'inactive_top_cells'): + offset = config.getint('vertical_grid', + 'inactive_top_cells') + minLevelCell = ds.minLevelCell+offset ds_out['minLevelCell'] = minLevelCell else: logger.info(" - Streams missing for inactive top cells") + remove_inactive_top_cells_output(in_filename, + inactive_top_cells=offset) + write_netcdf(ds_out, out_filename) logger.info(" - Complete") diff --git a/compass/ocean/tests/global_ocean/performance_test/__init__.py b/compass/ocean/tests/global_ocean/performance_test/__init__.py index 0487f4f9f3..a523c14aae 100644 --- a/compass/ocean/tests/global_ocean/performance_test/__init__.py +++ b/compass/ocean/tests/global_ocean/performance_test/__init__.py @@ -1,8 +1,10 @@ +import os from compass.ocean.tests.global_ocean.forward import ( ForwardStep, ForwardTestCase, ) from compass.validate import compare_timers, compare_variables +from compass.ocean.inactive_top_cells import remove_inactive_top_cells_output class PerformanceTest(ForwardTestCase): @@ -71,6 +73,18 @@ def validate(self): compare_variables(test_case=self, variables=variables, filename1=f'{step_subdir}/output.nc') + if self.config.has_option('vertical_grid', 'inactive_top_cells'): + offset = self.config.getint('vertical_grid', 'inactive_top_cells') + if offset > 0: + remove_inactive_top_cells_output('forward/output.nc', + inactive_top_cells=offset) + filename2=None + if os.path.exists('forward/output_comp.nc'): + filename2='forward/output_comp.nc' + compare_variables(test_case=self, variables=variables, + filename1='forward/output_crop.nc', + filename2=filename2) + if self.mesh.with_ice_shelf_cavities: variables = [ 'ssh', 'landIcePressure', 'landIceDraft', From 98061df607c96046bb9ecea3699819c70df27cdc Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 1 Aug 2021 11:59:54 +0200 Subject: [PATCH 09/16] Add mesh and out filename when cropping inactive top Use `minLevelCell` instead of an argument to determine the number of inactive top levels. --- compass/ocean/inactive_top_cells.py | 31 ++++++++++++++----- .../ocean/tests/global_ocean/init/__init__.py | 24 ++++++++++---- .../global_ocean/performance_test/__init__.py | 20 +++++++----- 3 files changed, 54 insertions(+), 21 deletions(-) diff --git a/compass/ocean/inactive_top_cells.py b/compass/ocean/inactive_top_cells.py index fff8a311e1..c28599e68e 100644 --- a/compass/ocean/inactive_top_cells.py +++ b/compass/ocean/inactive_top_cells.py @@ -4,7 +4,8 @@ from mpas_tools.io import write_netcdf -def remove_inactive_top_cells_output(in_filename, inactive_top_cells=1): +def remove_inactive_top_cells_output(in_filename, out_filename=None, + mesh_filename=None): """ Remove inactive top cells from the output netCDF file @@ -13,16 +14,32 @@ def remove_inactive_top_cells_output(in_filename, inactive_top_cells=1): in_filename : str Filename for the netCDF file to be cropped - inactive_top_cells : str, optional - The number of inactive top cell layers. It should be equal to - ``config.inactive_top_cells``. + out_filename : str, optional + Filename for the netCDF file after cropping. Tbe default name is the + original file name with ``_crop`` appended before the extension + + mesh_filename : str, optional + Filename for an MPAS mesh if not included in the file to be cropped + """ if not os.path.exists(in_filename): - raise OSError('File {} does not exist.'.format(in_filename)) + raise OSError(f'File {in_filename} does not exist.') - out_filename = in_filename.split('.')[0] + '_crop.nc' + if out_filename is None: + basename, ext = os.path.splitext(in_filename) + out_filename = f'{basename}_crop{ext}' with xarray.open_dataset(in_filename) as ds_in: - ds_out = ds_in.isel(nVertLevels=slice(1, None)) + if mesh_filename is None: + ds_mesh = ds_in + else: + ds_mesh = xarray.open_dataset(mesh_filename) + minLevelCell = ds_mesh.minLevelCell + minval = minLevelCell.min().values + maxval = minLevelCell.max().values + if minval != maxval: + raise ValueError('Expected minLevelCell to have a constant ' + 'value for inactive top cell tests') + ds_out = ds_in.isel(nVertLevels=slice(minval-1, None)) write_netcdf(ds_out, out_filename) diff --git a/compass/ocean/tests/global_ocean/init/__init__.py b/compass/ocean/tests/global_ocean/init/__init__.py index f26063ca49..233dfa1c5f 100644 --- a/compass/ocean/tests/global_ocean/init/__init__.py +++ b/compass/ocean/tests/global_ocean/init/__init__.py @@ -100,12 +100,24 @@ def validate(self): filename1='initial_state/initial_state.nc') if self.with_inactive_top_cells: - variables = [ - 'temperature', 'salinity', 'layerThickness', 'normalVelocity'] - compare_variables(test_case=self, variables=variables, - filename1='initial_state/initial_state_crop.nc', - filename2='initial_state/initial_state_comp.nc', - quiet=False, check_outputs=False) + # construct the work directory for the other test + filename2 = os.path.join(self.base_work_dir, self.mpas_core.name, + self.test_group.name, + self.inactive_top_comp_subdir, + 'init/initial_state/initial_state.nc') + if os.path.exists(filename2): + variables = ['temperature', 'salinity', 'layerThickness', + 'normalVelocity'] + compare_variables(test_case=self, variables=variables, + filename1='initial_state/initial_state_crop.nc' + filename2=filename2, + quiet=False, check_outputs=False, + skip_if_step_not_run=False) + + else: + self.logger.warn('The version of "init" without inactive top ' + 'cells was not run. Skipping\n' + 'validation.') if self.mesh.with_ice_shelf_cavities: variables = ['ssh', 'landIcePressure'] diff --git a/compass/ocean/tests/global_ocean/performance_test/__init__.py b/compass/ocean/tests/global_ocean/performance_test/__init__.py index a523c14aae..59fcfd81d4 100644 --- a/compass/ocean/tests/global_ocean/performance_test/__init__.py +++ b/compass/ocean/tests/global_ocean/performance_test/__init__.py @@ -73,17 +73,21 @@ def validate(self): compare_variables(test_case=self, variables=variables, filename1=f'{step_subdir}/output.nc') - if self.config.has_option('vertical_grid', 'inactive_top_cells'): - offset = self.config.getint('vertical_grid', 'inactive_top_cells') - if offset > 0: - remove_inactive_top_cells_output('forward/output.nc', - inactive_top_cells=offset) - filename2=None - if os.path.exists('forward/output_comp.nc'): - filename2='forward/output_comp.nc' + if self.init.with_inactive_top_cells: + # construct the work directory for the other test + subdir = get_forward_subdir(self.init.inactive_top_comp_subdir, + self.time_integrator, self.name) + filename2 = os.path.join(self.base_work_dir, self.mpas_core.name, + self.test_group.name, subdir, + 'forward/output.nc') + if os.path.exists(filename2): compare_variables(test_case=self, variables=variables, filename1='forward/output_crop.nc', filename2=filename2) + else: + self.logger.warn('The version of "performance_test" without ' + 'inactive top cells was not run.\n' + 'Skipping validation.') if self.mesh.with_ice_shelf_cavities: variables = [ From ba3b04a364e1f494f7adee9d0271b263b24676f9 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 1 Aug 2021 12:02:30 +0200 Subject: [PATCH 10/16] Update the init test case Store the `inactive_top_comp_subdir` so we know what files to compare with during validation. Use `ds.load()` when modifying `minLevelCell` in the initial state to prevent issues with reading from and writing to the same file. Crop inactive top cells after writing out the modified initial state (because we now use the modified `minLevelCell`). Set the default to `inactive_top_cells = 0` in a new `init.cfg` instead of the config file for the QU240 mesh. --- .../ocean/tests/global_ocean/init/__init__.py | 8 +++++ .../ocean/tests/global_ocean/init/init.cfg | 5 +++ .../tests/global_ocean/init/initial_state.py | 36 +++++++++++-------- 3 files changed, 34 insertions(+), 15 deletions(-) create mode 100644 compass/ocean/tests/global_ocean/init/init.cfg diff --git a/compass/ocean/tests/global_ocean/init/__init__.py b/compass/ocean/tests/global_ocean/init/__init__.py index 233dfa1c5f..31309a1933 100644 --- a/compass/ocean/tests/global_ocean/init/__init__.py +++ b/compass/ocean/tests/global_ocean/init/__init__.py @@ -24,6 +24,11 @@ class Init(TestCase): init_subdir : str The subdirectory within the test group for all test cases with this initial condition + + inactive_top_comp_subdir : str + If ``with_inactive_top_cells == True``, the subdirectory equivalent to + ``init_subdir`` for test cases without inactive top cells for use for + validation between tests with and without inactive top cells """ def __init__(self, test_group, mesh, initial_condition, with_inactive_top_cells=False): @@ -46,6 +51,9 @@ def __init__(self, test_group, mesh, initial_condition, ic_dir = initial_condition init_subdir = os.path.join(mesh_name, ic_dir) if with_inactive_top_cells: + # Add the name of the subdir without inactive top cells whether or + # not is has or will be run + self.inactive_top_comp_subdir = init_subdir init_subdir = os.path.join(init_subdir, inactive_top) self.init_subdir = init_subdir subdir = os.path.join(self.init_subdir, name) diff --git a/compass/ocean/tests/global_ocean/init/init.cfg b/compass/ocean/tests/global_ocean/init/init.cfg new file mode 100644 index 0000000000..3688d8510d --- /dev/null +++ b/compass/ocean/tests/global_ocean/init/init.cfg @@ -0,0 +1,5 @@ +# Options related to the vertical grid +[vertical_grid] + +# no inactive top cells by default +inactive_top_cells = 0 diff --git a/compass/ocean/tests/global_ocean/init/initial_state.py b/compass/ocean/tests/global_ocean/init/initial_state.py index 94854d3e0d..4277e0fa82 100644 --- a/compass/ocean/tests/global_ocean/init/initial_state.py +++ b/compass/ocean/tests/global_ocean/init/initial_state.py @@ -1,18 +1,18 @@ import os from importlib.resources import contents + import xarray +from mpas_tools.io import write_netcdf from compass.model import run_model +from compass.ocean.inactive_top_cells import remove_inactive_top_cells_output from compass.ocean.plot import plot_initial_state, plot_vertical_grid from compass.ocean.tests.global_ocean.metadata import ( add_mesh_and_init_metadata, ) -from compass.ocean.inactive_top_cells import remove_inactive_top_cells_output from compass.ocean.vertical.grid_1d import generate_1d_grid, write_1d_grid from compass.step import Step -from mpas_tools.io import write_netcdf - class InitialState(Step): """ @@ -27,7 +27,7 @@ class InitialState(Step): initial_condition : {'WOA23', 'PHC', 'EN4_1900'} The initial condition dataset to use """ - def __init__(self, test_case, mesh, initial_condition, + def __init__(self, test_case, mesh, initial_condition, with_inactive_top_cells): """ Create the step @@ -188,7 +188,7 @@ def run(self): update_pio = config.getboolean('global_ocean', 'init_update_pio') run_model(self, update_pio=update_pio) - + if self.with_inactive_top_cells: logger.info(" * Updating minLevelCell for inactive top cells") @@ -197,26 +197,32 @@ def run(self): out_filename = in_filename with xarray.open_dataset(in_filename) as ds: + ds.load() # keep the data set with Time for output ds_out = ds ds = ds.isel(Time=0) - if ('minLevelCell' in ds): - if config.has_option('vertical_grid', - 'inactive_top_cells'): - offset = config.getint('vertical_grid', - 'inactive_top_cells') - minLevelCell = ds.minLevelCell+offset - ds_out['minLevelCell'] = minLevelCell + if config.has_option('vertical_grid', 'inactive_top_cells'): + offset = config.getint('vertical_grid', + 'inactive_top_cells') else: - logger.info(" - Streams missing for inactive top cells") + offset = 0 - remove_inactive_top_cells_output(in_filename, - inactive_top_cells=offset) + if 'minLevelCell' in ds: + minLevelCell = ds.minLevelCell + offset + ds_out['minLevelCell'] = minLevelCell + else: + logger.info(" - Variable minLevelCell, needed for " + "inactive top cells, is missing from the " + "initial condition") write_netcdf(ds_out, out_filename) + + remove_inactive_top_cells_output( + in_filename=in_filename, out_filename='initial_state_crop.nc') + logger.info(" - Complete") add_mesh_and_init_metadata(self.outputs, config, From 0bcab9a54bf8c7bab515b62f6c753a3743f8a209 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 1 Aug 2021 12:08:12 +0200 Subject: [PATCH 11/16] Update performance test Move cropping of inactive top cells to ForwardStep so it is explicitly included as an output of that step. Give `init.nc` as the mesh file for cropping. For validation, point to the output file from the version of the performance test without inactive top cells for `filename2`. --- compass/ocean/tests/global_ocean/forward.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/compass/ocean/tests/global_ocean/forward.py b/compass/ocean/tests/global_ocean/forward.py index bcb30eb0f3..2b11836c45 100644 --- a/compass/ocean/tests/global_ocean/forward.py +++ b/compass/ocean/tests/global_ocean/forward.py @@ -3,6 +3,7 @@ from importlib.resources import contents from compass.model import run_model +from compass.ocean.inactive_top_cells import remove_inactive_top_cells_output from compass.ocean.tests.global_ocean.metadata import ( add_mesh_and_init_metadata, ) @@ -155,6 +156,9 @@ def __init__(self, test_case, mesh, time_integrator, init=None, filename='graph.info', work_dir_target=f'{init.path}/initial_state/graph.info') + if init.with_inactive_top_cells: + self.add_output_file(filename='output_crop.nc') + self.add_model_as_input() def setup(self): @@ -214,6 +218,12 @@ def run(self): update_pio = self.config.getboolean('global_ocean', 'forward_update_pio') run_model(self, update_pio=update_pio) + + if self.init.with_inactive_top_cells: + remove_inactive_top_cells_output(in_filename='output.nc', + out_filename='output_crop.nc', + mesh_filename='init.nc') + add_mesh_and_init_metadata(self.outputs, self.config, init_filename='init.nc') From 2a15e8fbd8aa5ee479a7a0cc34d7e1449dd8b426 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 1 Aug 2021 12:11:15 +0200 Subject: [PATCH 12/16] Add a test suite for inactive top cells --- compass/ocean/suites/inactive_top.txt | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 compass/ocean/suites/inactive_top.txt diff --git a/compass/ocean/suites/inactive_top.txt b/compass/ocean/suites/inactive_top.txt new file mode 100644 index 0000000000..2f8b40d60c --- /dev/null +++ b/compass/ocean/suites/inactive_top.txt @@ -0,0 +1,7 @@ +ocean/global_ocean/QU240/mesh +ocean/global_ocean/QU240/WOA23/init +ocean/global_ocean/QU240/WOA23/init_inactive_top +ocean/global_ocean/QU240/WOA23/performance_test +ocean/global_ocean/QU240/WOA23/performance_test_inactive_top +ocean/global_ocean/QU240/WOA23/RK4/performance_test +ocean/global_ocean/QU240/WOA23/RK4/performance_test_inactive_top From 2a41fe5ac815f9bf05ff4ec691f6d53e430fe167 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 25 Aug 2023 14:09:00 -0500 Subject: [PATCH 13/16] Add inactive top cells as tests, not test group --- compass/ocean/tests/global_ocean/__init__.py | 28 +- compass/ocean/vertical/grid_1d.py | 290 ------------------- 2 files changed, 18 insertions(+), 300 deletions(-) delete mode 100644 compass/ocean/vertical/grid_1d.py diff --git a/compass/ocean/tests/global_ocean/__init__.py b/compass/ocean/tests/global_ocean/__init__.py index 7b2c83a539..d81e7b5ac9 100644 --- a/compass/ocean/tests/global_ocean/__init__.py +++ b/compass/ocean/tests/global_ocean/__init__.py @@ -56,13 +56,7 @@ def __init__(self, mpas_core): include_rk4=True, include_regression=True, include_phc=True, - include_en4_1900=True) - self._add_tests(mesh_names=['QU240', 'Icos240', 'QUwISC240'], - DynamicAdjustment=QU240DynamicAdjustment, - high_res_topography=False, - include_rk4=True, - include_regression=True, - include_phc=True, + include_en4_1900=True, with_inactive_top_cells=True) # for other meshes, we do fewer tests @@ -108,7 +102,7 @@ def _add_tests(self, mesh_names, DynamicAdjustment, init_test = Init(test_group=self, mesh=mesh_test, initial_condition=default_ic, - with_inactive_top_cells=with_inactive_top_cells) + with_inactive_top_cells=False) self.add_test_case(init_test) time_integrator = default_time_int @@ -171,6 +165,21 @@ def _add_tests(self, mesh_names, DynamicAdjustment, test_group=self, mesh=mesh_test, init=init_test, time_integrator=time_integrator)) + if with_inactive_top_cells: + init_test = Init(test_group=self, mesh=mesh_test, + initial_condition=default_ic, + with_inactive_top_cells=True) + self.add_test_case(init_test) + self.add_test_case( + PerformanceTest( + test_group=self, mesh=mesh_test, init=init_test, + time_integrator=default_time_int)) + if include_rk4: + self.add_test_case( + PerformanceTest( + test_group=self, mesh=mesh_test, init=init_test, + time_integrator='RK4')) + initial_conditions = [] if include_phc: initial_conditions.append('PHC') @@ -181,8 +190,7 @@ def _add_tests(self, mesh_names, DynamicAdjustment, # additional initial conditions (if any) time_integrator = default_time_int init_test = Init(test_group=self, mesh=mesh_test, - initial_condition=initial_condition, - with_inactive_top_cells=with_inactive_top_cells) + initial_condition=initial_condition) self.add_test_case(init_test) self.add_test_case( diff --git a/compass/ocean/vertical/grid_1d.py b/compass/ocean/vertical/grid_1d.py deleted file mode 100644 index e71d35f2a7..0000000000 --- a/compass/ocean/vertical/grid_1d.py +++ /dev/null @@ -1,290 +0,0 @@ -import numpy -from importlib import resources -import json -from netCDF4 import Dataset -import numpy as np -from scipy.optimize import root_scalar - - -def generate_1d_grid(config): - """ - Generate a vertical grid for a test case, using the config options in the - ``vertical_grid`` section - - Parameters - ---------- - config : configparser.ConfigParser - Configuration options with parameters used to construct the vertical - grid - - Returns - ------- - interfaces : numpy.ndarray - A 1D array of positive depths for layer interfaces in meters - """ - offset = 0 - if config.has_option('vertical_grid', 'inactive_top_cells'): - offset = config.getint('vertical_grid', 'inactive_top_cells') - - section = config['vertical_grid'] - grid_type = section.get('grid_type') - if grid_type == 'uniform': - vert_levels = section.getint('vert_levels') - interfaces = _generate_uniform(vert_levels - offset) - elif grid_type == 'tanh_dz': - vert_levels = section.getint('vert_levels') - min_layer_thickness = section.getfloat('min_layer_thickness') - max_layer_thickness = section.getfloat('max_layer_thickness') - bottom_depth = section.getfloat('bottom_depth') - interfaces = _create_tanh_dz_grid(vert_levels - offset, - bottom_depth, - min_layer_thickness, - max_layer_thickness) - - elif grid_type in ['60layerPHC', '100layerE3SMv1']: - interfaces = _read_json(grid_type) - else: - raise ValueError('Unexpected grid type: {}'.format(grid_type)) - - if config.has_option('vertical_grid', 'bottom_depth') and \ - grid_type != 'tanh_dz': - bottom_depth = section.getfloat('bottom_depth') - # renormalize to the requested range - interfaces = (bottom_depth/interfaces[-1]) * interfaces - - if config.has_option('vertical_grid', 'inactive_top_cells'): - interfaces = np.append(np.zeros((offset)), interfaces) - - return interfaces - - -def write_1d_grid(interfaces, out_filename): - """ - write the vertical grid to a file - - Parameters - ---------- - interfaces : numpy.ndarray - A 1D array of positive depths for layer interfaces in meters - - out_filename : str - MPAS file name for output of vertical grid - """ - - nz = len(interfaces) - 1 - - # open a new netCDF file for writing. - ncfile = Dataset(out_filename, 'w') - # create the depth_t dimension. - ncfile.createDimension('nVertLevels', nz) - - refBottomDepth = ncfile.createVariable( - 'refBottomDepth', np.dtype('float64').char, ('nVertLevels',)) - refMidDepth = ncfile.createVariable( - 'refMidDepth', np.dtype('float64').char, ('nVertLevels',)) - refLayerThickness = ncfile.createVariable( - 'refLayerThickness', np.dtype('float64').char, ('nVertLevels',)) - - botDepth = interfaces[1:] - midDepth = 0.5 * (interfaces[0:-1] + interfaces[1:]) - - refBottomDepth[:] = botDepth - refMidDepth[:] = midDepth - refLayerThickness[:] = interfaces[1:] - interfaces[0:-1] - ncfile.close() - - -def add_1d_grid(config, ds): - """ - Add a 1D vertical grid based on the config options in the ``vertical_grid`` - section to a mesh data set - - The following variables are added to the mesh: - * ``refTopDepth`` - the positive-down depth of the top of each ref. level - * ``refZMid`` - the positive-down depth of the middle of each ref. level - * ``refBottomDepth`` - the positive-down depth of the bottom of each ref. - level - * ``refInterfaces`` - the positive-down depth of the interfaces between - ref. levels (with ``nVertLevels`` + 1 elements). - There is considerable redundancy between these variables but each is - sometimes convenient. - - Parameters - ---------- - config : configparser.ConfigParser - Configuration options with parameters used to construct the vertical - grid - - ds : xarray.Dataset - A data set to add the grid variables to - """ - - interfaces = generate_1d_grid(config=config) - - ds['refTopDepth'] = ('nVertLevels', interfaces[0:-1]) - ds['refZMid'] = ('nVertLevels', -0.5 * (interfaces[1:] + interfaces[0:-1])) - ds['refBottomDepth'] = ('nVertLevels', interfaces[1:]) - ds['refInterfaces'] = ('nVertLevelsP1', interfaces) - - -def _generate_uniform(vert_levels): - """ Generate uniform layer interfaces between 0 and 1 """ - interfaces = numpy.linspace(0., 1., vert_levels+1) - return interfaces - - -def _read_json(grid_type): - """ Read the grid interfaces from a json file """ - - filename = '{}.json'.format(grid_type) - with resources.open_text("compass.ocean.vertical", filename) as data_file: - data = json.load(data_file) - interfaces = numpy.array(data) - - return interfaces - - -def _create_tanh_dz_grid(num_vert_levels, bottom_depth, min_layer_thickness, - max_layer_thickness): - """ - Creates the vertical grid for MPAS-Ocean and writes it to a NetCDF file - - Parameters - ---------- - num_vert_levels : int - Number of vertical levels for the grid - - bottom_depth : float - bottom depth for the chosen vertical coordinate [m] - - min_layer_thickness : float - Target thickness of the first layer [m] - - max_layer_thickness : float - Target maximum thickness in column [m] - - Returns - ------- - interfaces : numpy.ndarray - A 1D array of positive depths for layer interfaces in meters - """ - - nz = num_vert_levels - dz1 = min_layer_thickness - dz2 = max_layer_thickness - - # the bracket here is large enough that it should hopefully encompass any - # reasonable value of delta, the characteristic length scale over which - # dz varies. The args are passed on to the match_bottom function below, - # and the root finder will determine a value of delta (sol.root) such that - # match_bottom is within a tolerance of zero, meaning the bottom of the - # coordinate computed by cumsum_z hits bottom_depth almost exactly - sol = root_scalar(_match_bottom, method='brentq', - bracket=[dz1, 10 * bottom_depth], - args=(nz, dz1, dz2, bottom_depth)) - - delta = sol.root - layerThickness, z = _cumsum_z(delta, nz, dz1, dz2) - interfaces = -z - - return interfaces - - -def _match_bottom(delta, nz, dz1, dz2, bottom_depth): - """ - Compute the difference between the bottom depth computed with the given - parameters and the target ``bottom_depth``, used in the root finding - algorithm to determine which value of ``delta`` to use. - - Parameters - ---------- - delta : float - The characteristic length scale over which dz varies (this parameter - will be optimized to hit a target depth in a target number of layers) - - nz : int - The number of layers - - dz1 : float - The layer thickness at the top of the ocean (z = 0) - - dz2 : float - The layer thickness at z --> -infinity - - bottom_depth: float - depth of the bottom of the ocean that should match the bottom layer - interface. Note: the bottom_depth is positive, whereas the layer - interfaces are negative. - - Returns - ------- - diff : float - The computed bottom depth minus the target ``bottom_depth``. ``diff`` - should be zero when we have found the desired ``delta``. - """ - _, z = _cumsum_z(delta, nz, dz1, dz2) - diff = -bottom_depth - z[-1] - return diff - - -def _cumsum_z(delta, nz, dz1, dz2): - """ - Compute layer interface depths and layer thicknesses over ``nz`` layers - - Parameters - ---------- - delta : float - The characteristic length scale over which dz varies (this parameter - will be optimized to hit a target depth in a target number of layers) - - nz : int - The number of layers - - dz1 : float - The layer thickness at the top of the ocean (z = 0) - - dz2 : float - The layer thickness at z --> -infinity - - Returns - ------- - dz : numpy.ndarray - The layer thicknesses for each layer - - z : numpy.ndarray - The depth (positive up) of each layer interface (``nz + 1`` total - elements) - """ - dz = np.zeros(nz) - z = np.zeros(nz + 1) - for zindex in range(nz): - dz[zindex] = _dz_z(z[zindex], dz1, dz2, delta) - z[zindex + 1] = z[zindex] - dz[zindex] - return dz, z - - -def _dz_z(z, dz1, dz2, delta): - """ - layer thickness as a function of depth - - Parameters - ---------- - z : float - Depth coordinate (positive up) at which to find the layer thickness - - dz1 : float - The layer thickness at the top of the ocean (z = 0) - - dz2 : float - The layer thickness at z --> -infinity - - delta : float - The characteristic length scale over which dz varies (this parameter - will be optimized to hit a target depth in a target numer of layers) - - Returns - ------- - dz : float - The layer thickness - """ - return (dz2 - dz1) * np.tanh(-z * np.pi / delta) + dz1 From 84aef877eb5d469a9bb5a7d2c8f4dda95bc625fc Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 25 Aug 2023 14:14:39 -0500 Subject: [PATCH 14/16] Fixup output comparison paths --- .../ocean/tests/global_ocean/init/__init__.py | 15 ++++---- .../global_ocean/performance_test/__init__.py | 38 ++++++++++++------- 2 files changed, 33 insertions(+), 20 deletions(-) diff --git a/compass/ocean/tests/global_ocean/init/__init__.py b/compass/ocean/tests/global_ocean/init/__init__.py index 31309a1933..ea604224a1 100644 --- a/compass/ocean/tests/global_ocean/init/__init__.py +++ b/compass/ocean/tests/global_ocean/init/__init__.py @@ -54,7 +54,7 @@ def __init__(self, test_group, mesh, initial_condition, # Add the name of the subdir without inactive top cells whether or # not is has or will be run self.inactive_top_comp_subdir = init_subdir - init_subdir = os.path.join(init_subdir, inactive_top) + name = f'{name}_inactive_top' self.init_subdir = init_subdir subdir = os.path.join(self.init_subdir, name) super().__init__(test_group=test_group, name=name, subdir=subdir) @@ -66,7 +66,7 @@ def __init__(self, test_group, mesh, initial_condition, self.add_step( InitialState( test_case=self, mesh=mesh, - initial_condition=initial_condition, + initial_condition=initial_condition, with_inactive_top_cells=with_inactive_top_cells)) if mesh.with_ice_shelf_cavities: @@ -116,11 +116,12 @@ def validate(self): if os.path.exists(filename2): variables = ['temperature', 'salinity', 'layerThickness', 'normalVelocity'] - compare_variables(test_case=self, variables=variables, - filename1='initial_state/initial_state_crop.nc' - filename2=filename2, - quiet=False, check_outputs=False, - skip_if_step_not_run=False) + compare_variables( + test_case=self, variables=variables, + filename1='initial_state/initial_state_crop.nc', + filename2=filename2, + quiet=False, check_outputs=False, + skip_if_step_not_run=False) else: self.logger.warn('The version of "init" without inactive top ' diff --git a/compass/ocean/tests/global_ocean/performance_test/__init__.py b/compass/ocean/tests/global_ocean/performance_test/__init__.py index 59fcfd81d4..990fe397d7 100644 --- a/compass/ocean/tests/global_ocean/performance_test/__init__.py +++ b/compass/ocean/tests/global_ocean/performance_test/__init__.py @@ -1,10 +1,11 @@ import os + +from compass.ocean.inactive_top_cells import remove_inactive_top_cells_output from compass.ocean.tests.global_ocean.forward import ( ForwardStep, ForwardTestCase, ) from compass.validate import compare_timers, compare_variables -from compass.ocean.inactive_top_cells import remove_inactive_top_cells_output class PerformanceTest(ForwardTestCase): @@ -33,9 +34,18 @@ def __init__(self, test_group, mesh, init, time_integrator): time_integrator : {'split_explicit', 'RK4'} The time integrator to use for the forward run """ + name = 'performance_test' + if init.with_inactive_top_cells: + if time_integrator == 'RK4': + self.inactive_top_comp_subdir = os.path.join( + mesh.mesh_name, init.initial_condition, time_integrator, + name) + else: + self.inactive_top_comp_subdir = os.path.join( + mesh.mesh_name, init.initial_condition, name) + name = f'{name}_inactive_top' super().__init__(test_group=test_group, mesh=mesh, init=init, - time_integrator=time_integrator, - name='performance_test') + time_integrator=time_integrator, name=name) if mesh.with_ice_shelf_cavities: this_module = self.__module__ @@ -74,19 +84,21 @@ def validate(self): filename1=f'{step_subdir}/output.nc') if self.init.with_inactive_top_cells: - # construct the work directory for the other test - subdir = get_forward_subdir(self.init.inactive_top_comp_subdir, - self.time_integrator, self.name) - filename2 = os.path.join(self.base_work_dir, self.mpas_core.name, - self.test_group.name, subdir, + filename2 = os.path.join(self.base_work_dir, + self.mpas_core.name, + self.test_group.name, + self.inactive_top_comp_subdir, 'forward/output.nc') if os.path.exists(filename2): - compare_variables(test_case=self, variables=variables, - filename1='forward/output_crop.nc', - filename2=filename2) + compare_variables( + test_case=self, variables=variables, + filename1=f'{step_subdir}/output_crop.nc', + filename2=filename2, + quiet=False, check_outputs=False, + skip_if_step_not_run=False) else: - self.logger.warn('The version of "performance_test" without ' - 'inactive top cells was not run.\n' + self.logger.warn('The version of "performance_test" ' + 'without inactive top cells was not run. ' 'Skipping validation.') if self.mesh.with_ice_shelf_cavities: From 8308c57556eff461cf5cc2c252361e61fc009855 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 25 Aug 2023 14:15:28 -0500 Subject: [PATCH 15/16] Fixup vertical grid offset --- compass/ocean/vertical/grid_1d/__init__.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/compass/ocean/vertical/grid_1d/__init__.py b/compass/ocean/vertical/grid_1d/__init__.py index fab57b5f23..5731da0c90 100644 --- a/compass/ocean/vertical/grid_1d/__init__.py +++ b/compass/ocean/vertical/grid_1d/__init__.py @@ -28,16 +28,21 @@ def generate_1d_grid(config): A 1D array of positive depths for layer interfaces in meters """ section = config['vertical_grid'] + offset = 0 + if config.has_option('vertical_grid', 'inactive_top_cells'): + offset = section.getint('inactive_top_cells') + print(f'offset = {offset}') + grid_type = section.get('grid_type') if grid_type == 'uniform': vert_levels = section.getint('vert_levels') - interfaces = _generate_uniform(vert_levels) + interfaces = _generate_uniform(vert_levels - offset) elif grid_type == 'tanh_dz': vert_levels = section.getint('vert_levels') min_layer_thickness = section.getfloat('min_layer_thickness') max_layer_thickness = section.getfloat('max_layer_thickness') bottom_depth = section.getfloat('bottom_depth') - interfaces = create_tanh_dz_grid(vert_levels, + interfaces = create_tanh_dz_grid(vert_levels - offset, bottom_depth, min_layer_thickness, max_layer_thickness) @@ -66,6 +71,9 @@ def generate_1d_grid(config): # renormalize to the requested range interfaces = (bottom_depth / interfaces[-1]) * interfaces + if config.has_option('vertical_grid', 'inactive_top_cells'): + interfaces = np.append(np.zeros((offset)), interfaces) + return interfaces From c123941e0598c8b28047913aed62e58fea703753 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 25 Aug 2023 16:24:52 -0500 Subject: [PATCH 16/16] fixup inactive_top_cells config --- compass/ocean/inactive_top_cells.py | 4 ++-- compass/ocean/tests/global_ocean/init/__init__.py | 9 +++++++++ .../tests/global_ocean/init/initial_state.py | 15 ++------------- compass/ocean/tests/global_ocean/mesh/qu/icos.cfg | 6 ++++++ compass/ocean/tests/global_ocean/mesh/qu/qu.cfg | 2 ++ .../ocean/tests/global_ocean/mesh/qu240/qu240.cfg | 3 +++ 6 files changed, 24 insertions(+), 15 deletions(-) diff --git a/compass/ocean/inactive_top_cells.py b/compass/ocean/inactive_top_cells.py index c28599e68e..0f1ca7af62 100644 --- a/compass/ocean/inactive_top_cells.py +++ b/compass/ocean/inactive_top_cells.py @@ -1,6 +1,6 @@ import os -import xarray +import xarray from mpas_tools.io import write_netcdf @@ -40,6 +40,6 @@ def remove_inactive_top_cells_output(in_filename, out_filename=None, if minval != maxval: raise ValueError('Expected minLevelCell to have a constant ' 'value for inactive top cell tests') - ds_out = ds_in.isel(nVertLevels=slice(minval-1, None)) + ds_out = ds_in.isel(nVertLevels=slice(minval - 1, None)) write_netcdf(ds_out, out_filename) diff --git a/compass/ocean/tests/global_ocean/init/__init__.py b/compass/ocean/tests/global_ocean/init/__init__.py index ea604224a1..8d6081de6e 100644 --- a/compass/ocean/tests/global_ocean/init/__init__.py +++ b/compass/ocean/tests/global_ocean/init/__init__.py @@ -98,6 +98,15 @@ def configure(self, config=None): config.set('global_ocean', 'init_description', descriptions[initial_condition]) + if self.with_inactive_top_cells: + # Since we start at minLevelCell = 2, we need to increase the + # number of vertical levels in the cfg file to end up with the + # intended number in the initial state + vert_levels = config.getint('vertical_grid', 'vert_levels') + config.set('vertical_grid', 'vert_levels', f'{vert_levels + 1}', + comment='active vertical levels + inactive_top_cells') + config.set('vertical_grid', 'inactive_top_cells', '1') + def validate(self): """ Test cases can override this method to perform validation of variables diff --git a/compass/ocean/tests/global_ocean/init/initial_state.py b/compass/ocean/tests/global_ocean/init/initial_state.py index 4277e0fa82..a47b0d3e05 100644 --- a/compass/ocean/tests/global_ocean/init/initial_state.py +++ b/compass/ocean/tests/global_ocean/init/initial_state.py @@ -171,14 +171,6 @@ def run(self): Run this step of the testcase """ config = self.config - if self.with_inactive_top_cells: - # Since we start at minLevelCell = 2, we need to increase the - # number of vertical levels in the cfg file to end up with the - # intended number in the initial state - vert_levels = config.getint('vertical_grid', 'vert_levels') - config.set('vertical_grid', 'vert_levels', f'{vert_levels + 1}', - comment='the number of vertical levels + 1') - config.set('vertical_grid', 'inactive_top_cells', '1') logger = self.logger interfaces = generate_1d_grid(config=config) @@ -204,11 +196,8 @@ def run(self): ds = ds.isel(Time=0) - if config.has_option('vertical_grid', 'inactive_top_cells'): - offset = config.getint('vertical_grid', - 'inactive_top_cells') - else: - offset = 0 + offset = config.getint('vertical_grid', + 'inactive_top_cells') if 'minLevelCell' in ds: minLevelCell = ds.minLevelCell + offset diff --git a/compass/ocean/tests/global_ocean/mesh/qu/icos.cfg b/compass/ocean/tests/global_ocean/mesh/qu/icos.cfg index a1a448fea0..36774a8afc 100644 --- a/compass/ocean/tests/global_ocean/mesh/qu/icos.cfg +++ b/compass/ocean/tests/global_ocean/mesh/qu/icos.cfg @@ -9,3 +9,9 @@ prefix = Icos mesh_description = MPAS subdivided icosahedral mesh for E3SM version ${e3sm_version} at ${min_res}-km global resolution with <<>> vertical level + +# Options related to the vertical grid +[vertical_grid] + +# Number of inactive top cell layers +inactive_top_cells = 0 diff --git a/compass/ocean/tests/global_ocean/mesh/qu/qu.cfg b/compass/ocean/tests/global_ocean/mesh/qu/qu.cfg index eb4a478d49..7b20022a59 100644 --- a/compass/ocean/tests/global_ocean/mesh/qu/qu.cfg +++ b/compass/ocean/tests/global_ocean/mesh/qu/qu.cfg @@ -20,6 +20,8 @@ max_layer_thickness = 250.0 # the min and max occurs transition_levels = 28 +# Number of inactive top cell layers +inactive_top_cells = 0 # options for global ocean testcases [global_ocean] diff --git a/compass/ocean/tests/global_ocean/mesh/qu240/qu240.cfg b/compass/ocean/tests/global_ocean/mesh/qu240/qu240.cfg index de9bba14f9..4176455b76 100644 --- a/compass/ocean/tests/global_ocean/mesh/qu240/qu240.cfg +++ b/compass/ocean/tests/global_ocean/mesh/qu240/qu240.cfg @@ -16,6 +16,9 @@ min_layer_thickness = 3.0 # The maximum layer thickness max_layer_thickness = 500.0 +# Number of inactive top cell layers +inactive_top_cells = 0 + # options for spherical meshes [spherical_mesh]