Source code for compass.ocean.tests.gotm.default.analysis

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

from compass.step import Step


[docs] class Analysis(Step): """ A step for plotting the results of the default General Ocean Turbulence Model (GOTM) test case """
[docs] def __init__(self, test_case): """ Create a new test case Parameters ---------- test_case : compass.ocean.tests.gotm.default.Default The test case this step belongs to """ super().__init__(test_case=test_case, name='analysis', ntasks=1, min_tasks=1, openmp_threads=1) self.add_input_file(filename='output.nc', target='../forward/output.nc') self.add_output_file(filename='velocity_profile.png') self.add_output_file(filename='viscosity_profile.png')
[docs] def run(self): """ Run this step of the test case """ # render statically by default plt.switch_backend('agg') # constants kappa = 0.4 z0b = 1.5e-3 gssh = 1e-5 g = 9.81 h = 15 # load output ds = xr.open_dataset('output.nc') # velocity u = ds.velocityZonal.isel(Time=-1, nCells=0).values # viscosity nu = ds.vertViscTopOfCell.isel(Time=-1, nCells=0).values # depth bottom_depth = ds.refBottomDepth.values z = np.zeros_like(bottom_depth) z[0] = -0.5*bottom_depth[0] z[1:] = -0.5*(bottom_depth[0:-1]+bottom_depth[1:]) zi = np.zeros(bottom_depth.size+1) zi[0] = 0.0 zi[1:] = -bottom_depth[0:] # analytical solution ustarb = np.sqrt(g*h*gssh) u_a = ustarb/kappa*np.log((z0b+z+h)/z0b) nu_a = -ustarb/h*kappa*(z0b+z+h)*z # infered drag coefficient cd = ustarb**2/u_a[-1]**2 self.logger.info('C_d = {:6.4g}'.format(cd)) # plot velocity plt.figure() plt.plot(u_a, z, 'k--', label='Analytical') plt.plot(u, z, 'k-', label='GOTM') plt.xlabel('Velocity (m/s)') plt.ylabel('Depth (m)') plt.legend() plt.savefig('velocity_profile.png') # plot viscosity plt.figure() plt.plot(nu_a, z, 'k--', label='Analytical') plt.plot(nu, zi, 'k-', label='GOTM') plt.xlabel('Viscosity (m$^2$/s)') plt.ylabel('Depth (m)') plt.legend() plt.savefig('viscosity_profile.png')