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()