From 22b0bbd4e819587c4b8e8ae3e1c6372f1496a410 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Mon, 4 Oct 2021 13:02:46 -0500 Subject: [PATCH 01/21] First commit of turbulence closure test case --- compass/ocean/__init__.py | 2 + .../tests/turbulence_closure/__init__.py | 54 +++++ .../turbulence_closure/boundary_values.py | 193 ++++++++++++++++++ .../decomp_test/__init__.py | 63 ++++++ .../turbulence_closure/default/__init__.py | 48 +++++ .../ocean/tests/turbulence_closure/forward.py | 84 ++++++++ .../tests/turbulence_closure/initial_state.py | 149 ++++++++++++++ .../tests/turbulence_closure/namelist.forward | 11 + .../restart_test/__init__.py | 71 +++++++ .../restart_test/namelist.full | 3 + .../restart_test/namelist.restart | 4 + .../restart_test/streams.full | 12 ++ .../restart_test/streams.restart | 12 ++ .../tests/turbulence_closure/streams.forward | 27 +++ .../turbulence_closure/turbulence_closure.cfg | 80 ++++++++ 15 files changed, 813 insertions(+) create mode 100644 compass/ocean/tests/turbulence_closure/__init__.py create mode 100644 compass/ocean/tests/turbulence_closure/boundary_values.py create mode 100644 compass/ocean/tests/turbulence_closure/decomp_test/__init__.py create mode 100644 compass/ocean/tests/turbulence_closure/default/__init__.py create mode 100644 compass/ocean/tests/turbulence_closure/forward.py create mode 100644 compass/ocean/tests/turbulence_closure/initial_state.py create mode 100644 compass/ocean/tests/turbulence_closure/namelist.forward create mode 100644 compass/ocean/tests/turbulence_closure/restart_test/__init__.py create mode 100644 compass/ocean/tests/turbulence_closure/restart_test/namelist.full create mode 100644 compass/ocean/tests/turbulence_closure/restart_test/namelist.restart create mode 100644 compass/ocean/tests/turbulence_closure/restart_test/streams.full create mode 100644 compass/ocean/tests/turbulence_closure/restart_test/streams.restart create mode 100644 compass/ocean/tests/turbulence_closure/streams.forward create mode 100644 compass/ocean/tests/turbulence_closure/turbulence_closure.cfg diff --git a/compass/ocean/__init__.py b/compass/ocean/__init__.py index 3c50d5125e..5d86c11962 100644 --- a/compass/ocean/__init__.py +++ b/compass/ocean/__init__.py @@ -1,4 +1,5 @@ from compass.mpas_core import MpasCore +from compass.ocean.tests.turbulence_closure import TurbulenceClosure from compass.ocean.tests.baroclinic_channel import BaroclinicChannel from compass.ocean.tests.dam_break import DamBreak from compass.ocean.tests.drying_slope import DryingSlope @@ -32,6 +33,7 @@ def __init__(self): """ super().__init__(name='ocean') + self.add_test_group(TurbulenceClosure(mpas_core=self)) self.add_test_group(BaroclinicChannel(mpas_core=self)) self.add_test_group(DamBreak(mpas_core=self)) self.add_test_group(DryingSlope(mpas_core=self)) diff --git a/compass/ocean/tests/turbulence_closure/__init__.py b/compass/ocean/tests/turbulence_closure/__init__.py new file mode 100644 index 0000000000..8f42731998 --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/__init__.py @@ -0,0 +1,54 @@ +from compass.testgroup import TestGroup +from compass.ocean.tests.turbulence_closure.decomp_test import DecompTest +from compass.ocean.tests.turbulence_closure.default import Default +from compass.ocean.tests.turbulence_closure.restart_test import RestartTest + + +class TurbulenceClosure(TestGroup): + """ + A test group for turbulence closure 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='turbulence_closure') + + for resolution in ['10km']: + self.add_test_case( + DecompTest(test_group=self, resolution=resolution)) + self.add_test_case( + Default(test_group=self, resolution=resolution)) + self.add_test_case( + RestartTest(test_group=self, resolution=resolution)) + + +def configure(resolution, config): + """ + Modify the configuration options for one of the turbulence closure test cases + + Parameters + ---------- + resolution : str + The resolution of the test case + + config : configparser.ConfigParser + Configuration options for this test case + """ + res_params = {'10km': {'nx': 16, + 'ny': 50, + 'dc': 10e3}, + '4km': {'nx': 40, + 'ny': 126, + 'dc': 4e3}, + '1km': {'nx': 160, + 'ny': 500, + 'dc': 1e3}} + + if resolution not in res_params: + raise ValueError('Unsupported resolution {}. Supported values are: ' + '{}'.format(resolution, list(res_params))) + res_params = res_params[resolution] + for param in res_params: + config.set('turbulence_closure', param, '{}'.format(res_params[param])) diff --git a/compass/ocean/tests/turbulence_closure/boundary_values.py b/compass/ocean/tests/turbulence_closure/boundary_values.py new file mode 100644 index 0000000000..1da8b94589 --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/boundary_values.py @@ -0,0 +1,193 @@ +import numpy + +def add_boundary_arrays(mesh, x_nonperiodic, y_nonperiodic): + """ + Computes boundary arrays necessary for LES test cases + Values are appended to a MPAS mesh dataset + + if the domain is doubly periodic, returns array of all zeros + as arrays are meaningless in that case. + + Parameters + ---------- + mesh : xarray dataset containing an MPAS mesh + + x_nonperiodic : logical if x direction is non periodic + + y_nonperiodic : logical if y direction is non periodic + """ + + nCells = mesh.dims['nCells'] + nEdges = mesh.dims['nEdges'] + neonc = mesh.nEdgesOnCell.values + eonc = mesh.edgesOnCell.values - 1 + conc = mesh.cellsOnCell.values - 1 + nVertLevels = mesh.dims['nVertLevels'] + mlc = mesh.maxLevelCell.values - 1 + cone = mesh.cellsOnEdge.values - 1 + dv = mesh.dvEdge.values + dc = mesh.dcEdge.values + + if (not x_nonperiodic) and (not y_nonperiodic): + mesh['distanceToBoundary'] = (['nCells', 'nVertLevels'], numpy.zeros((nCells,nVertLevels))) + mesh['boundaryNormalAngle'] = (['nCells', 'nVertLevels'], numpy.zeros((nCells,nVertLevels))) + mesh['boundaryCellMask'] = (['nCells', 'nVertLevels'], numpy.zeros((nCells,nVertLevels)).astype(int)) + return mesh + + #surface first and save until later + indVals = [] + angle = mesh.angleEdge.values + xc = mesh.xEdge.values + yc = mesh.yEdge.values + xCell = mesh.xCell.values + yCell = mesh.yCell.values + aCell = numpy.zeros(nCells) + bCell = numpy.zeros(nCells) + xEdge = [] + yEdge = [] + for i in range(nCells): + if conc[i,:].min() < 0: + aEdge = 0 + count = 0 + for j in range(neonc[i]): + if conc[i,j] < 0: + indVals.append(i) + xEdge.append(xc[eonc[i,j]]) + yEdge.append(yc[eonc[i,j]]) + #compute dx, dy + #check angle edge + #add pi if needed + #check if over 2pi correct if needed + dx = xc[eonc[i,j]] - xCell[i] + dy = yc[eonc[i,j]] - yCell[i] + if dx<=0 and dy<=0: # first quadrant + if angle[eonc[i,j]] >= numpy.pi/2.: + aEdge += angle[eonc[i,j]] - numpy.pi + count += 1 + else: + aEdge += angle[eonc[i,j]] + count += 1 + elif dx<=0 and dy>=0: # fourth quadrant, but reference to 0 not 2pi + if angle[eonc[i,j]] < 3.0*numpy.pi/2.0: + aEdge += angle[eonc[i,j]] + numpy.pi + count += 1 + else: + aEdge += angle[eonc[i,j]] + count += 1 + elif dx>=0 and dy>=0: #third quadrant + if angle[eonc[i,j]] < numpy.pi/2.0: # not in correct place + aEdge += angle[eonc[i,j]] + numpy.pi + count += 1 + else: + aEdge += angle[eonc[i,j]] + count += 1 + else: #quadrant 2 + if angle[eonc[i,j]] > 3.0*numpy.pi/2.0: #wrong place + aEdge += angle[eonc[i,j]] - numpy.pi + count += 1 + else: + aEdge += angle[eonc[i,j]] + count += 1 + if count > 0: + aCellTemp = aEdge / count + if aCellTemp > numpy.pi: + aCellTemp = 2.0*numpy.pi - aCellTemp + aCell[i] = aCellTemp + bCell[i] = 1 +#with angle and index arrays, now can do distance + + dist = numpy.zeros((nCells,nVertLevels)) + angleCells = numpy.zeros((nCells,nVertLevels)) + boundaryCells = numpy.zeros((nCells,nVertLevels)) + for i in range(nCells): + d = numpy.sqrt((xCell[i] - xEdge)**2 + (yCell[i] - yEdge)**2) + ind = d.argmin() + angleCells[i,0] = aCell[indVals[ind]] + boundaryCells[i,0] = bCell[i] + dist[i,0] = d.min() + + #now to the rest + for k in range(1,nVertLevels): + inds = numpy.where(k==mlc)[0] + edgeList = [] + cellList = [] + if len(inds) > 0: + #First step is forming list of edges bordering maxLC + for i in range(len(inds)): + for j in range(neonc[inds[i]]): + if conc[inds[i],j] not in inds: + edgeList.append(eonc[inds[i],j]) + for i in range(len(edgeList)): + c1 = cone[edgeList[i],0] + c2 = cone[edgeList[i],1] + if c1 in inds: + if c2 not in cellList: + cellList.append(c2) + if c2 in inds: + if c1 not in cellList: + cellList.append(c1) + + for i in range(len(cellList)): + aEdge = 0 + count = 0 + for j in range(neonc[cellList[i]]): + if eonc[cellList[i],j] in edgeList: + indVals.append(cellList[i]) + xEdge.append(xc[eonc[cellList[i],j]]) + yEdge.append(yc[eonc[cellList[i],j]]) + #compute dx, dy + #check angle edge + #add pi if needed + #check if over 2pi correct if needed + dx = xc[eonc[cellList[i],j]] - xCell[cellList[i]] + dy = yc[eonc[cellList[i],j]] - yCell[cellList[i]] + if dx<=0 and dy<=0: # first quadrant + if angle[eonc[cellList[i],j]] >= numpy.pi/2.: + aEdge += angle[eonc[cellList[i],j]] - numpy.pi + count += 1 + else: + aEdge += angle[eonc[cellList[i],j]] + count += 1 + elif dx<=0 and dy>=0: # fourth quadrant, but reference to 0 not 2pi + if angle[eonc[cellList[i],j]] < 3.0*numpy.pi/2.0: + aEdge += angle[eonc[cellList[i],j]] + numpy.pi + count += 1 + else: + aEdge += angle[eonc[cellList[i],j]] + count += 1 + elif dx>=0 and dy>=0: #third quadrant + if angle[eonc[cellList[i],j]] < numpy.pi/2.0: # not in correct place + aEdge += angle[eonc[cellList[i],j]] + numpy.pi + count += 1 + else: + aEdge += angle[eonc[cellList[i],j]] + count += 1 + else: #quadrant 2 + if angle[eonc[cellList[i],j]] > 3.0*numpy.pi/2.0: #wrong place + aEdge += angle[eonc[cellList[i],j]] - numpy.pi + count += 1 + else: + aEdge += angle[eonc[cellList[i],j]] + count += 1 + if count > 0: + aCellTemp = aEdge / count + if aCellTemp > numpy.pi: + aCellTemp = 2.0*numpy.pi - aCellTemp + aCell[cellList[i]] = aCellTemp + bCell[cellList[i]] = 1 + + for ii in range(nCells): + d = numpy.sqrt((xCell[ii] - xEdge)**2 + (yCell[ii] - yEdge)**2) + ind = d.argmin() + angleCells[ii,k] = aCell[indVals[ind]] + boundaryCells[ii,k] = bCell[ii] + dist[ii,k] = d.min() + else: + dist[:,k] = dist[:,0] + angleCells[:,k] = angleCells[:,0] + boundaryCells[:,k] = boundaryCells[:,0] + mesh['distanceToBoundary'] = (['nCells', 'nVertLevels'], dist) + mesh['boundaryNormalAngle'] = (['nCells', 'nVertLevels'], angleCells) + mesh['boundaryCellMask'] = (['nCells', 'nVertLevels'], boundaryCells.astype(int)) + + return mesh diff --git a/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py b/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py new file mode 100644 index 0000000000..6c3487faff --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py @@ -0,0 +1,63 @@ +from compass.testcase import TestCase +from compass.ocean.tests.turbulence_closure.initial_state import InitialState +from compass.ocean.tests.turbulence_closure.forward import Forward +from compass.ocean.tests import turbulence_closure +from compass.validate import compare_variables + + +class DecompTest(TestCase): + """ + A decomposition test case for the turbulence closure test group, which + makes sure the model produces identical results on 1 and 4 cores. + + Attributes + ---------- + resolution : str + The resolution of the test case + """ + + def __init__(self, test_group, resolution): + """ + Create the test case + + Parameters + ---------- + test_group : compass.ocean.tests.turbulence_closure.TurbulenceClosure + The test group that this test case belongs to + + resolution : str + The resolution of the test case + """ + name = 'decomp_test' + self.resolution = resolution + subdir = '{}/{}'.format(resolution, name) + super().__init__(test_group=test_group, name=name, + subdir=subdir) + + self.add_step( + InitialState(test_case=self, resolution=resolution)) + + for procs in [4, 8]: + name = '{}proc'.format(procs) + self.add_step( + Forward(test_case=self, name=name, subdir=name, cores=procs, + threads=1, resolution=resolution)) + + def configure(self): + """ + Modify the configuration options for this test case. + """ + turbulence_closure.configure(self.resolution, self.config) + + # no run() method is needed + + def validate(self): + """ + Test cases can override this method to perform validation of variables + and timers + """ + variables = ['temperature', 'salinity', 'layerThickness', + 'normalVelocity'] + compare_variables(test_case=self, variables=variables, + filename1='4proc/output.nc', + filename2='8proc/output.nc') diff --git a/compass/ocean/tests/turbulence_closure/default/__init__.py b/compass/ocean/tests/turbulence_closure/default/__init__.py new file mode 100644 index 0000000000..87b2958208 --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/default/__init__.py @@ -0,0 +1,48 @@ +from compass.testcase import TestCase +from compass.ocean.tests.turbulence_closure.initial_state import InitialState +from compass.ocean.tests.turbulence_closure.forward import Forward +from compass.ocean.tests import turbulence_closure + + +class Default(TestCase): + """ + The default test case for the turbulence closure test group simply creates + the mesh and initial condition, then performs a short forward run on 4 + cores. + + Attributes + ---------- + resolution : str + The resolution of the test case + """ + + def __init__(self, test_group, resolution): + """ + Create the test case + + Parameters + ---------- + test_group : compass.ocean.tests.turbulence_closure.TurbulenceClosure + The test group that this test case belongs to + + resolution : str + The resolution of the test case + """ + name = 'default' + self.resolution = resolution + subdir = '{}/{}'.format(resolution, name) + super().__init__(test_group=test_group, name=name, + subdir=subdir) + + self.add_step( + InitialState(test_case=self, resolution=resolution)) + self.add_step( + Forward(test_case=self, cores=4, threads=1, resolution=resolution)) + + def configure(self): + """ + Modify the configuration options for this test case. + """ + turbulence_closure.configure(self.resolution, self.config) + + # no run() is needed because we're doing the default: running all steps diff --git a/compass/ocean/tests/turbulence_closure/forward.py b/compass/ocean/tests/turbulence_closure/forward.py new file mode 100644 index 0000000000..33055a3098 --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/forward.py @@ -0,0 +1,84 @@ +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 baroclinic + channel test cases. + + Attributes + ---------- + resolution : str + The resolution of the test case + """ + def __init__(self, test_case, resolution, 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 + + resolution : str + The resolution of the test case + + 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) + """ + self.resolution = resolution + 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.baroclinic_channel', + 'namelist.forward') + self.add_namelist_file('compass.ocean.tests.baroclinic_channel', + 'namelist.{}.forward'.format(resolution)) + if nu is not None: + # update the viscosity to the requested value + options = {'config_mom_del2': '{}'.format(nu)} + self.add_namelist_options(options) + + # make sure output is double precision + self.add_streams_file('compass.ocean.streams', 'streams.output') + + self.add_streams_file('compass.ocean.tests.baroclinic_channel', + 'streams.forward') + + self.add_input_file(filename='init.nc', + target='../initial_state/ocean.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) diff --git a/compass/ocean/tests/turbulence_closure/initial_state.py b/compass/ocean/tests/turbulence_closure/initial_state.py new file mode 100644 index 0000000000..fd692cd852 --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/initial_state.py @@ -0,0 +1,149 @@ +import xarray +import numpy + +from compass.ocean.tests.turbulence_closure.boundary_values import add_boundary_arrays +from random import random, seed + +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 turbulence closure + test cases + + Attributes + ---------- + resolution : str + The resolution of the test case + """ + def __init__(self, test_case, resolution): + """ + Create the step + + Parameters + ---------- + test_case : compass.TestCase + The test case this step belongs to + + resolution : str + The resolution of the test case + """ + super().__init__(test_case=test_case, name='initial_state') + self.resolution = resolution + + for file in ['base_mesh.nc', 'culled_graph.info', 'culled_mesh.nc', + 'init_mode_forcing_data.nc', 'ocean.nc']: + self.add_output_file(file) + + def run(self): + """ + Run this step of the test case + """ + config = self.config + logger = self.logger + + section = config['turbulence_closure'] + nx = section.getint('nx') + ny = section.getint('ny') + dc = section.getfloat('dc') + + x_nonperiodic = section.getboolean('x_nonperiodic') + y_nonperiodic = section.getboolean('y_nonperiodic') + disturbance_amplitude = section.getfloat('disturbance_amplitude') + surface_temperature = section.getfloat('surface_temperature') + surface_salinity = section.getfloat('surface_salinity') + mixed_layer_salinity_gradient = section.getfloat('mixed_layer_salinity_gradient') + mixed_layer_temperature_gradient = section.getfloat('mixed_layer_temperature_gradient') + mixed_layer_depth_salinity = section.getfloat('mixed_layer_depth_salinity') + mixed_layer_depth_temperature = section.getfloat('mixed_layer_depth_temperature') + interior_salinity_gradient = section.getfloat('interior_salinity_gradient') + interior_temperature_gradient = section.getfloat('interior_temperature_gradient') + surface_heat_flux = section.getfloat('surface_heat_flux') + surface_freshwater_flux = section.getfloat('surface_freshwater_flux') + wind_stress_zonal = section.getfloat('wind_stress_zonal') + wind_stress_meridional = section.getfloat('wind_stress_meridional') + coriolis_parameter = section.getfloat('coriolis_parameter') + bottom_depth = config.getfloat('vertical_grid', 'bottom_depth') + + dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, nonperiodic_x=x_nonperiodic, + nonperiodic_y=y_nonperiodic) + + write_netcdf(dsMesh, 'base_mesh.nc') + + dsMesh = cull(dsMesh, logger=logger) + dsMesh = convert(dsMesh, graphInfoFileName='culled_graph.info', + logger=logger) + write_netcdf(dsMesh, 'culled_mesh.nc') + + ds = dsMesh.copy() + dsForcing = dsMesh.copy() + xCell = ds.xCell + yCell = ds.yCell + + ds['bottomDepth'] = bottom_depth * xarray.ones_like(xCell) + ds['ssh'] = xarray.zeros_like(xCell) + + init_vertical_coord(config, ds) + + dsForcing['windStressZonal'] = wind_stress_zonal * xarray.ones_like(xCell) + dsForcing['windStressMeridional'] = wind_stress_meridional * xarray.ones_like(xCell) + dsForcing['sensibleHeatFlux'] = surface_heat_flux * xarray.ones_like(xCell) + dsForcing['rainFlux'] = surface_freshwater_flux * xarray.ones_like(xCell) + + write_netcdf(dsForcing, 'init_mode_forcing_data.nc') + + normalVelocity = xarray.zeros_like(ds.xEdge) + normalVelocity, _ = xarray.broadcast(normalVelocity, ds.refBottomDepth) + normalVelocity = normalVelocity.transpose('nEdges', 'nVertLevels') + + temperature = xarray.zeros_like(ds.layerThickness) + + salinity = xarray.zeros_like(ds.layerThickness) + + surf_indices = numpy.where(ds.zMid[0,:,:].values >= -mixed_layer_depth_temperature) + + if len(surf_indices[0]) > 0: + temperature[0,surf_indices[0],surf_indices[1]] = surface_temperature + mixed_layer_temperature_gradient* \ + ds.zMid[0,surf_indices[0],surf_indices[1]].values + + int_indices = numpy.where(ds.zMid[0,:,:].values < -mixed_layer_depth_temperature) + + if len(int_indices[0]) > 0: + temperature[0,int_indices[0],int_indices[1]] = surface_temperature - mixed_layer_temperature_gradient* \ + mixed_layer_depth_temperature + interior_temperature_gradient* \ + (ds.zMid[0,int_indices[0],int_indices[1]] + \ + mixed_layer_depth_temperature) + + temperature[0,:,0] += disturbance_amplitude*2*(numpy.random.rand(len(ds.xCell.values)) - 0.5) + + surf_indices = numpy.where(ds.zMid[0,:,:].values >= -mixed_layer_depth_salinity) + if len(surf_indices[0]) > 0: + salinity[0,surf_indices[0],surf_indices[1]] = surface_salinity - mixed_layer_salinity_gradient * \ + ds.zMid[0,surf_indices[0],surf_indices[1]] + + int_indices = numpy.where(ds.zMid[0,:,:].values < -mixed_layer_depth_salinity) + if len(int_indices[0]) > 0: + salinity[0,int_indices[0],int_indices[1]] = surface_salinity - mixed_layer_salinity_gradient* \ + mixed_layer_depth_salinity + interior_salinity_gradient* \ + (ds.zMid[0,int_indices[0],int_indices[1]] + \ + mixed_layer_depth_salinity) + + normalVelocity = normalVelocity.expand_dims(dim='Time', axis=0) + + ds['normalVelocity'] = normalVelocity + ds['temperature'] = temperature + ds['salinity'] = salinity + ds['fCell'] = coriolis_parameter * xarray.ones_like(xCell) + ds['fEdge'] = coriolis_parameter * xarray.ones_like(ds.xEdge) + ds['fVertex'] = coriolis_parameter * xarray.ones_like(ds.xVertex) + + ds = add_boundary_arrays(ds, x_nonperiodic, y_nonperiodic) + + write_netcdf(ds, 'ocean.nc') + diff --git a/compass/ocean/tests/turbulence_closure/namelist.forward b/compass/ocean/tests/turbulence_closure/namelist.forward new file mode 100644 index 0000000000..c3f3d3780e --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/namelist.forward @@ -0,0 +1,11 @@ +config_write_output_on_startup = .false. +config_dt = '00:01:00' +config_time_integrator = 'split_implicit' +config_run_duration = '0001_00:00:00' +config_implicit_bottom_drag_coeff = 1.0e-2 +config_zonal_ssh_grad = 0.0 +config_meridional_ssh_grad = 0.0 +config_enable_nonhydrostatic_mode = .true. +config_use_constant_forced_pgf = .false. +config_use_two_equation_turbulence_Model = .true. +config_two_equation_model_choice = 'omega' diff --git a/compass/ocean/tests/turbulence_closure/restart_test/__init__.py b/compass/ocean/tests/turbulence_closure/restart_test/__init__.py new file mode 100644 index 0000000000..e8c64030af --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/restart_test/__init__.py @@ -0,0 +1,71 @@ +from compass.testcase import TestCase +from compass.ocean.tests.baroclinic_channel.initial_state import InitialState +from compass.ocean.tests.baroclinic_channel.forward import Forward +from compass.ocean.tests import baroclinic_channel +from compass.validate import compare_variables + + +class RestartTest(TestCase): + """ + A restart test case for the baroclinic channel test group, which makes sure + the model produces identical results with one longer run and two shorter + runs with a restart in between. + + Attributes + ---------- + resolution : str + The resolution of the test case + """ + + def __init__(self, test_group, resolution): + """ + Create the test case + + Parameters + ---------- + test_group : compass.ocean.tests.baroclinic_channel.BaroclinicChannel + The test group that this test case belongs to + + resolution : str + The resolution of the test case + """ + name = 'restart_test' + self.resolution = resolution + subdir = '{}/{}'.format(resolution, name) + super().__init__(test_group=test_group, name=name, + subdir=subdir) + + self.add_step( + InitialState(test_case=self, resolution=resolution)) + + for part in ['full', 'restart']: + name = '{}_run'.format(part) + step = Forward(test_case=self, name=name, subdir=name, cores=4, + threads=1, resolution=resolution) + + step.add_namelist_file( + 'compass.ocean.tests.baroclinic_channel.restart_test', + 'namelist.{}'.format(part)) + step.add_streams_file( + 'compass.ocean.tests.baroclinic_channel.restart_test', + 'streams.{}'.format(part)) + self.add_step(step) + + def configure(self): + """ + Modify the configuration options for this test case. + """ + baroclinic_channel.configure(self.resolution, self.config) + + # no run() method is needed + + def validate(self): + """ + Test cases can override this method to perform validation of variables + and timers + """ + variables = ['temperature', 'salinity', 'layerThickness', + 'normalVelocity'] + compare_variables(test_case=self, variables=variables, + filename1='full_run/output.nc', + filename2='restart_run/output.nc') diff --git a/compass/ocean/tests/turbulence_closure/restart_test/namelist.full b/compass/ocean/tests/turbulence_closure/restart_test/namelist.full new file mode 100644 index 0000000000..d99b25960d --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/restart_test/namelist.full @@ -0,0 +1,3 @@ +config_start_time = '0001-01-01_00:00:00' +config_run_duration = '0000_00:10:00' +config_write_output_on_startup = .false. diff --git a/compass/ocean/tests/turbulence_closure/restart_test/namelist.restart b/compass/ocean/tests/turbulence_closure/restart_test/namelist.restart new file mode 100644 index 0000000000..620cd94498 --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/restart_test/namelist.restart @@ -0,0 +1,4 @@ +config_start_time = '0001-01-01_00:05:00' +config_run_duration = '0000_00:05:00' +config_write_output_on_startup = .false. +config_do_restart = .true. diff --git a/compass/ocean/tests/turbulence_closure/restart_test/streams.full b/compass/ocean/tests/turbulence_closure/restart_test/streams.full new file mode 100644 index 0000000000..5d80d1abb5 --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/restart_test/streams.full @@ -0,0 +1,12 @@ + + + + + + + + diff --git a/compass/ocean/tests/turbulence_closure/restart_test/streams.restart b/compass/ocean/tests/turbulence_closure/restart_test/streams.restart new file mode 100644 index 0000000000..c99ec09db0 --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/restart_test/streams.restart @@ -0,0 +1,12 @@ + + + + + + + + diff --git a/compass/ocean/tests/turbulence_closure/streams.forward b/compass/ocean/tests/turbulence_closure/streams.forward new file mode 100644 index 0000000000..93d0b2a722 --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/streams.forward @@ -0,0 +1,27 @@ + + + + + + + + + + + + + + + + + + + + + diff --git a/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg b/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg new file mode 100644 index 0000000000..c276e32b66 --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg @@ -0,0 +1,80 @@ +# 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 = 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 + +# config options for baroclinic channel testcases +[turbulence_closure] + +# Grid cell count in x direction +nx = 4 + +# Grid cell count in y direction +ny = 4 + +# Grid cell size +dc = 50000 + +# Logical determining if grid is not periodic in x-direction +x_nonperiodic = False + +# Logical determining if grid is not periodic in y-direction +y_nonperiodic = False + +# Temperature at surface (^oC) +surface_temperature = 20 + +# Salinity at surface (PSU) +surface_salinity = 35 + +# Depth of mixed layer in temperature (m) +mixed_layer_depth_temperature = 0 + +# Depth of mixed layer in salinity (m) +mixed_layer_depth_salinity = 0 + +# Temperature gradient in the mixed layer (^oC/m) +mixed_layer_temperature_gradient = 0 + +# Salinity gradient in the mixed layer (PSU/m) +mixed_layer_salinity_gradient = 0 + +# Temperature gradient in the interior (^oC/m) +interior_temperature_gradient = 0.1 + +# Salinity gradient in the interior (PSU/m) +interior_salinity_gradient = 0 + +# Disturbance amplitude (added to initiate convection faster; ^oC) +disturbance_amplitude = 0.01 + +# Surface heat flux (W/m^2) +surface_heat_flux = -100 + +# Surface freshwater flux (kg / (m^2 s)) +surface_freshwater_flux = 0 + +# Surface zonal windstress (N/m^2) +wind_stress_zonal = 0 + +# Surface meiridional windstress (N/m^2) +wind_stress_meridional = 0 + +# Coriolis parameter for entire domain. +coriolis_parameter = -1.2e-4 From 19450c3c17c5d4d6289317b0ba15ffc9a0e10d0b Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Mon, 4 Oct 2021 13:11:06 -0500 Subject: [PATCH 02/21] fixes some default parameters --- compass/ocean/tests/turbulence_closure/forward.py | 10 +++++----- .../tests/turbulence_closure/namelist.10km.forward | 11 +++++++++++ 2 files changed, 16 insertions(+), 5 deletions(-) create mode 100644 compass/ocean/tests/turbulence_closure/namelist.10km.forward diff --git a/compass/ocean/tests/turbulence_closure/forward.py b/compass/ocean/tests/turbulence_closure/forward.py index 33055a3098..b20cb7c6d7 100644 --- a/compass/ocean/tests/turbulence_closure/forward.py +++ b/compass/ocean/tests/turbulence_closure/forward.py @@ -4,8 +4,8 @@ class Forward(Step): """ - A step for performing forward MPAS-Ocean runs as part of baroclinic - channel test cases. + A step for performing forward MPAS-Ocean runs as part of turbulence + closure test cases. Attributes ---------- @@ -51,9 +51,9 @@ def __init__(self, test_case, resolution, name='forward', subdir=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.baroclinic_channel', + self.add_namelist_file('compass.ocean.tests.turbulence_closure', 'namelist.forward') - self.add_namelist_file('compass.ocean.tests.baroclinic_channel', + self.add_namelist_file('compass.ocean.tests.turbulence_closure', 'namelist.{}.forward'.format(resolution)) if nu is not None: # update the viscosity to the requested value @@ -63,7 +63,7 @@ def __init__(self, test_case, resolution, name='forward', subdir=None, # make sure output is double precision self.add_streams_file('compass.ocean.streams', 'streams.output') - self.add_streams_file('compass.ocean.tests.baroclinic_channel', + self.add_streams_file('compass.ocean.tests.turbulence_closure', 'streams.forward') self.add_input_file(filename='init.nc', diff --git a/compass/ocean/tests/turbulence_closure/namelist.10km.forward b/compass/ocean/tests/turbulence_closure/namelist.10km.forward new file mode 100644 index 0000000000..c3f3d3780e --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/namelist.10km.forward @@ -0,0 +1,11 @@ +config_write_output_on_startup = .false. +config_dt = '00:01:00' +config_time_integrator = 'split_implicit' +config_run_duration = '0001_00:00:00' +config_implicit_bottom_drag_coeff = 1.0e-2 +config_zonal_ssh_grad = 0.0 +config_meridional_ssh_grad = 0.0 +config_enable_nonhydrostatic_mode = .true. +config_use_constant_forced_pgf = .false. +config_use_two_equation_turbulence_Model = .true. +config_two_equation_model_choice = 'omega' From 2ccb8f7afda269c4bc04b5b564355f5b92f94e98 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Mon, 11 Oct 2021 09:49:13 -0500 Subject: [PATCH 03/21] simplifies z dependence --- .../tests/turbulence_closure/initial_state.py | 30 +++++++++---------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/initial_state.py b/compass/ocean/tests/turbulence_closure/initial_state.py index fd692cd852..81c87687cf 100644 --- a/compass/ocean/tests/turbulence_closure/initial_state.py +++ b/compass/ocean/tests/turbulence_closure/initial_state.py @@ -106,32 +106,32 @@ def run(self): salinity = xarray.zeros_like(ds.layerThickness) - surf_indices = numpy.where(ds.zMid[0,:,:].values >= -mixed_layer_depth_temperature) + surf_indices = numpy.where(ds.refZMid.values >= -mixed_layer_depth_temperature)[0] - if len(surf_indices[0]) > 0: - temperature[0,surf_indices[0],surf_indices[1]] = surface_temperature + mixed_layer_temperature_gradient* \ - ds.zMid[0,surf_indices[0],surf_indices[1]].values + if len(surf_indices) > 0: + temperature[0,:,surf_indices] = surface_temperature + mixed_layer_temperature_gradient* \ + ds.zMid[0,:,surf_indices].values - int_indices = numpy.where(ds.zMid[0,:,:].values < -mixed_layer_depth_temperature) + int_indices = numpy.where(ds.refZMid.values < -mixed_layer_depth_temperature)[0] - if len(int_indices[0]) > 0: - temperature[0,int_indices[0],int_indices[1]] = surface_temperature - mixed_layer_temperature_gradient* \ + if len(int_indices) > 0: + temperature[0,:,int_indices] = surface_temperature - mixed_layer_temperature_gradient* \ mixed_layer_depth_temperature + interior_temperature_gradient* \ - (ds.zMid[0,int_indices[0],int_indices[1]] + \ + (ds.zMid[0,:,int_indices] + \ mixed_layer_depth_temperature) temperature[0,:,0] += disturbance_amplitude*2*(numpy.random.rand(len(ds.xCell.values)) - 0.5) - surf_indices = numpy.where(ds.zMid[0,:,:].values >= -mixed_layer_depth_salinity) + surf_indices = numpy.where(ds.refZMid.values >= -mixed_layer_depth_salinity)[0] if len(surf_indices[0]) > 0: - salinity[0,surf_indices[0],surf_indices[1]] = surface_salinity - mixed_layer_salinity_gradient * \ - ds.zMid[0,surf_indices[0],surf_indices[1]] + salinity[0,:,surf_indices] = surface_salinity - mixed_layer_salinity_gradient * \ + ds.zMid[0,:,surf_indices] - int_indices = numpy.where(ds.zMid[0,:,:].values < -mixed_layer_depth_salinity) - if len(int_indices[0]) > 0: - salinity[0,int_indices[0],int_indices[1]] = surface_salinity - mixed_layer_salinity_gradient* \ + int_indices = numpy.where(ds.refZMid.values < -mixed_layer_depth_salinity)[0] + if len(int_indices) > 0: + salinity[0,:,int_indices] = surface_salinity - mixed_layer_salinity_gradient* \ mixed_layer_depth_salinity + interior_salinity_gradient* \ - (ds.zMid[0,int_indices[0],int_indices[1]] + \ + (ds.zMid[0,:,int_indices] + \ mixed_layer_depth_salinity) normalVelocity = normalVelocity.expand_dims(dim='Time', axis=0) From a13156759b4b6688f9975ef3634276680c6b2223 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Thu, 10 Mar 2022 14:13:57 -0600 Subject: [PATCH 04/21] Fixes bug in initial_state --- .../ocean/tests/turbulence_closure/initial_state.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/initial_state.py b/compass/ocean/tests/turbulence_closure/initial_state.py index 81c87687cf..6eb1310a73 100644 --- a/compass/ocean/tests/turbulence_closure/initial_state.py +++ b/compass/ocean/tests/turbulence_closure/initial_state.py @@ -109,13 +109,13 @@ def run(self): surf_indices = numpy.where(ds.refZMid.values >= -mixed_layer_depth_temperature)[0] if len(surf_indices) > 0: - temperature[0,:,surf_indices] = surface_temperature + mixed_layer_temperature_gradient* \ + temperature[0,:,surf_indices] = surface_temperature + mixed_layer_temperature_gradient* \ ds.zMid[0,:,surf_indices].values int_indices = numpy.where(ds.refZMid.values < -mixed_layer_depth_temperature)[0] if len(int_indices) > 0: - temperature[0,:,int_indices] = surface_temperature - mixed_layer_temperature_gradient* \ + temperature[0,:,int_indices] = surface_temperature - mixed_layer_temperature_gradient* \ mixed_layer_depth_temperature + interior_temperature_gradient* \ (ds.zMid[0,:,int_indices] + \ mixed_layer_depth_temperature) @@ -123,13 +123,13 @@ def run(self): temperature[0,:,0] += disturbance_amplitude*2*(numpy.random.rand(len(ds.xCell.values)) - 0.5) surf_indices = numpy.where(ds.refZMid.values >= -mixed_layer_depth_salinity)[0] - if len(surf_indices[0]) > 0: - salinity[0,:,surf_indices] = surface_salinity - mixed_layer_salinity_gradient * \ + if len(surf_indices) > 0: + salinity[0,:,surf_indices] = surface_salinity - mixed_layer_salinity_gradient * \ ds.zMid[0,:,surf_indices] int_indices = numpy.where(ds.refZMid.values < -mixed_layer_depth_salinity)[0] if len(int_indices) > 0: - salinity[0,:,int_indices] = surface_salinity - mixed_layer_salinity_gradient* \ + salinity[0,:,int_indices] = surface_salinity - mixed_layer_salinity_gradient* \ mixed_layer_depth_salinity + interior_salinity_gradient* \ (ds.zMid[0,:,int_indices] + \ mixed_layer_depth_salinity) From efcd76206cf1787d13b087565a16877282a66d02 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 17 Nov 2022 17:26:23 -0600 Subject: [PATCH 05/21] Update task-related arguments --- .../decomp_test/__init__.py | 4 ++-- .../turbulence_closure/default/__init__.py | 2 +- .../ocean/tests/turbulence_closure/forward.py | 24 +++++++++---------- .../restart_test/__init__.py | 4 ++-- 4 files changed, 17 insertions(+), 17 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py b/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py index 6c3487faff..1cb18b63ae 100644 --- a/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py +++ b/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py @@ -40,8 +40,8 @@ def __init__(self, test_group, resolution): for procs in [4, 8]: name = '{}proc'.format(procs) self.add_step( - Forward(test_case=self, name=name, subdir=name, cores=procs, - threads=1, resolution=resolution)) + Forward(test_case=self, name=name, subdir=name, ntasks=procs, + openmp_threads=1, resolution=resolution)) def configure(self): """ diff --git a/compass/ocean/tests/turbulence_closure/default/__init__.py b/compass/ocean/tests/turbulence_closure/default/__init__.py index 87b2958208..e2fa0e0a22 100644 --- a/compass/ocean/tests/turbulence_closure/default/__init__.py +++ b/compass/ocean/tests/turbulence_closure/default/__init__.py @@ -37,7 +37,7 @@ def __init__(self, test_group, resolution): self.add_step( InitialState(test_case=self, resolution=resolution)) self.add_step( - Forward(test_case=self, cores=4, threads=1, resolution=resolution)) + Forward(test_case=self, ntasks=4, openmp_threads=1, resolution=resolution)) def configure(self): """ diff --git a/compass/ocean/tests/turbulence_closure/forward.py b/compass/ocean/tests/turbulence_closure/forward.py index b20cb7c6d7..c13478ec7d 100644 --- a/compass/ocean/tests/turbulence_closure/forward.py +++ b/compass/ocean/tests/turbulence_closure/forward.py @@ -13,7 +13,7 @@ class Forward(Step): The resolution of the test case """ def __init__(self, test_case, resolution, name='forward', subdir=None, - cores=1, min_cores=None, threads=1, nu=None): + ntasks=1, min_tasks=None, openmp_threads=1, nu=None): """ Create a new test case @@ -31,26 +31,26 @@ def __init__(self, test_case, resolution, name='forward', subdir=None, 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 + ntasks: int, optional + the number of tasks the step would ideally use. If fewer tasks are available on the system, the step will run on all available - cores as long as this is not below ``min_cores`` + cores as long as this is not below ``min_tasks`` - 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 + min_tasks : int, optional + the number of tasks the step requires. If the system has fewer + than this number of tasks, the step will fail - threads : int, optional - the number of threads the step will use + openmp_threads : int, optional + the number of OpenMP threads the step will use nu : float, optional the viscosity (if different from the default for the test group) """ self.resolution = resolution - if min_cores is None: - min_cores = cores + if min_tasks is None: + min_tasks = ntasks super().__init__(test_case=test_case, name=name, subdir=subdir, - cores=cores, min_cores=min_cores, threads=threads) + ntasks=ntasks, min_tasks=min_tasks, openmp_threads=openmp_threads) self.add_namelist_file('compass.ocean.tests.turbulence_closure', 'namelist.forward') self.add_namelist_file('compass.ocean.tests.turbulence_closure', diff --git a/compass/ocean/tests/turbulence_closure/restart_test/__init__.py b/compass/ocean/tests/turbulence_closure/restart_test/__init__.py index e8c64030af..6e075545e9 100644 --- a/compass/ocean/tests/turbulence_closure/restart_test/__init__.py +++ b/compass/ocean/tests/turbulence_closure/restart_test/__init__.py @@ -40,8 +40,8 @@ def __init__(self, test_group, resolution): for part in ['full', 'restart']: name = '{}_run'.format(part) - step = Forward(test_case=self, name=name, subdir=name, cores=4, - threads=1, resolution=resolution) + step = Forward(test_case=self, name=name, subdir=name, ntasks=4, + openmp_threads=1, resolution=resolution) step.add_namelist_file( 'compass.ocean.tests.baroclinic_channel.restart_test', From 8426c7305d3fe88879851a761a7964b412f46ed9 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 17 Nov 2022 17:44:08 -0600 Subject: [PATCH 06/21] Add a few vars to streams and change restart interval --- .../tests/turbulence_closure/streams.forward | 47 ++++++++++++++++--- 1 file changed, 41 insertions(+), 6 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/streams.forward b/compass/ocean/tests/turbulence_closure/streams.forward index 93d0b2a722..fb3a781eb6 100644 --- a/compass/ocean/tests/turbulence_closure/streams.forward +++ b/compass/ocean/tests/turbulence_closure/streams.forward @@ -6,22 +6,57 @@ - + + + + + + + + + + + + + + + + - - - - + + + + + + + + + + + + + + + + + + + + + + + From 5c5371e7d82357b2615e80b2673d2b9a0f6ce809 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 18 Nov 2022 17:27:16 -0600 Subject: [PATCH 07/21] Add viz --- .../turbulence_closure/default/__init__.py | 2 + compass/ocean/tests/turbulence_closure/viz.py | 219 ++++++++++++++++++ 2 files changed, 221 insertions(+) create mode 100644 compass/ocean/tests/turbulence_closure/viz.py diff --git a/compass/ocean/tests/turbulence_closure/default/__init__.py b/compass/ocean/tests/turbulence_closure/default/__init__.py index e2fa0e0a22..6238fe7500 100644 --- a/compass/ocean/tests/turbulence_closure/default/__init__.py +++ b/compass/ocean/tests/turbulence_closure/default/__init__.py @@ -1,6 +1,7 @@ from compass.testcase import TestCase from compass.ocean.tests.turbulence_closure.initial_state import InitialState from compass.ocean.tests.turbulence_closure.forward import Forward +from compass.ocean.tests.turbulence_closure.viz import Viz from compass.ocean.tests import turbulence_closure @@ -38,6 +39,7 @@ def __init__(self, test_group, resolution): InitialState(test_case=self, resolution=resolution)) self.add_step( Forward(test_case=self, ntasks=4, openmp_threads=1, resolution=resolution)) + self.add_step(Viz(test_case=self)) def configure(self): """ diff --git a/compass/ocean/tests/turbulence_closure/viz.py b/compass/ocean/tests/turbulence_closure/viz.py new file mode 100644 index 0000000000..5b0a37d4b5 --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/viz.py @@ -0,0 +1,219 @@ +import xarray +import numpy +import matplotlib.pyplot as plt + +from compass.step import Step +import cmocean + + +class Viz(Step): + """ + A step for visualizing a cross-section through the domain + """ + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.TestCase + The test case this step belongs to + """ + super().__init__(test_case=test_case, name='viz') + + self.add_input_file(filename='output.nc', + target='../forward/output.nc') + self.add_output_file('uNormal_depth_section_t0.png') + self.add_output_file('pt_depth_section_t0.png') + self.add_output_file('sa_depth_section_t0.png') + self.add_output_file('layerThickness_depth_section_t0.png') + + def run(self): + """ + Run this step of the test case + """ + + ds = xarray.open_dataset('output.nc') + ds = ds.sortby('yEdge') + + dsInit = xarray.open_dataset('../forward/init.nc') + ds0 = ds.isel(Time=0) + figsize = [6.4, 4.8] + markersize = 20 + if 'Time' not in ds.dims: + print('Dataset missing time dimension') + return + nSteps = ds.sizes['Time'] # number of timesteps + tend = nSteps - 1 + + nCells = dsInit.sizes['nCells'] + nEdges = dsInit.sizes['nEdges'] + nVertLevels = dsInit.sizes['nVertLevels'] + + xEdge = dsInit.xEdge + yEdge = dsInit.yEdge + xCell = dsInit.xCell + yCell = dsInit.yCell + + xEdge_mid = numpy.median(xEdge) + edgeMask_x = numpy.equal(xEdge, xEdge_mid) + + # Solve for lateral boundaries of uNormal at cell centers for + # x-section + cellsOnEdge = dsInit.cellsOnEdge + cellsOnEdge_x = cellsOnEdge[edgeMask_x, :] + yEdges = numpy.zeros((len(cellsOnEdge_x)+1)) + for i in range(len(cellsOnEdge_x)): + if cellsOnEdge[i, 1] == 0: + yEdges[i] = yCell[cellsOnEdge_x[i, 0] - 1] + yEdges[i+1] = yCell[cellsOnEdge_x[i, 0] - 1] + elif cellsOnEdge[i, 1] == 0: + yEdges[i] = yCell[cellsOnEdge_x[i, 1] - 1] + yEdges[i+1] = yCell[cellsOnEdge_x[i, 1] - 1] + else: + yEdges[i] = min(yCell[cellsOnEdge_x[i, 0] - 1], + yCell[cellsOnEdge_x[i, 1] - 1]) + yEdges[i+1] = max(yCell[cellsOnEdge_x[i, 0] - 1], + yCell[cellsOnEdge_x[i, 1] - 1]) + + # Prep variables for cell quantities + cellIndex = numpy.subtract(cellsOnEdge_x[1:, 0], 1) + yEdge_x = yEdge[edgeMask_x] + + # prep all variables for uNormal plot + zInterface = dsInit.refInterfaces + zMid = dsInit.refZMid + + zInterfaces_edge_mesh, yEdges_mesh = numpy.meshgrid(zInterface, + yEdges) + zInterfaces_cell_mesh, yCells_mesh = numpy.meshgrid(zInterface, + yEdge_x) + temperature0 = ds0.temperature + temperature0_z = temperature0.mean(dim='nCells') + + for j in [tend]: + ds1 = ds.isel(Time=j) + + normalVelocity = ds1.normalVelocity + normalVelocity_xmesh = normalVelocity[edgeMask_x, :] + + + # Import cell quantities + layerThickness = ds1.layerThickness + layerThickness_x = layerThickness[cellIndex, :] + temperature = ds1.temperature + temperature_z = temperature.mean(dim='nCells') + temperature_x = temperature[cellIndex, :] + salinity = ds1.salinity + salinity_z = salinity.mean(dim='nCells') + salinity_x = salinity[cellIndex, :] + w = ds1.vertVelocityTop + w_x = w[cellIndex, :] + w_zTop = w[:, 0] + + # Figures + plt.figure(figsize=figsize, dpi=100) + cmax = numpy.max(numpy.abs(normalVelocity_xmesh.values)) + plt.pcolormesh(numpy.divide(yEdges_mesh, 1e3), + zInterfaces_edge_mesh, + normalVelocity_xmesh.values, + cmap='cmo.balance', vmin=-1.*cmax, vmax=cmax) + cbar = plt.colorbar() + cbar.ax.set_title('uNormal (m/s)') + plt.savefig('uNormal_depth_section_t{}.png'.format(j), + bbox_inches='tight', dpi=200) + plt.close() + + # ------------------------------------------------------------------ + # Plot cell-centered variables + # ------------------------------------------------------------------ + # Figures + plt.figure(figsize=figsize, dpi=100) + plt.plot(temperature_z.values, zMid) + plt.xlabel('PT (C)') + plt.ylabel('z (m)') + plt.savefig('pt_depth_t{}.png'.format(j), + bbox_inches='tight', dpi=200) + plt.close() + + plt.figure(figsize=figsize, dpi=100) + plt.pcolormesh(numpy.divide(yCells_mesh, 1e3), + zInterfaces_cell_mesh, + temperature_x.values, cmap='cmo.thermal') + plt.xlabel('y (km)') + plt.ylabel('z (m)') + cbar = plt.colorbar() + cbar.ax.set_title('PT (C)') + plt.savefig('pt_depth_section_t{}.png'.format(j), + bbox_inches='tight', dpi=200) + plt.close() + + plt.figure(figsize=figsize, dpi=100) + plt.plot(numpy.subtract(temperature_z.values, + temperature0_z.values), + zMid) + plt.xlabel('delta PT (C)') + plt.ylabel('z (m)') + plt.savefig('dpt_depth_t{}.png'.format(j), + bbox_inches='tight', dpi=200) + plt.close() + + plt.figure(figsize=figsize, dpi=100) + plt.plot(salinity_z.values, zMid) + plt.xlabel('SA (g/kg)') + plt.ylabel('z (m)') + plt.savefig('sa_depth_t{}.png'.format(j), + bbox_inches='tight', dpi=200) + plt.close() + + plt.figure(figsize=figsize, dpi=100) + plt.pcolormesh(numpy.divide(yCells_mesh, 1e3), + zInterfaces_cell_mesh, + salinity_x.values, cmap='cmo.haline') + plt.xlabel('y (km)') + plt.ylabel('z (m)') + cbar = plt.colorbar() + cbar.ax.set_title('SA (g/kg)') + plt.savefig('sa_depth_section_t{}.png'.format(j), + bbox_inches='tight', dpi=200) + plt.close() + + plt.figure(figsize=figsize, dpi=100) + cmax = numpy.max(numpy.abs(w_x.values)) + plt.pcolormesh(numpy.divide(yCells_mesh[:-1,:], 1e3), + zInterfaces_cell_mesh[:-1,:], + w_x.values, cmap='cmo.balance', + vmin=-1.*cmax, vmax=cmax) + plt.xlabel('y (km)') + plt.ylabel('z (m)') + cbar = plt.colorbar() + cbar.ax.set_title('w (m/s)') + plt.savefig('w_depth_section_t{}.png'.format(j), + bbox_inches='tight', dpi=200) + plt.close() + + plt.figure(figsize=figsize, dpi=100) + cmax = numpy.max(numpy.abs(w_zTop.values)) + plt.scatter(numpy.divide(xCell, 1e3), + numpy.divide(yCell, 1e3), + s=5, c=w_zTop.values, + cmap='cmo.balance', vmin=-1.*cmax, vmax=cmax) + plt.xlabel('x (km)') + plt.ylabel('y (km)') + cbar = plt.colorbar() + cbar.ax.set_title('w (m/s)') + plt.savefig('w_top_section_t{}.png'.format(j), + bbox_inches='tight', dpi=200) + plt.close() + + plt.figure(figsize=figsize, dpi=100) + plt.pcolormesh(numpy.divide(yCells_mesh, 1e3), + zInterfaces_cell_mesh, + layerThickness_x.values, cmap='viridis') + plt.xlabel('y (km)') + plt.ylabel('z (m)') + cbar = plt.colorbar() + cbar.ax.set_title('h (m)') + plt.savefig('layerThickness_depth_section_t{}.png'.format(j), + bbox_inches='tight', dpi=200) + plt.close() From 6ba55624b3455f2b4a18d3f5e9420985583cf441 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 22 Nov 2022 16:46:26 -0600 Subject: [PATCH 08/21] Fixup forcing --- compass/ocean/tests/turbulence_closure/forward.py | 2 ++ compass/ocean/tests/turbulence_closure/namelist.forward | 1 + compass/ocean/tests/turbulence_closure/streams.forward | 7 +++++++ 3 files changed, 10 insertions(+) diff --git a/compass/ocean/tests/turbulence_closure/forward.py b/compass/ocean/tests/turbulence_closure/forward.py index c13478ec7d..97d15a721f 100644 --- a/compass/ocean/tests/turbulence_closure/forward.py +++ b/compass/ocean/tests/turbulence_closure/forward.py @@ -68,6 +68,8 @@ def __init__(self, test_case, resolution, name='forward', subdir=None, self.add_input_file(filename='init.nc', target='../initial_state/ocean.nc') + self.add_input_file(filename='forcing.nc', + target='../initial_state/init_mode_forcing_data.nc') self.add_input_file(filename='graph.info', target='../initial_state/culled_graph.info') diff --git a/compass/ocean/tests/turbulence_closure/namelist.forward b/compass/ocean/tests/turbulence_closure/namelist.forward index c3f3d3780e..10a1a577fe 100644 --- a/compass/ocean/tests/turbulence_closure/namelist.forward +++ b/compass/ocean/tests/turbulence_closure/namelist.forward @@ -9,3 +9,4 @@ config_enable_nonhydrostatic_mode = .true. config_use_constant_forced_pgf = .false. config_use_two_equation_turbulence_Model = .true. config_two_equation_model_choice = 'omega' +config_use_activeTracers_surface_bulk_forcing = .true. diff --git a/compass/ocean/tests/turbulence_closure/streams.forward b/compass/ocean/tests/turbulence_closure/streams.forward index fb3a781eb6..dc6a77fe5f 100644 --- a/compass/ocean/tests/turbulence_closure/streams.forward +++ b/compass/ocean/tests/turbulence_closure/streams.forward @@ -9,6 +9,13 @@ + + + + Date: Tue, 22 Nov 2022 16:46:44 -0600 Subject: [PATCH 09/21] Add fine resolution case --- .../ocean/tests/turbulence_closure/__init__.py | 7 +++++-- .../ocean/tests/turbulence_closure/forward.py | 18 +++++++++--------- .../turbulence_closure/namelist.10km.forward | 11 ----------- .../tests/turbulence_closure/namelist.forward | 14 +++++++++++++- 4 files changed, 27 insertions(+), 23 deletions(-) delete mode 100644 compass/ocean/tests/turbulence_closure/namelist.10km.forward diff --git a/compass/ocean/tests/turbulence_closure/__init__.py b/compass/ocean/tests/turbulence_closure/__init__.py index 8f42731998..705ab306ea 100644 --- a/compass/ocean/tests/turbulence_closure/__init__.py +++ b/compass/ocean/tests/turbulence_closure/__init__.py @@ -15,7 +15,7 @@ def __init__(self, mpas_core): """ super().__init__(mpas_core=mpas_core, name='turbulence_closure') - for resolution in ['10km']: + for resolution in ['50m','10km']: self.add_test_case( DecompTest(test_group=self, resolution=resolution)) self.add_test_case( @@ -44,7 +44,10 @@ def configure(resolution, config): 'dc': 4e3}, '1km': {'nx': 160, 'ny': 500, - 'dc': 1e3}} + 'dc': 1e3}, + '50m': {'nx': 3200, + 'ny': 10000, + 'dc': 5e1}} if resolution not in res_params: raise ValueError('Unsupported resolution {}. Supported values are: ' diff --git a/compass/ocean/tests/turbulence_closure/forward.py b/compass/ocean/tests/turbulence_closure/forward.py index 97d15a721f..ca98307692 100644 --- a/compass/ocean/tests/turbulence_closure/forward.py +++ b/compass/ocean/tests/turbulence_closure/forward.py @@ -42,9 +42,6 @@ def __init__(self, test_case, resolution, name='forward', subdir=None, openmp_threads : int, optional the number of OpenMP threads the step will use - - nu : float, optional - the viscosity (if different from the default for the test group) """ self.resolution = resolution if min_tasks is None: @@ -53,12 +50,6 @@ def __init__(self, test_case, resolution, name='forward', subdir=None, ntasks=ntasks, min_tasks=min_tasks, openmp_threads=openmp_threads) self.add_namelist_file('compass.ocean.tests.turbulence_closure', 'namelist.forward') - self.add_namelist_file('compass.ocean.tests.turbulence_closure', - 'namelist.{}.forward'.format(resolution)) - if nu is not None: - # update the viscosity to the requested value - options = {'config_mom_del2': '{}'.format(nu)} - self.add_namelist_options(options) # make sure output is double precision self.add_streams_file('compass.ocean.streams', 'streams.output') @@ -83,4 +74,13 @@ def run(self): """ Run this step of the test case """ + # update the time step + resolution = self.resolution + if resolution == '10km': + self.update_namelist_at_runtime({'config_dt': + "'0000_00:01:00'"}) + elif resolution == '1m': + self.update_namelist_at_runtime({'config_dt': + "'0000_00:00:00.5'"}) + run_model(self) diff --git a/compass/ocean/tests/turbulence_closure/namelist.10km.forward b/compass/ocean/tests/turbulence_closure/namelist.10km.forward deleted file mode 100644 index c3f3d3780e..0000000000 --- a/compass/ocean/tests/turbulence_closure/namelist.10km.forward +++ /dev/null @@ -1,11 +0,0 @@ -config_write_output_on_startup = .false. -config_dt = '00:01:00' -config_time_integrator = 'split_implicit' -config_run_duration = '0001_00:00:00' -config_implicit_bottom_drag_coeff = 1.0e-2 -config_zonal_ssh_grad = 0.0 -config_meridional_ssh_grad = 0.0 -config_enable_nonhydrostatic_mode = .true. -config_use_constant_forced_pgf = .false. -config_use_two_equation_turbulence_Model = .true. -config_two_equation_model_choice = 'omega' diff --git a/compass/ocean/tests/turbulence_closure/namelist.forward b/compass/ocean/tests/turbulence_closure/namelist.forward index 10a1a577fe..659db76d9e 100644 --- a/compass/ocean/tests/turbulence_closure/namelist.forward +++ b/compass/ocean/tests/turbulence_closure/namelist.forward @@ -1,12 +1,24 @@ config_write_output_on_startup = .false. -config_dt = '00:01:00' config_time_integrator = 'split_implicit' config_run_duration = '0001_00:00:00' config_implicit_bottom_drag_coeff = 1.0e-2 config_zonal_ssh_grad = 0.0 config_meridional_ssh_grad = 0.0 +config_aust_scale_factor = 0.70710678 +config_use_mom_div2 = .true. +config_mom_div2 = .00002 +config_use_mom_div4 = .false. +config_mom_div4 = 1.0e12 config_enable_nonhydrostatic_mode = .true. +config_nonhydrostatic_preconditioner = 'bjacobi' +config_nonhydrostatic_solver_type = 'cg' +config_nonhydrostatic_solve_surface_boundary_condition = 'pressureTopGradientBottom' config_use_constant_forced_pgf = .false. config_use_two_equation_turbulence_Model = .true. config_two_equation_model_choice = 'omega' config_use_activeTracers_surface_bulk_forcing = .true. +config_LES_mode = .true. +config_mom_del2 = .true. +config_mom_del2 = 2e-6 +config_use_cvmix = .false. +config_btr_si_partition_match_mode = .true. From 4a93fa91b0134c24a7586a6425f4b285d5b79f2f Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 22 Dec 2022 14:12:01 -0600 Subject: [PATCH 10/21] Remove grid specs from cfg file --- .../tests/turbulence_closure/namelist.forward | 2 +- .../turbulence_closure/turbulence_closure.cfg | 15 --------------- 2 files changed, 1 insertion(+), 16 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/namelist.forward b/compass/ocean/tests/turbulence_closure/namelist.forward index 659db76d9e..a003b266ab 100644 --- a/compass/ocean/tests/turbulence_closure/namelist.forward +++ b/compass/ocean/tests/turbulence_closure/namelist.forward @@ -14,7 +14,7 @@ config_nonhydrostatic_preconditioner = 'bjacobi' config_nonhydrostatic_solver_type = 'cg' config_nonhydrostatic_solve_surface_boundary_condition = 'pressureTopGradientBottom' config_use_constant_forced_pgf = .false. -config_use_two_equation_turbulence_Model = .true. +config_use_two_equation_turbulence_model = .true. config_two_equation_model_choice = 'omega' config_use_activeTracers_surface_bulk_forcing = .true. config_LES_mode = .true. diff --git a/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg b/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg index c276e32b66..f29a1bc4f8 100644 --- a/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg +++ b/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg @@ -4,12 +4,6 @@ # the type of vertical grid grid_type = uniform -# Number of vertical levels -vert_levels = 20 - -# 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 @@ -22,15 +16,6 @@ min_pc_fraction = 0.1 # config options for baroclinic channel testcases [turbulence_closure] -# Grid cell count in x direction -nx = 4 - -# Grid cell count in y direction -ny = 4 - -# Grid cell size -dc = 50000 - # Logical determining if grid is not periodic in x-direction x_nonperiodic = False From ffa460df3c94b2e2c51d83783f492e73ab60b8f8 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 22 Dec 2022 14:13:22 -0600 Subject: [PATCH 11/21] Add forcing options --- .../tests/turbulence_closure/__init__.py | 49 +++++++++++++------ .../decomp_test/__init__.py | 4 +- .../turbulence_closure/default/__init__.py | 19 +++++-- .../restart_test/__init__.py | 4 +- 4 files changed, 51 insertions(+), 25 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/__init__.py b/compass/ocean/tests/turbulence_closure/__init__.py index 705ab306ea..957edca7b4 100644 --- a/compass/ocean/tests/turbulence_closure/__init__.py +++ b/compass/ocean/tests/turbulence_closure/__init__.py @@ -15,16 +15,18 @@ def __init__(self, mpas_core): """ super().__init__(mpas_core=mpas_core, name='turbulence_closure') - for resolution in ['50m','10km']: + for resolution in ['10km']: self.add_test_case( DecompTest(test_group=self, resolution=resolution)) - self.add_test_case( - Default(test_group=self, resolution=resolution)) self.add_test_case( RestartTest(test_group=self, resolution=resolution)) + for resolution in ['1m', '10km']: + for forcing in ['cooling', 'evaporation']: + self.add_test_case( + Default(test_group=self, resolution=resolution, forcing=forcing)) -def configure(resolution, config): +def configure(resolution, forcing, config): """ Modify the configuration options for one of the turbulence closure test cases @@ -39,19 +41,34 @@ def configure(resolution, config): res_params = {'10km': {'nx': 16, 'ny': 50, 'dc': 10e3}, - '4km': {'nx': 40, - 'ny': 126, - 'dc': 4e3}, - '1km': {'nx': 160, - 'ny': 500, - 'dc': 1e3}, - '50m': {'nx': 3200, - 'ny': 10000, - 'dc': 5e1}} + '1m': {'nx': 128, + 'ny': 128, + 'dc': 1}} + vert_params = {'10km': {'vert_levels': 20, + 'bottom_depth': 1e3}, + '1m': {'vert_levels': 128, + 'bottom_depth': 128.0}} if resolution not in res_params: - raise ValueError('Unsupported resolution {}. Supported values are: ' - '{}'.format(resolution, list(res_params))) + raise ValueError(f'Unsupported resolution {resolution}. ' + f'Supported values are: {list(res_params)}') + res_params = res_params[resolution] for param in res_params: - config.set('turbulence_closure', param, '{}'.format(res_params[param])) + config.set('turbulence_closure', param, f'{res_params[param]}') + vert_params = vert_params[resolution] + for param in vert_params: + config.set('vertical_grid', param, f'{vert_params[param]}') + + if forcing == 'cooling': + config.set('turbulence_closure', 'surface_heat_flux', '-100') + config.set('turbulence_closure', 'surface_freshwater_flux', '0') + config.set('turbulence_closure', 'interior_temperature_gradient', '0.1') + config.set('turbulence_closure', 'interior_salinity_gradient', '0') + config.set('turbulence_closure', 'wind_stress_zonal', '0') + if forcing == 'evaporation': + config.set('turbulence_closure', 'surface_heat_flux', '0') + config.set('turbulence_closure', 'surface_freshwater_flux', '0.429') + config.set('turbulence_closure', 'interior_temperature_gradient', '0') + config.set('turbulence_closure', 'interior_salinity_gradient', '-0.025') + config.set('turbulence_closure', 'wind_stress_zonal', '0') diff --git a/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py b/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py index 1cb18b63ae..8d26a51148 100644 --- a/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py +++ b/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py @@ -16,7 +16,7 @@ class DecompTest(TestCase): The resolution of the test case """ - def __init__(self, test_group, resolution): + def __init__(self, test_group, resolution, forcing='cooling'): """ Create the test case @@ -30,7 +30,7 @@ def __init__(self, test_group, resolution): """ name = 'decomp_test' self.resolution = resolution - subdir = '{}/{}'.format(resolution, name) + subdir = f'{resolution}/{forcing}/{name}' super().__init__(test_group=test_group, name=name, subdir=subdir) diff --git a/compass/ocean/tests/turbulence_closure/default/__init__.py b/compass/ocean/tests/turbulence_closure/default/__init__.py index 6238fe7500..67fc5015aa 100644 --- a/compass/ocean/tests/turbulence_closure/default/__init__.py +++ b/compass/ocean/tests/turbulence_closure/default/__init__.py @@ -17,7 +17,7 @@ class Default(TestCase): The resolution of the test case """ - def __init__(self, test_group, resolution): + def __init__(self, test_group, resolution, forcing='cooling'): """ Create the test case @@ -28,23 +28,32 @@ def __init__(self, test_group, resolution): resolution : str The resolution of the test case + + forcing: str + The forcing applied to the test case """ name = 'default' self.resolution = resolution - subdir = '{}/{}'.format(resolution, name) + self.forcing = forcing + subdir = f'{resolution}/{forcing}/{name}' super().__init__(test_group=test_group, name=name, subdir=subdir) + if resolution == '1m': + ntasks = 128 + else: + ntasks = 4 + self.add_step( InitialState(test_case=self, resolution=resolution)) self.add_step( - Forward(test_case=self, ntasks=4, openmp_threads=1, resolution=resolution)) - self.add_step(Viz(test_case=self)) + Forward(test_case=self, ntasks=ntasks, openmp_threads=1, resolution=resolution)) + self.add_step(Viz(test_case=self, resolution=resolution, forcing=forcing)) def configure(self): """ Modify the configuration options for this test case. """ - turbulence_closure.configure(self.resolution, self.config) + turbulence_closure.configure(self.resolution, self.forcing, self.config) # no run() is needed because we're doing the default: running all steps diff --git a/compass/ocean/tests/turbulence_closure/restart_test/__init__.py b/compass/ocean/tests/turbulence_closure/restart_test/__init__.py index 6e075545e9..c3e434b1dc 100644 --- a/compass/ocean/tests/turbulence_closure/restart_test/__init__.py +++ b/compass/ocean/tests/turbulence_closure/restart_test/__init__.py @@ -17,7 +17,7 @@ class RestartTest(TestCase): The resolution of the test case """ - def __init__(self, test_group, resolution): + def __init__(self, test_group, resolution, forcing='cooling'): """ Create the test case @@ -31,7 +31,7 @@ def __init__(self, test_group, resolution): """ name = 'restart_test' self.resolution = resolution - subdir = '{}/{}'.format(resolution, name) + subdir = f'{resolution}/{forcing}/{name}' super().__init__(test_group=test_group, name=name, subdir=subdir) From 2112c55d477f285e9d61918c0bd070ea417d8a3f Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 22 Dec 2022 14:16:59 -0600 Subject: [PATCH 12/21] Add comparison to palm output --- .../turbulence_closure/default/__init__.py | 6 ++- compass/ocean/tests/turbulence_closure/viz.py | 42 ++++++++++++++++--- 2 files changed, 41 insertions(+), 7 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/default/__init__.py b/compass/ocean/tests/turbulence_closure/default/__init__.py index 67fc5015aa..7bc580c312 100644 --- a/compass/ocean/tests/turbulence_closure/default/__init__.py +++ b/compass/ocean/tests/turbulence_closure/default/__init__.py @@ -39,16 +39,18 @@ def __init__(self, test_group, resolution, forcing='cooling'): super().__init__(test_group=test_group, name=name, subdir=subdir) - if resolution == '1m': + if resolution == '1m' or resolution == '2m': ntasks = 128 + plot_comparison=True else: ntasks = 4 + plot_comparison=False self.add_step( InitialState(test_case=self, resolution=resolution)) self.add_step( Forward(test_case=self, ntasks=ntasks, openmp_threads=1, resolution=resolution)) - self.add_step(Viz(test_case=self, resolution=resolution, forcing=forcing)) + self.add_step(Viz(test_case=self, resolution=resolution, forcing=forcing, do_comparison=plot_comparison)) def configure(self): """ diff --git a/compass/ocean/tests/turbulence_closure/viz.py b/compass/ocean/tests/turbulence_closure/viz.py index 5b0a37d4b5..e550db68d2 100644 --- a/compass/ocean/tests/turbulence_closure/viz.py +++ b/compass/ocean/tests/turbulence_closure/viz.py @@ -10,7 +10,7 @@ class Viz(Step): """ A step for visualizing a cross-section through the domain """ - def __init__(self, test_case): + def __init__(self, test_case, resolution, forcing, do_comparison=False): """ Create the step @@ -23,10 +23,27 @@ def __init__(self, test_case): self.add_input_file(filename='output.nc', target='../forward/output.nc') - self.add_output_file('uNormal_depth_section_t0.png') - self.add_output_file('pt_depth_section_t0.png') - self.add_output_file('sa_depth_section_t0.png') - self.add_output_file('layerThickness_depth_section_t0.png') + + if resolution == '1m': + suffix = 'g128_l128' + + if do_comparison: + if forcing == 'cooling': + forcing_name = 'c02' + elif forcing == 'evaporation': + forcing_name = 'e04' + else: + #TODO change to log + print('Comparison simulation not available for this configuration') + do_comparison = False + + if do_comparison: + filename = f'case_{forcing_name}_{suffix}.nc' + print(f'Compare to {filename}') + self.add_input_file(filename='palm.nc', target=filename, + database='turbulence_closure') + + self.do_comparison = do_comparison def run(self): """ @@ -36,6 +53,8 @@ def run(self): ds = xarray.open_dataset('output.nc') ds = ds.sortby('yEdge') + if self.do_comparison: + dsPalm = xarray.open_dataset('palm.nc') dsInit = xarray.open_dataset('../forward/init.nc') ds0 = ds.isel(Time=0) figsize = [6.4, 4.8] @@ -97,6 +116,10 @@ def run(self): normalVelocity = ds1.normalVelocity normalVelocity_xmesh = normalVelocity[edgeMask_x, :] + velocityZonal = ds1.velocityZonal + velocityZonal_z = velocityZonal.mean(dim='nCells') + velocityMeridional = ds1.velocityMeridional + velocityMeridional_z = velocityMeridional.mean(dim='nCells') # Import cell quantities layerThickness = ds1.layerThickness @@ -178,6 +201,15 @@ def run(self): bbox_inches='tight', dpi=200) plt.close() + plt.figure(figsize=figsize, dpi=100) + plt.plot(velocityZonal_z.values, zMid) + plt.plot(velocityMeridional_z.values, zMid, '--') + plt.xlabel('u,v (m/s)') + plt.ylabel('z (m)') + plt.savefig('uv_depth_t{}.png'.format(j), + bbox_inches='tight', dpi=200) + plt.close() + plt.figure(figsize=figsize, dpi=100) cmax = numpy.max(numpy.abs(w_x.values)) plt.pcolormesh(numpy.divide(yCells_mesh[:-1,:], 1e3), From a82f797967cdcd12a727bddf25213f4281ae7829 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 24 Jan 2023 14:08:26 -0600 Subject: [PATCH 13/21] Add 2m case --- compass/ocean/tests/turbulence_closure/__init__.py | 7 ++++++- compass/ocean/tests/turbulence_closure/forward.py | 7 +++++-- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/__init__.py b/compass/ocean/tests/turbulence_closure/__init__.py index 957edca7b4..2fa36b37dd 100644 --- a/compass/ocean/tests/turbulence_closure/__init__.py +++ b/compass/ocean/tests/turbulence_closure/__init__.py @@ -20,7 +20,7 @@ def __init__(self, mpas_core): DecompTest(test_group=self, resolution=resolution)) self.add_test_case( RestartTest(test_group=self, resolution=resolution)) - for resolution in ['1m', '10km']: + for resolution in ['1m', '2m', '10km']: for forcing in ['cooling', 'evaporation']: self.add_test_case( Default(test_group=self, resolution=resolution, forcing=forcing)) @@ -41,11 +41,16 @@ def configure(resolution, forcing, config): res_params = {'10km': {'nx': 16, 'ny': 50, 'dc': 10e3}, + '2m': {'nx': 150, + 'ny': 150, + 'dc': 2}, '1m': {'nx': 128, 'ny': 128, 'dc': 1}} vert_params = {'10km': {'vert_levels': 20, 'bottom_depth': 1e3}, + '2m': {'vert_levels': 50, + 'bottom_depth': 100.0}, '1m': {'vert_levels': 128, 'bottom_depth': 128.0}} diff --git a/compass/ocean/tests/turbulence_closure/forward.py b/compass/ocean/tests/turbulence_closure/forward.py index ca98307692..9e4210cc6e 100644 --- a/compass/ocean/tests/turbulence_closure/forward.py +++ b/compass/ocean/tests/turbulence_closure/forward.py @@ -78,9 +78,12 @@ def run(self): resolution = self.resolution if resolution == '10km': self.update_namelist_at_runtime({'config_dt': - "'0000_00:01:00'"}) + "'0000_00:00:01'"}) + elif resolution == '2m': + self.update_namelist_at_runtime({'config_dt': + "'0000_00:00:00.1'"}) elif resolution == '1m': self.update_namelist_at_runtime({'config_dt': - "'0000_00:00:00.5'"}) + "'0000_00:00:00.1'"}) run_model(self) From e238374724ea0e665da3cec303845533ca76e058 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 24 Jan 2023 14:08:49 -0600 Subject: [PATCH 14/21] Match cfg, namelist options with successful run --- .../tests/turbulence_closure/namelist.forward | 22 ++++++++++++++++--- .../turbulence_closure/turbulence_closure.cfg | 4 ++-- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/namelist.forward b/compass/ocean/tests/turbulence_closure/namelist.forward index a003b266ab..6209aa0d6b 100644 --- a/compass/ocean/tests/turbulence_closure/namelist.forward +++ b/compass/ocean/tests/turbulence_closure/namelist.forward @@ -1,7 +1,8 @@ config_write_output_on_startup = .false. config_time_integrator = 'split_implicit' -config_run_duration = '0001_00:00:00' -config_implicit_bottom_drag_coeff = 1.0e-2 +config_run_duration = '0004_00:00:00' +config_dt = '00:00:00.1' +config_use_implicit_bottom_drag = .false. config_zonal_ssh_grad = 0.0 config_meridional_ssh_grad = 0.0 config_aust_scale_factor = 0.70710678 @@ -10,6 +11,7 @@ config_mom_div2 = .00002 config_use_mom_div4 = .false. config_mom_div4 = 1.0e12 config_enable_nonhydrostatic_mode = .true. +config_triskCV = .false. config_nonhydrostatic_preconditioner = 'bjacobi' config_nonhydrostatic_solver_type = 'cg' config_nonhydrostatic_solve_surface_boundary_condition = 'pressureTopGradientBottom' @@ -18,7 +20,21 @@ config_use_two_equation_turbulence_model = .true. config_two_equation_model_choice = 'omega' config_use_activeTracers_surface_bulk_forcing = .true. config_LES_mode = .true. +config_vertMom_del2 = 2e-6 +config_tracer_del2 = 2e-6 config_mom_del2 = .true. config_mom_del2 = 2e-6 config_use_cvmix = .false. -config_btr_si_partition_match_mode = .true. +config_cvmix_background_diffusion = 1.0e-2 +config_cvmix_background_viscosity = 1.0e-2 +config_btr_si_partition_match_mode = .false. +config_two_equation_length_min = 0.1 +config_two_equation_tke_min = 1e-10 +config_two_equation_psi_min = 1e-10 +config_tke_ini = 1.0e-30 +config_thickness_flux_type = 'upwind' +config_LES_noise_enable = .true. +config_LES_noise_max_time = 150 +config_LES_noise_min_level = 5 +config_LES_noise_max_level = 45 +config_LES_noise_perturbation_magnitude = 1e-4 diff --git a/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg b/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg index f29a1bc4f8..b1ad3076ad 100644 --- a/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg +++ b/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg @@ -47,7 +47,7 @@ interior_temperature_gradient = 0.1 interior_salinity_gradient = 0 # Disturbance amplitude (added to initiate convection faster; ^oC) -disturbance_amplitude = 0.01 +disturbance_amplitude = 0.005 # Surface heat flux (W/m^2) surface_heat_flux = -100 @@ -62,4 +62,4 @@ wind_stress_zonal = 0 wind_stress_meridional = 0 # Coriolis parameter for entire domain. -coriolis_parameter = -1.2e-4 +coriolis_parameter = 1.0e-4 From 2a1ee3aa34d9c55aa94e141b03c882cca23d8d41 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 31 Jan 2023 18:16:11 -0600 Subject: [PATCH 15/21] Clean-up viz --- compass/ocean/tests/turbulence_closure/viz.py | 335 ++++++++---------- 1 file changed, 141 insertions(+), 194 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/viz.py b/compass/ocean/tests/turbulence_closure/viz.py index e550db68d2..2f34651c5a 100644 --- a/compass/ocean/tests/turbulence_closure/viz.py +++ b/compass/ocean/tests/turbulence_closure/viz.py @@ -1,9 +1,8 @@ -import xarray -import numpy +import xarray as xr +import numpy as np import matplotlib.pyplot as plt from compass.step import Step -import cmocean class Viz(Step): @@ -23,9 +22,9 @@ def __init__(self, test_case, resolution, forcing, do_comparison=False): self.add_input_file(filename='output.nc', target='../forward/output.nc') + self.add_input_file(filename='mesh.nc', + target='../initial_state/culled_mesh.nc') - if resolution == '1m': - suffix = 'g128_l128' if do_comparison: if forcing == 'cooling': @@ -33,13 +32,12 @@ def __init__(self, test_case, resolution, forcing, do_comparison=False): elif forcing == 'evaporation': forcing_name = 'e04' else: - #TODO change to log - print('Comparison simulation not available for this configuration') do_comparison = False if do_comparison: + # Compare all cases with 1m resolution PALM output + suffix = 'g128_l128' filename = f'case_{forcing_name}_{suffix}.nc' - print(f'Compare to {filename}') self.add_input_file(filename='palm.nc', target=filename, database='turbulence_closure') @@ -50,202 +48,151 @@ def run(self): Run this step of the test case """ - ds = xarray.open_dataset('output.nc') - ds = ds.sortby('yEdge') + figsize = [6.4, 4.8] + markersize = 5 + dsInit = xr.open_dataset('../forward/init.nc') + ds = xr.open_dataset('output.nc') if self.do_comparison: - dsPalm = xarray.open_dataset('palm.nc') - dsInit = xarray.open_dataset('../forward/init.nc') - ds0 = ds.isel(Time=0) - figsize = [6.4, 4.8] - markersize = 20 + ds_palm = xr.open_dataset('palm.nc') + if 'Time' not in ds.dims: print('Dataset missing time dimension') return - nSteps = ds.sizes['Time'] # number of timesteps - tend = nSteps - 1 - nCells = dsInit.sizes['nCells'] - nEdges = dsInit.sizes['nEdges'] - nVertLevels = dsInit.sizes['nVertLevels'] + nt = ds.sizes['Time'] # number of timesteps + tidx = nt - 1 + seconds_per_day = 24.0 * 3600.0 + if 'daysSinceStartOfSim' in ds.keys(): + time = ds.daysSinceStartOfSim + else: + # This routine is not generic but should not be used as + # daysSinceStartOfSim is included in the streams file + time0 = 2.0 + (7.0/24.0) + dt = 0.5/24.0 + time = np.linspace(time0 + dt, time0 + nt*dt, num=nt) + if self.do_comparison: + time_palm = ds_palm.time + time_palm_day = (time_palm.astype('float64') / + (seconds_per_day * 1e9)) + tidx_palm = np.argmin(np.abs(np.subtract( + time_palm_day.values, time[tidx]))) + + ds = ds.isel(Time=tidx) + ds_palm = ds_palm.isel(time=tidx_palm) + + if 'yEdge' not in ds.keys(): + ds['yEdge'] = dsInit.yEdge + ds = ds.sortby('yEdge') + + # Get mesh variables xEdge = dsInit.xEdge yEdge = dsInit.yEdge xCell = dsInit.xCell yCell = dsInit.yCell - - xEdge_mid = numpy.median(xEdge) - edgeMask_x = numpy.equal(xEdge, xEdge_mid) - - # Solve for lateral boundaries of uNormal at cell centers for - # x-section - cellsOnEdge = dsInit.cellsOnEdge - cellsOnEdge_x = cellsOnEdge[edgeMask_x, :] - yEdges = numpy.zeros((len(cellsOnEdge_x)+1)) - for i in range(len(cellsOnEdge_x)): - if cellsOnEdge[i, 1] == 0: - yEdges[i] = yCell[cellsOnEdge_x[i, 0] - 1] - yEdges[i+1] = yCell[cellsOnEdge_x[i, 0] - 1] - elif cellsOnEdge[i, 1] == 0: - yEdges[i] = yCell[cellsOnEdge_x[i, 1] - 1] - yEdges[i+1] = yCell[cellsOnEdge_x[i, 1] - 1] - else: - yEdges[i] = min(yCell[cellsOnEdge_x[i, 0] - 1], - yCell[cellsOnEdge_x[i, 1] - 1]) - yEdges[i+1] = max(yCell[cellsOnEdge_x[i, 0] - 1], - yCell[cellsOnEdge_x[i, 1] - 1]) - - # Prep variables for cell quantities - cellIndex = numpy.subtract(cellsOnEdge_x[1:, 0], 1) - yEdge_x = yEdge[edgeMask_x] - - # prep all variables for uNormal plot - zInterface = dsInit.refInterfaces zMid = dsInit.refZMid - zInterfaces_edge_mesh, yEdges_mesh = numpy.meshgrid(zInterface, - yEdges) - zInterfaces_cell_mesh, yCells_mesh = numpy.meshgrid(zInterface, - yEdge_x) - temperature0 = ds0.temperature - temperature0_z = temperature0.mean(dim='nCells') - - for j in [tend]: - ds1 = ds.isel(Time=j) - - normalVelocity = ds1.normalVelocity - normalVelocity_xmesh = normalVelocity[edgeMask_x, :] - - velocityZonal = ds1.velocityZonal - velocityZonal_z = velocityZonal.mean(dim='nCells') - velocityMeridional = ds1.velocityMeridional - velocityMeridional_z = velocityMeridional.mean(dim='nCells') - - # Import cell quantities - layerThickness = ds1.layerThickness - layerThickness_x = layerThickness[cellIndex, :] - temperature = ds1.temperature - temperature_z = temperature.mean(dim='nCells') - temperature_x = temperature[cellIndex, :] - salinity = ds1.salinity - salinity_z = salinity.mean(dim='nCells') - salinity_x = salinity[cellIndex, :] - w = ds1.vertVelocityTop - w_x = w[cellIndex, :] - w_zTop = w[:, 0] - - # Figures - plt.figure(figsize=figsize, dpi=100) - cmax = numpy.max(numpy.abs(normalVelocity_xmesh.values)) - plt.pcolormesh(numpy.divide(yEdges_mesh, 1e3), - zInterfaces_edge_mesh, - normalVelocity_xmesh.values, - cmap='cmo.balance', vmin=-1.*cmax, vmax=cmax) - cbar = plt.colorbar() - cbar.ax.set_title('uNormal (m/s)') - plt.savefig('uNormal_depth_section_t{}.png'.format(j), - bbox_inches='tight', dpi=200) - plt.close() - - # ------------------------------------------------------------------ - # Plot cell-centered variables - # ------------------------------------------------------------------ - # Figures - plt.figure(figsize=figsize, dpi=100) - plt.plot(temperature_z.values, zMid) - plt.xlabel('PT (C)') - plt.ylabel('z (m)') - plt.savefig('pt_depth_t{}.png'.format(j), - bbox_inches='tight', dpi=200) - plt.close() - - plt.figure(figsize=figsize, dpi=100) - plt.pcolormesh(numpy.divide(yCells_mesh, 1e3), - zInterfaces_cell_mesh, - temperature_x.values, cmap='cmo.thermal') - plt.xlabel('y (km)') - plt.ylabel('z (m)') - cbar = plt.colorbar() - cbar.ax.set_title('PT (C)') - plt.savefig('pt_depth_section_t{}.png'.format(j), - bbox_inches='tight', dpi=200) - plt.close() - - plt.figure(figsize=figsize, dpi=100) - plt.plot(numpy.subtract(temperature_z.values, - temperature0_z.values), - zMid) - plt.xlabel('delta PT (C)') - plt.ylabel('z (m)') - plt.savefig('dpt_depth_t{}.png'.format(j), - bbox_inches='tight', dpi=200) - plt.close() - - plt.figure(figsize=figsize, dpi=100) - plt.plot(salinity_z.values, zMid) - plt.xlabel('SA (g/kg)') - plt.ylabel('z (m)') - plt.savefig('sa_depth_t{}.png'.format(j), - bbox_inches='tight', dpi=200) - plt.close() - - plt.figure(figsize=figsize, dpi=100) - plt.pcolormesh(numpy.divide(yCells_mesh, 1e3), - zInterfaces_cell_mesh, - salinity_x.values, cmap='cmo.haline') - plt.xlabel('y (km)') - plt.ylabel('z (m)') - cbar = plt.colorbar() - cbar.ax.set_title('SA (g/kg)') - plt.savefig('sa_depth_section_t{}.png'.format(j), - bbox_inches='tight', dpi=200) - plt.close() - - plt.figure(figsize=figsize, dpi=100) - plt.plot(velocityZonal_z.values, zMid) - plt.plot(velocityMeridional_z.values, zMid, '--') - plt.xlabel('u,v (m/s)') - plt.ylabel('z (m)') - plt.savefig('uv_depth_t{}.png'.format(j), - bbox_inches='tight', dpi=200) - plt.close() - - plt.figure(figsize=figsize, dpi=100) - cmax = numpy.max(numpy.abs(w_x.values)) - plt.pcolormesh(numpy.divide(yCells_mesh[:-1,:], 1e3), - zInterfaces_cell_mesh[:-1,:], - w_x.values, cmap='cmo.balance', - vmin=-1.*cmax, vmax=cmax) - plt.xlabel('y (km)') - plt.ylabel('z (m)') - cbar = plt.colorbar() - cbar.ax.set_title('w (m/s)') - plt.savefig('w_depth_section_t{}.png'.format(j), - bbox_inches='tight', dpi=200) - plt.close() - - plt.figure(figsize=figsize, dpi=100) - cmax = numpy.max(numpy.abs(w_zTop.values)) - plt.scatter(numpy.divide(xCell, 1e3), - numpy.divide(yCell, 1e3), - s=5, c=w_zTop.values, - cmap='cmo.balance', vmin=-1.*cmax, vmax=cmax) - plt.xlabel('x (km)') - plt.ylabel('y (km)') - cbar = plt.colorbar() - cbar.ax.set_title('w (m/s)') - plt.savefig('w_top_section_t{}.png'.format(j), - bbox_inches='tight', dpi=200) - plt.close() - - plt.figure(figsize=figsize, dpi=100) - plt.pcolormesh(numpy.divide(yCells_mesh, 1e3), - zInterfaces_cell_mesh, - layerThickness_x.values, cmap='viridis') - plt.xlabel('y (km)') - plt.ylabel('z (m)') - cbar = plt.colorbar() - cbar.ax.set_title('h (m)') - plt.savefig('layerThickness_depth_section_t{}.png'.format(j), - bbox_inches='tight', dpi=200) - plt.close() + # Import cell quantities + temperature = ds.temperature + temperature_z = temperature.mean(dim='nCells') + salinity = ds.salinity + salinity_z = salinity.mean(dim='nCells') + w = ds.verticalVelocity + w_zTop = w[:, 0] + + velocityZonal = ds.velocityZonal + velocityZonal_z = velocityZonal.mean(dim='nCells') + velocityMeridional = ds.velocityMeridional + velocityMeridional_z = velocityMeridional.mean(dim='nCells') + buoyancyProduction = ds.buoyancyProduction + buoyancyProduction_z = buoyancyProduction.mean(dim='nCells') + wpt = ds.temperatureVerticalAdvectionTendency + wpt_z = wpt.mean(dim='nCells') + + if self.do_comparison: + alpha_T = ds_palm.alpha_T + if 'beta_S' in ds_palm.keys(): + beta_S = ds_palm.beta_S + else: + beta_S = 7.8e-4 + pt_palm = ds_palm.pt - 273.15 + sa_palm = ds_palm.sa + u_palm = ds_palm.u + v_palm = ds_palm.v + wpt_palm = np.add(ds_palm['w*pt*'].values, ds_palm['w"pt"'].values) + wsa_palm = np.add(ds_palm['w*sa*'].values, ds_palm['w"sa"'].values) + temp1 = np.multiply(alpha_T, wpt_palm) + temp2 = np.multiply(beta_S, wsa_palm) + buoyancyProduction_palm = np.multiply(9.81, + np.subtract(temp1, temp2)) + zu_palm = ds_palm.zu + z_palm = ds_palm.zprho + + # Figures + + plt.figure(figsize=figsize, dpi=100) + plt.plot(temperature_z.values, zMid, 'k-') + if self.do_comparison: + plt.plot(pt_palm.values, z_palm.values, 'b-') + plt.xlabel('PT (C)') + plt.ylabel('z (m)') + plt.savefig(f'pt_depth_t{int(time[tidx]*24.0)}h.png', + bbox_inches='tight', dpi=200) + plt.close() + + plt.figure(figsize=figsize, dpi=100) + plt.plot(wpt_z, zMid, 'k-') + if self.do_comparison: + plt.plot(wpt_palm, z_palm.values, 'b-') + plt.xlabel('wpt (C m/s)') + plt.ylabel('z (m)') + plt.savefig(f'wpt_depth_t{int(time[tidx]*24.0)}h.png', + bbox_inches='tight', dpi=200) + plt.close() + + plt.figure(figsize=figsize, dpi=100) + plt.plot(buoyancyProduction_z, zMid, 'k-') + if self.do_comparison: + plt.plot(buoyancyProduction_palm, z_palm.values, 'b-') + plt.xlabel('bouyancy production') + plt.ylabel('z (m)') + plt.savefig(f'buoy_depth_t{int(time[tidx]*24.0)}h.png', + bbox_inches='tight', dpi=200) + plt.close() + + plt.figure(figsize=figsize, dpi=100) + plt.plot(salinity_z.values, zMid, 'k-') + if self.do_comparison: + plt.plot(sa_palm.values, z_palm.values, 'b-') + plt.xlabel('SA (g/kg)') + plt.ylabel('z (m)') + plt.savefig(f'sa_depth_t{int(time[tidx]*24.0)}h.png', + bbox_inches='tight', dpi=200) + plt.close() + + plt.figure(figsize=figsize, dpi=100) + plt.plot(velocityZonal_z.values, zMid, 'k-') + plt.plot(velocityMeridional_z.values, zMid, 'k--') + if self.do_comparison: + plt.plot(u_palm.values, zu_palm.values, 'b-') + plt.plot(v_palm.values, zu_palm.values, 'b--') + plt.xlabel('u,v (m/s)') + plt.ylabel('z (m)') + plt.savefig(f'uv_depth_t{int(time[tidx]*24.0)}h.png', + bbox_inches='tight', dpi=200) + plt.close() + + plt.figure(figsize=figsize, dpi=100) + cmax = np.max(np.abs(w_zTop.values)) + plt.scatter(np.divide(xCell, 1e3), + np.divide(yCell, 1e3), + s=markersize, c=w_zTop.values, + cmap='cmo.balance', vmin=-1.*cmax, vmax=cmax) + plt.xlabel('x (km)') + plt.ylabel('y (km)') + cbar = plt.colorbar() + cbar.ax.set_title('w (m/s)') + plt.savefig(f'w_top_section_t{int(time[tidx]*24.0)}h.png', + bbox_inches='tight', dpi=200) + plt.close() From 4319588385fd0f0a9de6f62770e1578c5b6db756 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 2 Feb 2023 15:30:34 -0600 Subject: [PATCH 16/21] Differentiate namelists --- .../ocean/tests/turbulence_closure/forward.py | 6 +++++ .../tests/turbulence_closure/namelist.forward | 24 +------------------ .../turbulence_closure/namelist.les.forward | 13 ++++++++++ .../namelist.nonhydro.forward | 8 +++++++ 4 files changed, 28 insertions(+), 23 deletions(-) create mode 100644 compass/ocean/tests/turbulence_closure/namelist.les.forward create mode 100644 compass/ocean/tests/turbulence_closure/namelist.nonhydro.forward diff --git a/compass/ocean/tests/turbulence_closure/forward.py b/compass/ocean/tests/turbulence_closure/forward.py index 9e4210cc6e..2bbb838a5a 100644 --- a/compass/ocean/tests/turbulence_closure/forward.py +++ b/compass/ocean/tests/turbulence_closure/forward.py @@ -50,6 +50,12 @@ def __init__(self, test_case, resolution, name='forward', subdir=None, ntasks=ntasks, min_tasks=min_tasks, openmp_threads=openmp_threads) self.add_namelist_file('compass.ocean.tests.turbulence_closure', 'namelist.forward') + if resolution < 500: + self.add_namelist_file('compass.ocean.tests.turbulence_closure', + 'namelist.les.forward') + if resolution < 100: + self.add_namelist_file('compass.ocean.tests.turbulence_closure', + 'namelist.nonhydro.forward') # make sure output is double precision self.add_streams_file('compass.ocean.streams', 'streams.output') diff --git a/compass/ocean/tests/turbulence_closure/namelist.forward b/compass/ocean/tests/turbulence_closure/namelist.forward index 6209aa0d6b..d1b9c76235 100644 --- a/compass/ocean/tests/turbulence_closure/namelist.forward +++ b/compass/ocean/tests/turbulence_closure/namelist.forward @@ -1,40 +1,18 @@ config_write_output_on_startup = .false. config_time_integrator = 'split_implicit' +config_btr_si_partition_match_mode = .false. config_run_duration = '0004_00:00:00' -config_dt = '00:00:00.1' config_use_implicit_bottom_drag = .false. config_zonal_ssh_grad = 0.0 config_meridional_ssh_grad = 0.0 -config_aust_scale_factor = 0.70710678 config_use_mom_div2 = .true. config_mom_div2 = .00002 config_use_mom_div4 = .false. config_mom_div4 = 1.0e12 -config_enable_nonhydrostatic_mode = .true. -config_triskCV = .false. -config_nonhydrostatic_preconditioner = 'bjacobi' -config_nonhydrostatic_solver_type = 'cg' -config_nonhydrostatic_solve_surface_boundary_condition = 'pressureTopGradientBottom' -config_use_constant_forced_pgf = .false. -config_use_two_equation_turbulence_model = .true. -config_two_equation_model_choice = 'omega' config_use_activeTracers_surface_bulk_forcing = .true. -config_LES_mode = .true. config_vertMom_del2 = 2e-6 config_tracer_del2 = 2e-6 config_mom_del2 = .true. config_mom_del2 = 2e-6 -config_use_cvmix = .false. config_cvmix_background_diffusion = 1.0e-2 config_cvmix_background_viscosity = 1.0e-2 -config_btr_si_partition_match_mode = .false. -config_two_equation_length_min = 0.1 -config_two_equation_tke_min = 1e-10 -config_two_equation_psi_min = 1e-10 -config_tke_ini = 1.0e-30 -config_thickness_flux_type = 'upwind' -config_LES_noise_enable = .true. -config_LES_noise_max_time = 150 -config_LES_noise_min_level = 5 -config_LES_noise_max_level = 45 -config_LES_noise_perturbation_magnitude = 1e-4 diff --git a/compass/ocean/tests/turbulence_closure/namelist.les.forward b/compass/ocean/tests/turbulence_closure/namelist.les.forward new file mode 100644 index 0000000000..00ea3f1591 --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/namelist.les.forward @@ -0,0 +1,13 @@ +config_LES_mode = .true. +config_use_two_equation_turbulence_model = .true. +config_two_equation_model_choice = 'omega' +config_use_cvmix = .false. +config_two_equation_length_min = 0.1 +config_two_equation_tke_min = 1e-10 +config_two_equation_psi_min = 1e-10 +config_tke_ini = 1.0e-30 +config_LES_noise_enable = .true. +config_LES_noise_max_time = 150 +config_LES_noise_min_level = 5 +config_LES_noise_max_level = 45 +config_LES_noise_perturbation_magnitude = 1e-4 diff --git a/compass/ocean/tests/turbulence_closure/namelist.nonhydro.forward b/compass/ocean/tests/turbulence_closure/namelist.nonhydro.forward new file mode 100644 index 0000000000..281f2de6ce --- /dev/null +++ b/compass/ocean/tests/turbulence_closure/namelist.nonhydro.forward @@ -0,0 +1,8 @@ +config_enable_nonhydrostatic_mode = .true. +config_triskCV = .false. +config_nonhydrostatic_preconditioner = 'bjacobi' +config_nonhydrostatic_solver_type = 'cg' +config_nonhydrostatic_solve_surface_boundary_condition = 'pressureTopGradientBottom' +config_use_constant_forced_pgf = .false. +config_thickness_flux_type = 'upwind' +config_aust_scale_factor = 0.70710678 From 511ff3aef5f35e8d287714d965266fc09416bf81 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 2 Feb 2023 15:46:12 -0600 Subject: [PATCH 17/21] Replace resolution str with float --- .../tests/turbulence_closure/__init__.py | 59 ++++++++++--------- .../decomp_test/__init__.py | 14 +++-- .../turbulence_closure/default/__init__.py | 16 +++-- .../restart_test/__init__.py | 14 +++-- 4 files changed, 59 insertions(+), 44 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/__init__.py b/compass/ocean/tests/turbulence_closure/__init__.py index 2fa36b37dd..41ef11321f 100644 --- a/compass/ocean/tests/turbulence_closure/__init__.py +++ b/compass/ocean/tests/turbulence_closure/__init__.py @@ -15,12 +15,12 @@ def __init__(self, mpas_core): """ super().__init__(mpas_core=mpas_core, name='turbulence_closure') - for resolution in ['10km']: + for resolution in [1e4]: self.add_test_case( DecompTest(test_group=self, resolution=resolution)) self.add_test_case( RestartTest(test_group=self, resolution=resolution)) - for resolution in ['1m', '2m', '10km']: + for resolution in [1, 2, 1e4]: for forcing in ['cooling', 'evaporation']: self.add_test_case( Default(test_group=self, resolution=resolution, forcing=forcing)) @@ -32,38 +32,41 @@ def configure(resolution, forcing, config): Parameters ---------- - resolution : str - The resolution of the test case + resolution : float + The resolution of the test case in meters config : configparser.ConfigParser Configuration options for this test case """ - res_params = {'10km': {'nx': 16, - 'ny': 50, - 'dc': 10e3}, - '2m': {'nx': 150, - 'ny': 150, - 'dc': 2}, - '1m': {'nx': 128, - 'ny': 128, - 'dc': 1}} - vert_params = {'10km': {'vert_levels': 20, - 'bottom_depth': 1e3}, - '2m': {'vert_levels': 50, - 'bottom_depth': 100.0}, - '1m': {'vert_levels': 128, - 'bottom_depth': 128.0}} + # The resolution parameters are different for different resolutions + # to match existing simulations + if resolution > 1e3: + nx = 16 + ny = 50 + vert_levels = 20 + bottom_depth = 1e3 + elif resolution <= 1e3 and resolution > 5: + nx = 50 + ny = 50 + vert_levels = 50 + bottom_depth = 100.0 + elif resolution <= 5 and resolution > 1: + nx = 150 + ny = 150 + vert_levels = 50 + bottom_depth = 100.0 + elif resolution <= 1: + nx = 128 + ny = 128 + vert_levels = 128 + bottom_depth = 128.0 - if resolution not in res_params: - raise ValueError(f'Unsupported resolution {resolution}. ' - f'Supported values are: {list(res_params)}') - res_params = res_params[resolution] - for param in res_params: - config.set('turbulence_closure', param, f'{res_params[param]}') - vert_params = vert_params[resolution] - for param in vert_params: - config.set('vertical_grid', param, f'{vert_params[param]}') + config.set('turbulence_closure', 'nx', f'{nx}') + config.set('turbulence_closure', 'ny', f'{ny}') + config.set('turbulence_closure', 'dc', f'{resolution}') + config.set('vertical_grid', 'vert_levels', f'{vert_levels}') + config.set('vertical_grid', 'bottom_depth', f'{bottom_depth}') if forcing == 'cooling': config.set('turbulence_closure', 'surface_heat_flux', '-100') diff --git a/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py b/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py index 8d26a51148..e0f7ae809b 100644 --- a/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py +++ b/compass/ocean/tests/turbulence_closure/decomp_test/__init__.py @@ -12,8 +12,8 @@ class DecompTest(TestCase): Attributes ---------- - resolution : str - The resolution of the test case + resolution : float + The resolution of the test case in meters """ def __init__(self, test_group, resolution, forcing='cooling'): @@ -25,12 +25,16 @@ def __init__(self, test_group, resolution, forcing='cooling'): test_group : compass.ocean.tests.turbulence_closure.TurbulenceClosure The test group that this test case belongs to - resolution : str - The resolution of the test case + resolution : float + The resolution of the test case in meters """ name = 'decomp_test' self.resolution = resolution - subdir = f'{resolution}/{forcing}/{name}' + if resolution >= 1e3: + res_name = f'{int(resolution/1e3)}km' + else: + res_name = f'{int(resolution)}m' + subdir = f'{res_name}/{forcing}/{name}' super().__init__(test_group=test_group, name=name, subdir=subdir) diff --git a/compass/ocean/tests/turbulence_closure/default/__init__.py b/compass/ocean/tests/turbulence_closure/default/__init__.py index 7bc580c312..7a20b310ce 100644 --- a/compass/ocean/tests/turbulence_closure/default/__init__.py +++ b/compass/ocean/tests/turbulence_closure/default/__init__.py @@ -13,8 +13,8 @@ class Default(TestCase): Attributes ---------- - resolution : str - The resolution of the test case + resolution : float + The resolution of the test case in meters """ def __init__(self, test_group, resolution, forcing='cooling'): @@ -26,8 +26,8 @@ def __init__(self, test_group, resolution, forcing='cooling'): test_group : compass.ocean.tests.turbulence_closure.TurbulenceClosure The test group that this test case belongs to - resolution : str - The resolution of the test case + resolution : float + The resolution of the test case in meters forcing: str The forcing applied to the test case @@ -35,11 +35,15 @@ def __init__(self, test_group, resolution, forcing='cooling'): name = 'default' self.resolution = resolution self.forcing = forcing - subdir = f'{resolution}/{forcing}/{name}' + if resolution >= 1e3: + res_name = f'{int(resolution/1e3)}km' + else: + res_name = f'{int(resolution)}m' + subdir = f'{res_name}/{forcing}/{name}' super().__init__(test_group=test_group, name=name, subdir=subdir) - if resolution == '1m' or resolution == '2m': + if resolution <= 5: ntasks = 128 plot_comparison=True else: diff --git a/compass/ocean/tests/turbulence_closure/restart_test/__init__.py b/compass/ocean/tests/turbulence_closure/restart_test/__init__.py index c3e434b1dc..25b1ec4f46 100644 --- a/compass/ocean/tests/turbulence_closure/restart_test/__init__.py +++ b/compass/ocean/tests/turbulence_closure/restart_test/__init__.py @@ -13,8 +13,8 @@ class RestartTest(TestCase): Attributes ---------- - resolution : str - The resolution of the test case + resolution : float + The resolution of the test case in meters """ def __init__(self, test_group, resolution, forcing='cooling'): @@ -26,12 +26,16 @@ def __init__(self, test_group, resolution, forcing='cooling'): test_group : compass.ocean.tests.baroclinic_channel.BaroclinicChannel The test group that this test case belongs to - resolution : str - The resolution of the test case + resolution : float + The resolution of the test case in meters """ name = 'restart_test' self.resolution = resolution - subdir = f'{resolution}/{forcing}/{name}' + if resolution >= 1e3: + res_name = f'{int(resolution/1e3)}km' + else: + res_name = f'{int(resolution)}m' + subdir = f'{res_name}/{forcing}/{name}' super().__init__(test_group=test_group, name=name, subdir=subdir) From 153f60575dde00735dcae9870312a44105e65e59 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 3 Mar 2023 17:01:48 -0600 Subject: [PATCH 18/21] WIP add turbulence_closure test to docs --- docs/users_guide/ocean/test_groups/index.rst | 3 +- .../ocean/test_groups/turbulence_closure.rst | 108 ++++++++++++++++++ 2 files changed, 110 insertions(+), 1 deletion(-) create mode 100644 docs/users_guide/ocean/test_groups/turbulence_closure.rst diff --git a/docs/users_guide/ocean/test_groups/index.rst b/docs/users_guide/ocean/test_groups/index.rst index 318d0b92ef..b80488f9f7 100644 --- a/docs/users_guide/ocean/test_groups/index.rst +++ b/docs/users_guide/ocean/test_groups/index.rst @@ -29,4 +29,5 @@ coming months. sphere_transport spherical_harmonic_transform tides - ziso \ No newline at end of file + turbulence_closure + ziso diff --git a/docs/users_guide/ocean/test_groups/turbulence_closure.rst b/docs/users_guide/ocean/test_groups/turbulence_closure.rst new file mode 100644 index 0000000000..4afec83860 --- /dev/null +++ b/docs/users_guide/ocean/test_groups/turbulence_closure.rst @@ -0,0 +1,108 @@ +.. _ocean_turbulence_closure: + +turbulence_closure +================== + +The ``turbulence_closure`` test group implements studies of the surface layer +evolution under different kinds of surface forcing for different turbulence +closures. + +The domain is periodic in the horizontal dimensions and no-slip at both +vertical boundaries. + +Initial temperature is 20 degC and salinity is 35 PSU. A linear gradient can be +given over the mixed layer or the interior. Furthermore, a polynomial function +can be used to describe the interior temperature and salinity. The default +parameters are based on a fit to a time average over the whole 1 year +deployment of ITP #118, located in the Beaufort Gyre (western Arctic). + +Variants of the test case are available at 1-m, 2-m and 10-km horizontal +resolution. The O(1)m resolution test cases are designed to be run in LES mode +whereas the O(1)km test cases are designed to be used to test turbulence +closures in "standard" MPAS-Ocean. + +The vertical resolution for the LES test cases is the same as the horizontal +resolution (1m or 2m) whereas the standard resolution test cases have a +vertical resolution of 50m. + +There are currently two surface forcing variants available, cooling and +evaporation. + +All test cases have 2 steps, +``initial_state``, which defines the mesh and initial conditions for the model, +and ``forward``, which performs time integration of the model. The ``default`` +test also has a viz step which performs visualization. + + +config options +-------------- + +namelist options +---------------- +When the horizontal resolution is less than 500m, LES mode is active with the +following config options: + +.. code-block:: cfg + + config_LES_mode = .true. + config_use_two_equation_turbulence_model = .true. + + # TODO Explain what 'omega' means + config_two_equation_model_choice = 'omega' + config_use_cvmix = .false. + + config_two_equation_length_min = 0.1 + config_two_equation_tke_min = 1e-10 + config_two_equation_psi_min = 1e-10 + config_tke_ini = 1.0e-30 + + config_LES_noise_enable = .true. + config_LES_noise_max_time = 150 + config_LES_noise_min_level = 5 + config_LES_noise_max_level = 45 + + # TODO Can you provide any guidance here about what appropriate choices + # might be? + config_LES_noise_perturbation_magnitude = 1e-4 + +When the horizontal resolution is less than 100m, nonhydrostatic mode is +active with the following config options: + +.. code-block:: cfg + + config_enable_nonhydrostatic_mode = .true. + config_triskCV = .false. + config_nonhydrostatic_preconditioner = 'bjacobi' + config_nonhydrostatic_solver_type = 'cg' + config_nonhydrostatic_solve_surface_boundary_condition = 'pressureTopGradientBottom' + config_use_constant_forced_pgf = .false. + config_thickness_flux_type = 'upwind' + config_aust_scale_factor = 0.70710678 + +default +------- + +The default version of the turbulence closure test case is available for +multiple resolutions: ``ocean/baroclinic_channel/1m/default``, +``ocean/baroclinic_channel/2m/default``, and +``ocean/baroclinic_channel/10km/default``. The default simulation time is 4 days. + +decomp_test +----------- + +``ocean/turbulence_closure/10km/decomp_test`` runs a short (15 min) integration +of the model forward in time on 4 (``4proc`` step) and then on 8 processors +(``8proc`` step) to make sure the resulting prognostic variables are +bit-for-bit identical between the two runs. Currently, only 10-km horizontal +resolution is supported. + +restart_test +------------ + +``ocean/turbulence_closure/10km/restart_test`` runs a short (10 min) +integration of the model forward in time (``full_run`` step), saving a restart +file every 5 minutes. Then, a second run (``restart_run`` step) is performed +from the restart file 5 minutes into the simulation and prognostic variables +are compared between the "full" and "restart" runs at minute 10 to make sure +they are bit-for-bit identical. Currently, only 10-km horizontal resolution is +supported. From 96a9311126bdc2ba9fc1a66b78c1d93177a0ce2a Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 13 Mar 2023 10:37:54 -0700 Subject: [PATCH 19/21] Add polynomial functions for Arctic ICs --- .../tests/turbulence_closure/initial_state.py | 163 +++++++++++++----- .../turbulence_closure/turbulence_closure.cfg | 22 +++ 2 files changed, 140 insertions(+), 45 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/initial_state.py b/compass/ocean/tests/turbulence_closure/initial_state.py index 6eb1310a73..571238e90e 100644 --- a/compass/ocean/tests/turbulence_closure/initial_state.py +++ b/compass/ocean/tests/turbulence_closure/initial_state.py @@ -1,5 +1,6 @@ -import xarray -import numpy +import xarray as xr +import numpy as np +import matplotlib.pyplot as plt from compass.ocean.tests.turbulence_closure.boundary_values import add_boundary_arrays from random import random, seed @@ -57,6 +58,8 @@ def run(self): y_nonperiodic = section.getboolean('y_nonperiodic') disturbance_amplitude = section.getfloat('disturbance_amplitude') surface_temperature = section.getfloat('surface_temperature') + interior_temperature_shape = section.get('interior_temperature_shape') + interior_salinity_shape = section.get('interior_salinity_shape') surface_salinity = section.getfloat('surface_salinity') mixed_layer_salinity_gradient = section.getfloat('mixed_layer_salinity_gradient') mixed_layer_temperature_gradient = section.getfloat('mixed_layer_temperature_gradient') @@ -64,6 +67,18 @@ def run(self): mixed_layer_depth_temperature = section.getfloat('mixed_layer_depth_temperature') interior_salinity_gradient = section.getfloat('interior_salinity_gradient') interior_temperature_gradient = section.getfloat('interior_temperature_gradient') + interior_temperature_c0 = section.getfloat('interior_temperature_c0') + interior_temperature_c1 = section.getfloat('interior_temperature_c1') + interior_temperature_c2 = section.getfloat('interior_temperature_c2') + interior_temperature_c3 = section.getfloat('interior_temperature_c3') + interior_temperature_c4 = section.getfloat('interior_temperature_c4') + interior_temperature_c5 = section.getfloat('interior_temperature_c4') + interior_salinity_c0 = section.getfloat('interior_salinity_c0') + interior_salinity_c1 = section.getfloat('interior_salinity_c1') + interior_salinity_c2 = section.getfloat('interior_salinity_c2') + interior_salinity_c3 = section.getfloat('interior_salinity_c3') + interior_salinity_c4 = section.getfloat('interior_salinity_c4') + interior_salinity_c5 = section.getfloat('interior_salinity_c4') surface_heat_flux = section.getfloat('surface_heat_flux') surface_freshwater_flux = section.getfloat('surface_freshwater_flux') wind_stress_zonal = section.getfloat('wind_stress_zonal') @@ -86,64 +101,122 @@ def run(self): xCell = ds.xCell yCell = ds.yCell - ds['bottomDepth'] = bottom_depth * xarray.ones_like(xCell) - ds['ssh'] = xarray.zeros_like(xCell) + ds['bottomDepth'] = bottom_depth * xr.ones_like(xCell) + ds['ssh'] = xr.zeros_like(xCell) init_vertical_coord(config, ds) - dsForcing['windStressZonal'] = wind_stress_zonal * xarray.ones_like(xCell) - dsForcing['windStressMeridional'] = wind_stress_meridional * xarray.ones_like(xCell) - dsForcing['sensibleHeatFlux'] = surface_heat_flux * xarray.ones_like(xCell) - dsForcing['rainFlux'] = surface_freshwater_flux * xarray.ones_like(xCell) + dsForcing['windStressZonal'] = wind_stress_zonal * xr.ones_like(xCell) + dsForcing['windStressMeridional'] = wind_stress_meridional * xr.ones_like(xCell) + dsForcing['sensibleHeatFlux'] = surface_heat_flux * xr.ones_like(xCell) + dsForcing['rainFlux'] = surface_freshwater_flux * xr.ones_like(xCell) write_netcdf(dsForcing, 'init_mode_forcing_data.nc') - normalVelocity = xarray.zeros_like(ds.xEdge) - normalVelocity, _ = xarray.broadcast(normalVelocity, ds.refBottomDepth) + normalVelocity = xr.zeros_like(ds.xEdge) + normalVelocity, _ = xr.broadcast(normalVelocity, ds.refBottomDepth) normalVelocity = normalVelocity.transpose('nEdges', 'nVertLevels') - temperature = xarray.zeros_like(ds.layerThickness) - - salinity = xarray.zeros_like(ds.layerThickness) - - surf_indices = numpy.where(ds.refZMid.values >= -mixed_layer_depth_temperature)[0] - - if len(surf_indices) > 0: - temperature[0,:,surf_indices] = surface_temperature + mixed_layer_temperature_gradient* \ - ds.zMid[0,:,surf_indices].values - - int_indices = numpy.where(ds.refZMid.values < -mixed_layer_depth_temperature)[0] - - if len(int_indices) > 0: - temperature[0,:,int_indices] = surface_temperature - mixed_layer_temperature_gradient* \ - mixed_layer_depth_temperature + interior_temperature_gradient* \ - (ds.zMid[0,:,int_indices] + \ - mixed_layer_depth_temperature) - - temperature[0,:,0] += disturbance_amplitude*2*(numpy.random.rand(len(ds.xCell.values)) - 0.5) - - surf_indices = numpy.where(ds.refZMid.values >= -mixed_layer_depth_salinity)[0] - if len(surf_indices) > 0: - salinity[0,:,surf_indices] = surface_salinity - mixed_layer_salinity_gradient * \ - ds.zMid[0,:,surf_indices] - - int_indices = numpy.where(ds.refZMid.values < -mixed_layer_depth_salinity)[0] - if len(int_indices) > 0: - salinity[0,:,int_indices] = surface_salinity - mixed_layer_salinity_gradient* \ - mixed_layer_depth_salinity + interior_salinity_gradient* \ - (ds.zMid[0,:,int_indices] + \ - mixed_layer_depth_salinity) + temperature = xr.zeros_like(ds.layerThickness) + + salinity = xr.zeros_like(ds.layerThickness) + + zMid = ds.refZMid.values + + temperature_at_mld = (surface_temperature - + mixed_layer_temperature_gradient * mixed_layer_depth_temperature) + if interior_temperature_shape == 'linear': + initial_temperature = np.where( + zMid >= -mixed_layer_depth_temperature, + surface_temperature + mixed_layer_temperature_gradient * zMid, + temperature_at_mld + interior_temperature_gradient * zMid) + elif interior_temperature_shape == 'polynomial': + initial_temperature = np.where( + zMid >= -mixed_layer_depth_temperature, + surface_temperature + mixed_layer_temperature_gradient * zMid, + interior_temperature_c0 + + interior_temperature_c1 * zMid + + interior_temperature_c2 * zMid**2. + + interior_temperature_c3 * zMid**3. + + interior_temperature_c4 * zMid**4. + + interior_temperature_c5 * zMid**5.) + else: + print('interior_temperature_shape is not supported') + + salinity_at_mld = (surface_salinity - + mixed_layer_salinity_gradient * mixed_layer_depth_salinity) + if interior_salinity_shape == 'linear': + initial_salinity = np.where( + zMid >= -mixed_layer_depth_salinity, + surface_salinity + mixed_layer_salinity_gradient * zMid, + salinity_at_mld + interior_salinity_gradient * zMid) + elif interior_salinity_shape == 'polynomial': + initial_salinity = np.where( + zMid >= -mixed_layer_depth_salinity, + surface_salinity + mixed_layer_salinity_gradient * zMid, + interior_salinity_c0 + + interior_salinity_c1 * zMid + + interior_salinity_c2 * zMid**2. + + interior_salinity_c3 * zMid**3. + + interior_salinity_c4 * zMid**4. + + interior_salinity_c5 * zMid**5.) + else: + print('interior_salinity_shape is not supported') +# surf_indices = np.where(ds.refZMid.values >= -mixed_layer_depth_temperature)[0] +# +# if len(surf_indices) > 0: +# temperature[0,:,surf_indices] = surface_temperature + mixed_layer_temperature_gradient* \ +# ds.zMid[0,:,surf_indices].values +# +# int_indices = np.where(ds.refZMid.values < -mixed_layer_depth_temperature)[0] +# +# if len(int_indices) > 0: +# temperature[0,:,int_indices] = surface_temperature - mixed_layer_temperature_gradient* \ +# mixed_layer_depth_temperature + interior_temperature_gradient* \ +# (ds.zMid[0,:,int_indices] + \ +# mixed_layer_depth_temperature) +# +# temperature[0,:,0] += disturbance_amplitude*2*(np.random.rand(len(ds.xCell.values)) - 0.5) +# +# surf_indices = np.where(ds.refZMid.values >= -mixed_layer_depth_salinity)[0] +# if len(surf_indices) > 0: +# salinity[0,:,surf_indices] = surface_salinity - mixed_layer_salinity_gradient * \ +# ds.zMid[0,:,surf_indices] +# +# int_indices = np.where(ds.refZMid.values < -mixed_layer_depth_salinity)[0] +# if len(int_indices) > 0: +# salinity[0,:,int_indices] = surface_salinity - mixed_layer_salinity_gradient* \ +# mixed_layer_depth_salinity + interior_salinity_gradient* \ +# (ds.zMid[0,:,int_indices] + \ +# mixed_layer_depth_salinity) normalVelocity = normalVelocity.expand_dims(dim='Time', axis=0) - ds['normalVelocity'] = normalVelocity + + temperature[0, :, :] = initial_temperature + salinity[0, :, :] = initial_salinity ds['temperature'] = temperature ds['salinity'] = salinity - ds['fCell'] = coriolis_parameter * xarray.ones_like(xCell) - ds['fEdge'] = coriolis_parameter * xarray.ones_like(ds.xEdge) - ds['fVertex'] = coriolis_parameter * xarray.ones_like(ds.xVertex) + ds['fCell'] = coriolis_parameter * xr.ones_like(xCell) + ds['fEdge'] = coriolis_parameter * xr.ones_like(ds.xEdge) + ds['fVertex'] = coriolis_parameter * xr.ones_like(ds.xVertex) ds = add_boundary_arrays(ds, x_nonperiodic, y_nonperiodic) write_netcdf(ds, 'ocean.nc') + plt.figure(dpi=100) + plt.plot(initial_temperature, zMid, 'k-') + plt.xlabel('PT (C)') + plt.ylabel('z (m)') + plt.savefig(f'pt_depth_t0h.png', + bbox_inches='tight', dpi=200) + plt.close() + + plt.figure(dpi=100) + plt.plot(initial_salinity, zMid, 'k-') + plt.xlabel('SA (g/kg)') + plt.ylabel('z (m)') + plt.savefig(f'sa_depth_t0h.png', + bbox_inches='tight', dpi=200) + plt.close() diff --git a/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg b/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg index b1ad3076ad..c755c8422d 100644 --- a/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg +++ b/compass/ocean/tests/turbulence_closure/turbulence_closure.cfg @@ -40,12 +40,34 @@ mixed_layer_temperature_gradient = 0 # Salinity gradient in the mixed layer (PSU/m) mixed_layer_salinity_gradient = 0 +# Interior temperature shape, linear or polynomial +interior_temperature_shape = linear + # Temperature gradient in the interior (^oC/m) interior_temperature_gradient = 0.1 +# Temperature polynomial coefficients +interior_temperature_c0 = 271.65 +interior_temperature_c1 = -0.42 +interior_temperature_c2 = -8.0e-3 +interior_temperature_c3 = -6.75e-5 +interior_temperature_c4 = -2.66e-7 +interior_temperature_c5 = -4.02e-10 + +# Interior salinity shape, linear or polynomial +interior_salinity_shape = linear + # Salinity gradient in the interior (PSU/m) interior_salinity_gradient = 0 +# Salinity polynomial coefficients +interior_salinity_c0 = 22.36 +interior_salinity_c1 = -0.3 +interior_salinity_c2 = -3.83e-3 +interior_salinity_c3 = -2.54e-5 +interior_salinity_c4 = -8.53e-8 +interior_salinity_c5 = -1.18e-10 + # Disturbance amplitude (added to initiate convection faster; ^oC) disturbance_amplitude = 0.005 From 07eac7dc6adba7390dcc54d1232e623d3c027bd8 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 6 Mar 2023 15:31:46 -0600 Subject: [PATCH 20/21] pre-commit changes --- .../tests/turbulence_closure/__init__.py | 15 +- .../turbulence_closure/boundary_values.py | 397 +++++++++--------- .../turbulence_closure/default/__init__.py | 19 +- .../ocean/tests/turbulence_closure/forward.py | 12 +- .../tests/turbulence_closure/initial_state.py | 77 ++-- .../turbulence_closure/namelist.les.forward | 2 +- compass/ocean/tests/turbulence_closure/viz.py | 15 +- 7 files changed, 271 insertions(+), 266 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/__init__.py b/compass/ocean/tests/turbulence_closure/__init__.py index 41ef11321f..21cd27881d 100644 --- a/compass/ocean/tests/turbulence_closure/__init__.py +++ b/compass/ocean/tests/turbulence_closure/__init__.py @@ -1,7 +1,7 @@ -from compass.testgroup import TestGroup from compass.ocean.tests.turbulence_closure.decomp_test import DecompTest from compass.ocean.tests.turbulence_closure.default import Default from compass.ocean.tests.turbulence_closure.restart_test import RestartTest +from compass.testgroup import TestGroup class TurbulenceClosure(TestGroup): @@ -23,12 +23,14 @@ def __init__(self, mpas_core): for resolution in [1, 2, 1e4]: for forcing in ['cooling', 'evaporation']: self.add_test_case( - Default(test_group=self, resolution=resolution, forcing=forcing)) + Default(test_group=self, resolution=resolution, + forcing=forcing)) def configure(resolution, forcing, config): """ - Modify the configuration options for one of the turbulence closure test cases + Modify the configuration options for one of the turbulence closure test + cases Parameters ---------- @@ -61,7 +63,6 @@ def configure(resolution, forcing, config): vert_levels = 128 bottom_depth = 128.0 - config.set('turbulence_closure', 'nx', f'{nx}') config.set('turbulence_closure', 'ny', f'{ny}') config.set('turbulence_closure', 'dc', f'{resolution}') @@ -71,12 +72,14 @@ def configure(resolution, forcing, config): if forcing == 'cooling': config.set('turbulence_closure', 'surface_heat_flux', '-100') config.set('turbulence_closure', 'surface_freshwater_flux', '0') - config.set('turbulence_closure', 'interior_temperature_gradient', '0.1') + config.set('turbulence_closure', 'interior_temperature_gradient', + '0.1') config.set('turbulence_closure', 'interior_salinity_gradient', '0') config.set('turbulence_closure', 'wind_stress_zonal', '0') if forcing == 'evaporation': config.set('turbulence_closure', 'surface_heat_flux', '0') config.set('turbulence_closure', 'surface_freshwater_flux', '0.429') config.set('turbulence_closure', 'interior_temperature_gradient', '0') - config.set('turbulence_closure', 'interior_salinity_gradient', '-0.025') + config.set('turbulence_closure', 'interior_salinity_gradient', + '-0.025') config.set('turbulence_closure', 'wind_stress_zonal', '0') diff --git a/compass/ocean/tests/turbulence_closure/boundary_values.py b/compass/ocean/tests/turbulence_closure/boundary_values.py index 1da8b94589..c44fb3c0f5 100644 --- a/compass/ocean/tests/turbulence_closure/boundary_values.py +++ b/compass/ocean/tests/turbulence_closure/boundary_values.py @@ -1,193 +1,210 @@ import numpy -def add_boundary_arrays(mesh, x_nonperiodic, y_nonperiodic): - """ - Computes boundary arrays necessary for LES test cases - Values are appended to a MPAS mesh dataset - - if the domain is doubly periodic, returns array of all zeros - as arrays are meaningless in that case. - - Parameters - ---------- - mesh : xarray dataset containing an MPAS mesh - - x_nonperiodic : logical if x direction is non periodic - - y_nonperiodic : logical if y direction is non periodic - """ - - nCells = mesh.dims['nCells'] - nEdges = mesh.dims['nEdges'] - neonc = mesh.nEdgesOnCell.values - eonc = mesh.edgesOnCell.values - 1 - conc = mesh.cellsOnCell.values - 1 - nVertLevels = mesh.dims['nVertLevels'] - mlc = mesh.maxLevelCell.values - 1 - cone = mesh.cellsOnEdge.values - 1 - dv = mesh.dvEdge.values - dc = mesh.dcEdge.values - - if (not x_nonperiodic) and (not y_nonperiodic): - mesh['distanceToBoundary'] = (['nCells', 'nVertLevels'], numpy.zeros((nCells,nVertLevels))) - mesh['boundaryNormalAngle'] = (['nCells', 'nVertLevels'], numpy.zeros((nCells,nVertLevels))) - mesh['boundaryCellMask'] = (['nCells', 'nVertLevels'], numpy.zeros((nCells,nVertLevels)).astype(int)) - return mesh - #surface first and save until later - indVals = [] - angle = mesh.angleEdge.values - xc = mesh.xEdge.values - yc = mesh.yEdge.values - xCell = mesh.xCell.values - yCell = mesh.yCell.values - aCell = numpy.zeros(nCells) - bCell = numpy.zeros(nCells) - xEdge = [] - yEdge = [] - for i in range(nCells): - if conc[i,:].min() < 0: - aEdge = 0 - count = 0 - for j in range(neonc[i]): - if conc[i,j] < 0: - indVals.append(i) - xEdge.append(xc[eonc[i,j]]) - yEdge.append(yc[eonc[i,j]]) - #compute dx, dy - #check angle edge - #add pi if needed - #check if over 2pi correct if needed - dx = xc[eonc[i,j]] - xCell[i] - dy = yc[eonc[i,j]] - yCell[i] - if dx<=0 and dy<=0: # first quadrant - if angle[eonc[i,j]] >= numpy.pi/2.: - aEdge += angle[eonc[i,j]] - numpy.pi - count += 1 - else: - aEdge += angle[eonc[i,j]] - count += 1 - elif dx<=0 and dy>=0: # fourth quadrant, but reference to 0 not 2pi - if angle[eonc[i,j]] < 3.0*numpy.pi/2.0: - aEdge += angle[eonc[i,j]] + numpy.pi - count += 1 - else: - aEdge += angle[eonc[i,j]] - count += 1 - elif dx>=0 and dy>=0: #third quadrant - if angle[eonc[i,j]] < numpy.pi/2.0: # not in correct place - aEdge += angle[eonc[i,j]] + numpy.pi - count += 1 - else: - aEdge += angle[eonc[i,j]] - count += 1 - else: #quadrant 2 - if angle[eonc[i,j]] > 3.0*numpy.pi/2.0: #wrong place - aEdge += angle[eonc[i,j]] - numpy.pi - count += 1 - else: - aEdge += angle[eonc[i,j]] - count += 1 - if count > 0: - aCellTemp = aEdge / count - if aCellTemp > numpy.pi: - aCellTemp = 2.0*numpy.pi - aCellTemp - aCell[i] = aCellTemp - bCell[i] = 1 -#with angle and index arrays, now can do distance - - dist = numpy.zeros((nCells,nVertLevels)) - angleCells = numpy.zeros((nCells,nVertLevels)) - boundaryCells = numpy.zeros((nCells,nVertLevels)) - for i in range(nCells): - d = numpy.sqrt((xCell[i] - xEdge)**2 + (yCell[i] - yEdge)**2) - ind = d.argmin() - angleCells[i,0] = aCell[indVals[ind]] - boundaryCells[i,0] = bCell[i] - dist[i,0] = d.min() - - #now to the rest - for k in range(1,nVertLevels): - inds = numpy.where(k==mlc)[0] - edgeList = [] - cellList = [] - if len(inds) > 0: - #First step is forming list of edges bordering maxLC - for i in range(len(inds)): - for j in range(neonc[inds[i]]): - if conc[inds[i],j] not in inds: - edgeList.append(eonc[inds[i],j]) - for i in range(len(edgeList)): - c1 = cone[edgeList[i],0] - c2 = cone[edgeList[i],1] - if c1 in inds: - if c2 not in cellList: - cellList.append(c2) - if c2 in inds: - if c1 not in cellList: - cellList.append(c1) - - for i in range(len(cellList)): - aEdge = 0 - count = 0 - for j in range(neonc[cellList[i]]): - if eonc[cellList[i],j] in edgeList: - indVals.append(cellList[i]) - xEdge.append(xc[eonc[cellList[i],j]]) - yEdge.append(yc[eonc[cellList[i],j]]) - #compute dx, dy - #check angle edge - #add pi if needed - #check if over 2pi correct if needed - dx = xc[eonc[cellList[i],j]] - xCell[cellList[i]] - dy = yc[eonc[cellList[i],j]] - yCell[cellList[i]] - if dx<=0 and dy<=0: # first quadrant - if angle[eonc[cellList[i],j]] >= numpy.pi/2.: - aEdge += angle[eonc[cellList[i],j]] - numpy.pi - count += 1 - else: - aEdge += angle[eonc[cellList[i],j]] - count += 1 - elif dx<=0 and dy>=0: # fourth quadrant, but reference to 0 not 2pi - if angle[eonc[cellList[i],j]] < 3.0*numpy.pi/2.0: - aEdge += angle[eonc[cellList[i],j]] + numpy.pi - count += 1 - else: - aEdge += angle[eonc[cellList[i],j]] - count += 1 - elif dx>=0 and dy>=0: #third quadrant - if angle[eonc[cellList[i],j]] < numpy.pi/2.0: # not in correct place - aEdge += angle[eonc[cellList[i],j]] + numpy.pi - count += 1 - else: - aEdge += angle[eonc[cellList[i],j]] - count += 1 - else: #quadrant 2 - if angle[eonc[cellList[i],j]] > 3.0*numpy.pi/2.0: #wrong place - aEdge += angle[eonc[cellList[i],j]] - numpy.pi - count += 1 - else: - aEdge += angle[eonc[cellList[i],j]] - count += 1 - if count > 0: - aCellTemp = aEdge / count - if aCellTemp > numpy.pi: - aCellTemp = 2.0*numpy.pi - aCellTemp - aCell[cellList[i]] = aCellTemp - bCell[cellList[i]] = 1 - - for ii in range(nCells): - d = numpy.sqrt((xCell[ii] - xEdge)**2 + (yCell[ii] - yEdge)**2) - ind = d.argmin() - angleCells[ii,k] = aCell[indVals[ind]] - boundaryCells[ii,k] = bCell[ii] - dist[ii,k] = d.min() - else: - dist[:,k] = dist[:,0] - angleCells[:,k] = angleCells[:,0] - boundaryCells[:,k] = boundaryCells[:,0] - mesh['distanceToBoundary'] = (['nCells', 'nVertLevels'], dist) - mesh['boundaryNormalAngle'] = (['nCells', 'nVertLevels'], angleCells) - mesh['boundaryCellMask'] = (['nCells', 'nVertLevels'], boundaryCells.astype(int)) - - return mesh +def add_boundary_arrays(mesh, x_nonperiodic, y_nonperiodic): # noqa: C901 + """ + Computes boundary arrays necessary for LES test cases + Values are appended to a MPAS mesh dataset + + if the domain is doubly periodic, returns array of all zeros + as arrays are meaningless in that case. + + Parameters + ---------- + mesh : xarray dataset containing an MPAS mesh + + x_nonperiodic : logical if x direction is non periodic + + y_nonperiodic : logical if y direction is non periodic + """ + + nCells = mesh.dims['nCells'] + neonc = mesh.nEdgesOnCell.values + eonc = mesh.edgesOnCell.values - 1 + conc = mesh.cellsOnCell.values - 1 + nVertLevels = mesh.dims['nVertLevels'] + mlc = mesh.maxLevelCell.values - 1 + cone = mesh.cellsOnEdge.values - 1 + + if (not x_nonperiodic) and (not y_nonperiodic): + mesh['distanceToBoundary'] = ( + ['nCells', 'nVertLevels'], + numpy.zeros((nCells, nVertLevels))) + mesh['boundaryNormalAngle'] = ( + ['nCells', 'nVertLevels'], + numpy.zeros((nCells, nVertLevels))) + mesh['boundaryCellMask'] = ( + ['nCells', 'nVertLevels'], + numpy.zeros((nCells, nVertLevels)).astype(int)) + return mesh + + # surface first and save until later + indVals = [] + angle = mesh.angleEdge.values + xc = mesh.xEdge.values + yc = mesh.yEdge.values + xCell = mesh.xCell.values + yCell = mesh.yCell.values + aCell = numpy.zeros(nCells) + bCell = numpy.zeros(nCells) + xEdge = [] + yEdge = [] + for i in range(nCells): + if conc[i, :].min() < 0: + aEdge = 0 + count = 0 + for j in range(neonc[i]): + if conc[i, j] < 0: + indVals.append(i) + xEdge.append(xc[eonc[i, j]]) + yEdge.append(yc[eonc[i, j]]) + # compute dx, dy + # check angle edge + # add pi if needed + # check if over 2pi correct if needed + dx = xc[eonc[i, j]] - xCell[i] + dy = yc[eonc[i, j]] - yCell[i] + if dx <= 0 and dy <= 0: # first quadrant + if angle[eonc[i, j]] >= numpy.pi / 2.: + aEdge += angle[eonc[i, j]] - numpy.pi + count += 1 + else: + aEdge += angle[eonc[i, j]] + count += 1 + + # fourth quadrant, but reference to 0 not 2pi + elif dx <= 0 and dy >= 0: + if angle[eonc[i, j]] < 3.0 * numpy.pi / 2.0: + aEdge += angle[eonc[i, j]] + numpy.pi + count += 1 + else: + aEdge += angle[eonc[i, j]] + count += 1 + + elif dx >= 0 and dy >= 0: # third quadrant + if angle[eonc[i, j]] < numpy.pi / 2.0: + # not in correct place + aEdge += angle[eonc[i, j]] + numpy.pi + count += 1 + else: + aEdge += angle[eonc[i, j]] + count += 1 + + else: # quadrant 2 + if angle[eonc[i, j]] > 3.0 * numpy.pi / 2.0: + # wrong place + aEdge += angle[eonc[i, j]] - numpy.pi + count += 1 + else: + aEdge += angle[eonc[i, j]] + count += 1 + if count > 0: + aCellTemp = aEdge / count + if aCellTemp > numpy.pi: + aCellTemp = 2.0 * numpy.pi - aCellTemp + aCell[i] = aCellTemp + bCell[i] = 1 + + # with angle and index arrays, now can do distance + dist = numpy.zeros((nCells, nVertLevels)) + angleCells = numpy.zeros((nCells, nVertLevels)) + boundaryCells = numpy.zeros((nCells, nVertLevels)) + for i in range(nCells): + d = numpy.sqrt((xCell[i] - xEdge)**2 + (yCell[i] - yEdge)**2) + ind = d.argmin() + angleCells[i, 0] = aCell[indVals[ind]] + boundaryCells[i, 0] = bCell[i] + dist[i, 0] = d.min() + + # now to the rest + for k in range(1, nVertLevels): + inds = numpy.where(k == mlc)[0] + edgeList = [] + cellList = [] + if len(inds) > 0: + # First step is forming list of edges bordering maxLC + for i in range(len(inds)): + for j in range(neonc[inds[i]]): + if conc[inds[i], j] not in inds: + edgeList.append(eonc[inds[i], j]) + for i in range(len(edgeList)): + c1 = cone[edgeList[i], 0] + c2 = cone[edgeList[i], 1] + if c1 in inds: + if c2 not in cellList: + cellList.append(c2) + if c2 in inds: + if c1 not in cellList: + cellList.append(c1) + + for i in range(len(cellList)): + aEdge = 0 + count = 0 + for j in range(neonc[cellList[i]]): + if eonc[cellList[i], j] in edgeList: + indVals.append(cellList[i]) + xEdge.append(xc[eonc[cellList[i], j]]) + yEdge.append(yc[eonc[cellList[i], j]]) + # compute dx, dy + # check angle edge + # add pi if needed + # check if over 2pi correct if needed + dx = xc[eonc[cellList[i], j]] - xCell[cellList[i]] + dy = yc[eonc[cellList[i], j]] - yCell[cellList[i]] + if dx <= 0 and dy <= 0: # first quadrant + if angle[eonc[cellList[i], j]] >= numpy.pi / 2.: + aEdge += angle[eonc[cellList[i], j]] - numpy.pi + count += 1 + else: + aEdge += angle[eonc[cellList[i], j]] + count += 1 + # fourth quadrant, but reference to 0 not 2pi + elif dx <= 0 and dy >= 0: + if {(angle[eonc[cellList[i], j]] < + 3.0 * numpy.pi / 2.0)}: + aEdge += angle[eonc[cellList[i], j]] + numpy.pi + count += 1 + else: + aEdge += angle[eonc[cellList[i], j]] + count += 1 + elif dx >= 0 and dy >= 0: # third quadrant + # not in correct place + if angle[eonc[cellList[i], j]] < numpy.pi / 2.0: + aEdge += angle[eonc[cellList[i], j]] + numpy.pi + count += 1 + else: + aEdge += angle[eonc[cellList[i], j]] + count += 1 + else: # quadrant 2 + # wrong place + if {(angle[eonc[cellList[i], j]] > + 3.0 * numpy.pi / 2.0)}: + aEdge += angle[eonc[cellList[i], j]] - numpy.pi + count += 1 + else: + aEdge += angle[eonc[cellList[i], j]] + count += 1 + if count > 0: + aCellTemp = aEdge / count + if aCellTemp > numpy.pi: + aCellTemp = 2.0 * numpy.pi - aCellTemp + aCell[cellList[i]] = aCellTemp + bCell[cellList[i]] = 1 + + for ii in range(nCells): + d = numpy.sqrt((xCell[ii] - xEdge)**2 + + (yCell[ii] - yEdge)**2) + ind = d.argmin() + angleCells[ii, k] = aCell[indVals[ind]] + boundaryCells[ii, k] = bCell[ii] + dist[ii, k] = d.min() + else: + dist[:, k] = dist[:, 0] + angleCells[:, k] = angleCells[:, 0] + boundaryCells[:, k] = boundaryCells[:, 0] + mesh['distanceToBoundary'] = (['nCells', 'nVertLevels'], dist) + mesh['boundaryNormalAngle'] = (['nCells', 'nVertLevels'], angleCells) + mesh['boundaryCellMask'] = (['nCells', 'nVertLevels'], + boundaryCells.astype(int)) + + return mesh diff --git a/compass/ocean/tests/turbulence_closure/default/__init__.py b/compass/ocean/tests/turbulence_closure/default/__init__.py index 7a20b310ce..70483f0418 100644 --- a/compass/ocean/tests/turbulence_closure/default/__init__.py +++ b/compass/ocean/tests/turbulence_closure/default/__init__.py @@ -1,8 +1,8 @@ -from compass.testcase import TestCase -from compass.ocean.tests.turbulence_closure.initial_state import InitialState +from compass.ocean.tests import turbulence_closure from compass.ocean.tests.turbulence_closure.forward import Forward +from compass.ocean.tests.turbulence_closure.initial_state import InitialState from compass.ocean.tests.turbulence_closure.viz import Viz -from compass.ocean.tests import turbulence_closure +from compass.testcase import TestCase class Default(TestCase): @@ -45,21 +45,24 @@ def __init__(self, test_group, resolution, forcing='cooling'): if resolution <= 5: ntasks = 128 - plot_comparison=True + plot_comparison = True else: ntasks = 4 - plot_comparison=False + plot_comparison = False self.add_step( InitialState(test_case=self, resolution=resolution)) self.add_step( - Forward(test_case=self, ntasks=ntasks, openmp_threads=1, resolution=resolution)) - self.add_step(Viz(test_case=self, resolution=resolution, forcing=forcing, do_comparison=plot_comparison)) + Forward(test_case=self, ntasks=ntasks, openmp_threads=1, + resolution=resolution)) + self.add_step(Viz(test_case=self, resolution=resolution, + forcing=forcing, do_comparison=plot_comparison)) def configure(self): """ Modify the configuration options for this test case. """ - turbulence_closure.configure(self.resolution, self.forcing, self.config) + turbulence_closure.configure(self.resolution, self.forcing, + self.config) # no run() is needed because we're doing the default: running all steps diff --git a/compass/ocean/tests/turbulence_closure/forward.py b/compass/ocean/tests/turbulence_closure/forward.py index 2bbb838a5a..aac73f1af6 100644 --- a/compass/ocean/tests/turbulence_closure/forward.py +++ b/compass/ocean/tests/turbulence_closure/forward.py @@ -4,7 +4,7 @@ class Forward(Step): """ - A step for performing forward MPAS-Ocean runs as part of turbulence + A step for performing forward MPAS-Ocean runs as part of turbulence closure test cases. Attributes @@ -32,7 +32,7 @@ def __init__(self, test_case, resolution, name='forward', subdir=None, the subdirectory for the step. The default is ``name`` ntasks: int, optional - the number of tasks the step would ideally use. If fewer tasks + the number of tasks the step would ideally use. If fewer tasks are available on the system, the step will run on all available cores as long as this is not below ``min_tasks`` @@ -47,7 +47,8 @@ def __init__(self, test_case, resolution, name='forward', subdir=None, if min_tasks is None: min_tasks = ntasks super().__init__(test_case=test_case, name=name, subdir=subdir, - ntasks=ntasks, min_tasks=min_tasks, openmp_threads=openmp_threads) + ntasks=ntasks, min_tasks=min_tasks, + openmp_threads=openmp_threads) self.add_namelist_file('compass.ocean.tests.turbulence_closure', 'namelist.forward') if resolution < 500: @@ -65,8 +66,9 @@ def __init__(self, test_case, resolution, name='forward', subdir=None, self.add_input_file(filename='init.nc', target='../initial_state/ocean.nc') - self.add_input_file(filename='forcing.nc', - target='../initial_state/init_mode_forcing_data.nc') + self.add_input_file( + filename='forcing.nc', + target='../initial_state/init_mode_forcing_data.nc') self.add_input_file(filename='graph.info', target='../initial_state/culled_graph.info') diff --git a/compass/ocean/tests/turbulence_closure/initial_state.py b/compass/ocean/tests/turbulence_closure/initial_state.py index 571238e90e..d02f6da0e5 100644 --- a/compass/ocean/tests/turbulence_closure/initial_state.py +++ b/compass/ocean/tests/turbulence_closure/initial_state.py @@ -1,14 +1,13 @@ -import xarray as xr -import numpy as np import matplotlib.pyplot as plt - -from compass.ocean.tests.turbulence_closure.boundary_values import add_boundary_arrays -from random import random, seed - -from mpas_tools.planar_hex import make_planar_hex_mesh +import numpy as np +import xarray as xr from mpas_tools.io import write_netcdf from mpas_tools.mesh.conversion import convert, cull +from mpas_tools.planar_hex import make_planar_hex_mesh +from compass.ocean.tests.turbulence_closure.boundary_values import ( + add_boundary_arrays, +) from compass.ocean.vertical import init_vertical_coord from compass.step import Step @@ -61,12 +60,18 @@ def run(self): interior_temperature_shape = section.get('interior_temperature_shape') interior_salinity_shape = section.get('interior_salinity_shape') surface_salinity = section.getfloat('surface_salinity') - mixed_layer_salinity_gradient = section.getfloat('mixed_layer_salinity_gradient') - mixed_layer_temperature_gradient = section.getfloat('mixed_layer_temperature_gradient') - mixed_layer_depth_salinity = section.getfloat('mixed_layer_depth_salinity') - mixed_layer_depth_temperature = section.getfloat('mixed_layer_depth_temperature') - interior_salinity_gradient = section.getfloat('interior_salinity_gradient') - interior_temperature_gradient = section.getfloat('interior_temperature_gradient') + mixed_layer_salinity_gradient = \ + section.getfloat('mixed_layer_salinity_gradient') + mixed_layer_temperature_gradient = \ + section.getfloat('mixed_layer_temperature_gradient') + mixed_layer_depth_salinity = \ + section.getfloat('mixed_layer_depth_salinity') + mixed_layer_depth_temperature = \ + section.getfloat('mixed_layer_depth_temperature') + interior_salinity_gradient = \ + section.getfloat('interior_salinity_gradient') + interior_temperature_gradient = \ + section.getfloat('interior_temperature_gradient') interior_temperature_c0 = section.getfloat('interior_temperature_c0') interior_temperature_c1 = section.getfloat('interior_temperature_c1') interior_temperature_c2 = section.getfloat('interior_temperature_c2') @@ -86,7 +91,8 @@ def run(self): coriolis_parameter = section.getfloat('coriolis_parameter') bottom_depth = config.getfloat('vertical_grid', 'bottom_depth') - dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, nonperiodic_x=x_nonperiodic, + dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, + nonperiodic_x=x_nonperiodic, nonperiodic_y=y_nonperiodic) write_netcdf(dsMesh, 'base_mesh.nc') @@ -99,7 +105,6 @@ def run(self): ds = dsMesh.copy() dsForcing = dsMesh.copy() xCell = ds.xCell - yCell = ds.yCell ds['bottomDepth'] = bottom_depth * xr.ones_like(xCell) ds['ssh'] = xr.zeros_like(xCell) @@ -107,7 +112,8 @@ def run(self): init_vertical_coord(config, ds) dsForcing['windStressZonal'] = wind_stress_zonal * xr.ones_like(xCell) - dsForcing['windStressMeridional'] = wind_stress_meridional * xr.ones_like(xCell) + dsForcing['windStressMeridional'] = \ + wind_stress_meridional * xr.ones_like(xCell) dsForcing['sensibleHeatFlux'] = surface_heat_flux * xr.ones_like(xCell) dsForcing['rainFlux'] = surface_freshwater_flux * xr.ones_like(xCell) @@ -124,7 +130,8 @@ def run(self): zMid = ds.refZMid.values temperature_at_mld = (surface_temperature - - mixed_layer_temperature_gradient * mixed_layer_depth_temperature) + mixed_layer_temperature_gradient * + mixed_layer_depth_temperature) if interior_temperature_shape == 'linear': initial_temperature = np.where( zMid >= -mixed_layer_depth_temperature, @@ -144,7 +151,8 @@ def run(self): print('interior_temperature_shape is not supported') salinity_at_mld = (surface_salinity - - mixed_layer_salinity_gradient * mixed_layer_depth_salinity) + mixed_layer_salinity_gradient * + mixed_layer_depth_salinity) if interior_salinity_shape == 'linear': initial_salinity = np.where( zMid >= -mixed_layer_depth_salinity, @@ -162,38 +170,13 @@ def run(self): interior_salinity_c5 * zMid**5.) else: print('interior_salinity_shape is not supported') -# surf_indices = np.where(ds.refZMid.values >= -mixed_layer_depth_temperature)[0] -# -# if len(surf_indices) > 0: -# temperature[0,:,surf_indices] = surface_temperature + mixed_layer_temperature_gradient* \ -# ds.zMid[0,:,surf_indices].values -# -# int_indices = np.where(ds.refZMid.values < -mixed_layer_depth_temperature)[0] -# -# if len(int_indices) > 0: -# temperature[0,:,int_indices] = surface_temperature - mixed_layer_temperature_gradient* \ -# mixed_layer_depth_temperature + interior_temperature_gradient* \ -# (ds.zMid[0,:,int_indices] + \ -# mixed_layer_depth_temperature) -# -# temperature[0,:,0] += disturbance_amplitude*2*(np.random.rand(len(ds.xCell.values)) - 0.5) -# -# surf_indices = np.where(ds.refZMid.values >= -mixed_layer_depth_salinity)[0] -# if len(surf_indices) > 0: -# salinity[0,:,surf_indices] = surface_salinity - mixed_layer_salinity_gradient * \ -# ds.zMid[0,:,surf_indices] -# -# int_indices = np.where(ds.refZMid.values < -mixed_layer_depth_salinity)[0] -# if len(int_indices) > 0: -# salinity[0,:,int_indices] = surface_salinity - mixed_layer_salinity_gradient* \ -# mixed_layer_depth_salinity + interior_salinity_gradient* \ -# (ds.zMid[0,:,int_indices] + \ -# mixed_layer_depth_salinity) normalVelocity = normalVelocity.expand_dims(dim='Time', axis=0) ds['normalVelocity'] = normalVelocity temperature[0, :, :] = initial_temperature + temperature[0, :, 0] += disturbance_amplitude * 2. * \ + (np.random.rand(len(ds.xCell.values)) - 0.5) salinity[0, :, :] = initial_salinity ds['temperature'] = temperature ds['salinity'] = salinity @@ -209,7 +192,7 @@ def run(self): plt.plot(initial_temperature, zMid, 'k-') plt.xlabel('PT (C)') plt.ylabel('z (m)') - plt.savefig(f'pt_depth_t0h.png', + plt.savefig('pt_depth_t0h.png', bbox_inches='tight', dpi=200) plt.close() @@ -217,6 +200,6 @@ def run(self): plt.plot(initial_salinity, zMid, 'k-') plt.xlabel('SA (g/kg)') plt.ylabel('z (m)') - plt.savefig(f'sa_depth_t0h.png', + plt.savefig('sa_depth_t0h.png', bbox_inches='tight', dpi=200) plt.close() diff --git a/compass/ocean/tests/turbulence_closure/namelist.les.forward b/compass/ocean/tests/turbulence_closure/namelist.les.forward index 00ea3f1591..f26b031187 100644 --- a/compass/ocean/tests/turbulence_closure/namelist.les.forward +++ b/compass/ocean/tests/turbulence_closure/namelist.les.forward @@ -9,5 +9,5 @@ config_tke_ini = 1.0e-30 config_LES_noise_enable = .true. config_LES_noise_max_time = 150 config_LES_noise_min_level = 5 -config_LES_noise_max_level = 45 +config_LES_noise_max_level = 45 config_LES_noise_perturbation_magnitude = 1e-4 diff --git a/compass/ocean/tests/turbulence_closure/viz.py b/compass/ocean/tests/turbulence_closure/viz.py index 2f34651c5a..405182ad9b 100644 --- a/compass/ocean/tests/turbulence_closure/viz.py +++ b/compass/ocean/tests/turbulence_closure/viz.py @@ -1,6 +1,6 @@ -import xarray as xr -import numpy as np import matplotlib.pyplot as plt +import numpy as np +import xarray as xr from compass.step import Step @@ -25,7 +25,6 @@ def __init__(self, test_case, resolution, forcing, do_comparison=False): self.add_input_file(filename='mesh.nc', target='../initial_state/culled_mesh.nc') - if do_comparison: if forcing == 'cooling': forcing_name = 'c02' @@ -68,9 +67,9 @@ def run(self): else: # This routine is not generic but should not be used as # daysSinceStartOfSim is included in the streams file - time0 = 2.0 + (7.0/24.0) - dt = 0.5/24.0 - time = np.linspace(time0 + dt, time0 + nt*dt, num=nt) + time0 = 2.0 + (7.0 / 24.0) + dt = 0.5 / 24.0 + time = np.linspace(time0 + dt, time0 + nt * dt, num=nt) if self.do_comparison: time_palm = ds_palm.time @@ -87,8 +86,6 @@ def run(self): ds = ds.sortby('yEdge') # Get mesh variables - xEdge = dsInit.xEdge - yEdge = dsInit.yEdge xCell = dsInit.xCell yCell = dsInit.yCell zMid = dsInit.refZMid @@ -188,7 +185,7 @@ def run(self): plt.scatter(np.divide(xCell, 1e3), np.divide(yCell, 1e3), s=markersize, c=w_zTop.values, - cmap='cmo.balance', vmin=-1.*cmax, vmax=cmax) + cmap='cmo.balance', vmin=-1. * cmax, vmax=cmax) plt.xlabel('x (km)') plt.ylabel('y (km)') cbar = plt.colorbar() From 7d5bcd388b4b72da2b99d25003312553cd27869f Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 17 Mar 2023 17:30:18 -0500 Subject: [PATCH 21/21] Update nonhydro config options --- .../turbulence_closure/namelist.nonhydro.forward | 3 --- .../ocean/test_groups/turbulence_closure.rst | 15 +++++++++++---- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/compass/ocean/tests/turbulence_closure/namelist.nonhydro.forward b/compass/ocean/tests/turbulence_closure/namelist.nonhydro.forward index 281f2de6ce..c5b9b0c098 100644 --- a/compass/ocean/tests/turbulence_closure/namelist.nonhydro.forward +++ b/compass/ocean/tests/turbulence_closure/namelist.nonhydro.forward @@ -1,8 +1,5 @@ config_enable_nonhydrostatic_mode = .true. -config_triskCV = .false. config_nonhydrostatic_preconditioner = 'bjacobi' config_nonhydrostatic_solver_type = 'cg' config_nonhydrostatic_solve_surface_boundary_condition = 'pressureTopGradientBottom' -config_use_constant_forced_pgf = .false. config_thickness_flux_type = 'upwind' -config_aust_scale_factor = 0.70710678 diff --git a/docs/users_guide/ocean/test_groups/turbulence_closure.rst b/docs/users_guide/ocean/test_groups/turbulence_closure.rst index 4afec83860..a4979a5269 100644 --- a/docs/users_guide/ocean/test_groups/turbulence_closure.rst +++ b/docs/users_guide/ocean/test_groups/turbulence_closure.rst @@ -70,14 +70,21 @@ active with the following config options: .. code-block:: cfg - config_enable_nonhydrostatic_mode = .true. - config_triskCV = .false. + # turns the nonhydro model on or off + config_enable_nonhydrostatic_mode = .true. + + # preconditioner for the linear solver. Other options can be used, like + # jacobi, sor and asm, but they were found to be slower than block jacobi. config_nonhydrostatic_preconditioner = 'bjacobi' + + # linear solver. Other options can be used, like gmres config_nonhydrostatic_solver_type = 'cg' + + # do not change! config_nonhydrostatic_solve_surface_boundary_condition = 'pressureTopGradientBottom' - config_use_constant_forced_pgf = .false. + + # Nonhydro will work with either 'centered' or 'upwind'. config_thickness_flux_type = 'upwind' - config_aust_scale_factor = 0.70710678 default -------