r"""
Generic Length Scale Turbulence Closure model
=============================================
This model solves two dynamic equations, for turbulent kinetic energy
(TKE, :math:`k`) and one for an additional variable, the generic length scale
:math:`\psi` [1]:
.. math::
\frac{\partial k}{\partial t} + \nabla_h \cdot (\textbf{u} k)
+ \frac{\partial (w k)}{\partial z}
= \frac{\partial}{\partial z}\left(\frac{\nu}{\sigma_k} \frac{\partial k}{\partial z}\right)
+ P + B - \varepsilon
:label: turb_tke_eq
.. math::
\frac{\partial \psi}{\partial t} + \nabla_h \cdot (\textbf{u} \psi)
+ \frac{\partial (w \psi)}{\partial z}
= \frac{\partial}{\partial z}\left(\frac{\nu}{\sigma_\psi} \frac{\partial \psi}{\partial z}\right)
+ \frac{\psi}{k} (c_1 P + c_3 B - c_2 \varepsilon f_{wall})
:label: turb_psi_eq
with the production :math:`P` and buoyancy production :math:`B` are
.. math::
P &= \nu M^2 \\
B &= -\mu N^2
where :math:`M` and :math:`N` are the shear and buoyancy frequencies
.. math::
M^2 &= \left(\frac{\partial u}{\partial z}\right)^2
+ \left(\frac{\partial v}{\partial z}\right)^2 \\
N^2 &= -\frac{g}{\rho_0}\frac{\partial \rho}{\partial z}
The generic length scale variable is defined as
.. math::
\psi = (c_\mu^0)^p k^m l^n
where :math:`p, m, n` are parameters and :math:`c_\mu^0` is an empirical constant.
The parameters :math:`c_1,c_2,c_3,f_{wall}` depend on the chosen closure.
The parameter :math:`c_3` takes two values: :math:`c_3^-` in stably stratified
regime, and :math:`c_3^+` in unstably stratified cases.
Turbulent length scale :math:`l`, and the TKE dissipation rate
:math:`\varepsilon` are obtained diagnostically as
.. math::
l &= (c_\mu^0)^3 k^{3/2} \varepsilon^{-1} \\
\varepsilon &= (c_\mu^0)^{3+p/n} k^{3/2 + m/n} \psi^{-1/n}
Finally the vertical eddy viscosity and diffusivity are also computed
diagnostically
.. math::
\nu &= \sqrt{2k} l S_m \\
\mu &= \sqrt{2k} l S_\rho
Stability functions :math:`S_m` and :math:`S_\rho` are defined in [2] or [3].
Implementation follows [4].
Supported closures
------------------
The GLS model parameters are controlled via the :class:`.GLSModelOptions` class.
The parameters can be accessed from the solver object:
.. code-block:: python
solver = FlowSolver(...)
solver.options.turbulence_model_type = 'gls' # activate GLS model (default)
turbulence_model_options = solver.options.turbulence_model_options
turbulence_model_options.closure_name = 'k-omega'
turbulence_model_options.stability_function_name = 'CB'
turbulence_model_options.compute_c3_minus = True
Currently the following closures have been implemented:
- :math:`k-\varepsilon` model
This closure is obtained with :math:`p=3, m=3/2, n=-1`, resulting in
:math:`\psi=\varepsilon`.
To use this option set ``closure_name = k-epsilon``
- :math:`k-\omega` model
This closure is obtained with :math:`p=-1, m=1/2, n=-1`, resulting in
:math:`\psi=\omega`.
To use this option set ``closure_name = k-omega``
- GLS model A
This closure is obtained with :math:`p=2, m=1, n=-2/3`, resulting in
:math:`\psi=\omega`.
To use this option set ``closure_name = gen``
The following stability functions have been implemented
- Canuto A [3]
To use this option set ``closure_name = CA``
- Canuto B [3]
To use this option set ``closure_name = CB``
- Kantha-Clayson [2]
To use this option set ``closure_name = KC``
- Cheng [6]
To use this option set ``closure_name = CH``
See :mod:`.stability_functions` for more information.
[1] Umlauf, L. and Burchard, H. (2003). A generic length-scale equation for
geophysical turbulence models. Journal of Marine Research, 61:235--265(31).
http://dx.doi.org/10.1357/002224003322005087
[2] Kantha, L. H. and Clayson, C. A. (1994). An improved mixed layer model for
geophysical applications. Journal of Geophysical Research: Oceans,
99(C12):25235--25266.
http://dx.doi.org/10.1029/94JC02257
[3] Canuto et al. (2001). Ocean Turbulence. Part I: One-Point Closure Model -
Momentum and Heat Vertical Diffusivities. Journal of Physical Oceanography,
31(6):1413-1426.
http://dx.doi.org/10.1175/1520-0485(2001)031
[4] Warner et al. (2005). Performance of four turbulence closure models
implemented using a generic length scale method. Ocean Modelling,
8(1-2):81--113.
http://dx.doi.org/10.1016/j.ocemod.2003.12.003
[5] Umlauf, L. and Burchard, H. (2005). Second-order turbulence closure models
for geophysical boundary layers. A review of recent work. Continental Shelf
Research, 25(7-8):795--827.
http://dx.doi.org/10.1016/j.csr.2004.08.004
[6] Cheng et al. (2002). An improved model for the turbulent PBL.
J. Atmos. Sci., 59:1550-1565.
http://dx.doi.org/10.1175/1520-0469(2002)059%3C1550:AIMFTT%3E2.0.CO;2
[7] Burchard and Petersen (1999). Models of turbulence in the marine
environment - a comparative study of two-equation turbulence models.
Journal of Marine Systems, 21(1-4):29-53.
http://dx.doi.org/10.1016/S0924-7963(99)00004-4
"""
from .utility import *
from .equation import Equation
from .tracer_eq import *
from .stability_functions import *
from .log import *
from .options import GLSModelOptions, PacanowskiPhilanderModelOptions
import numpy
from abc import ABC, abstractmethod
from pyop2.profiling import timed_stage
[docs]
def set_func_min_val(f, minval):
"""
Sets a minimum value to a :class:`Function`
"""
f.dat.data[f.dat.data < minval] = minval
[docs]
def set_func_max_val(f, maxval):
"""
Sets a minimum value to a :class:`Function`
"""
f.dat.data[f.dat.data > maxval] = maxval
[docs]
class VerticalGradSolver(object):
"""
Computes vertical gradient in the weak sense.
"""
@PETSc.Log.EventDecorator("thetis.VerticalGradSolver.__init__")
def __init__(self, source, solution, solver_parameters=None):
"""
:arg source: A :class:`Function` or expression to differentiate.
:arg solution: A :class:`Function` where the solution will be stored.
Must be in P0 space.
:kwarg dict solver_parameters: PETSc solver options
"""
if solver_parameters is None:
solver_parameters = {}
solver_parameters.setdefault('snes_type', 'ksponly')
solver_parameters.setdefault('ksp_type', 'preonly')
solver_parameters.setdefault('pc_type', 'bjacobi')
solver_parameters.setdefault('sub_ksp_type', 'preonly')
solver_parameters.setdefault('sub_pc_type', 'ilu')
self.source = source
self.solution = solution
self.fs = self.solution.function_space()
self.mesh = self.fs.mesh()
# weak gradient evaluator
test = TestFunction(self.fs)
tri = TrialFunction(self.fs)
normal = FacetNormal(self.mesh)
a = inner(test, tri)*dx
p = self.source
l = -inner(p, Dx(test, 2))*dx
l += avg(p)*jump(test, normal[2])*dS_h
l += p*test*normal[2]*(ds_t + ds_b)
prob = LinearVariationalProblem(a, l, self.solution, constant_jacobian=True)
self.weak_grad_solver = LinearVariationalSolver(prob, solver_parameters=solver_parameters)
[docs]
@PETSc.Log.EventDecorator("thetis.VerticalGradSolver.solve")
def solve(self):
"""Computes the gradient"""
self.weak_grad_solver.solve()
[docs]
class ShearFrequencySolver(object):
r"""
Computes vertical shear frequency squared form the given horizontal
velocity field.
.. math::
M^2 = \left(\frac{\partial u}{\partial z}\right)^2
+ \left(\frac{\partial v}{\partial z}\right)^2
"""
@PETSc.Log.EventDecorator("thetis.ShearFrequencySolver.__init__")
def __init__(self, uv, m2, mu, mv, mu_tmp, relaxation=1.0, minval=1e-12):
"""
:arg uv: horizontal velocity field
:type uv: :class:`Function`
:arg m2: :math:`M^2` field
:type m2: :class:`Function`
:arg mu: field for x component of :math:`M^2`
:type mu: :class:`Function`
:arg mv: field for y component of :math:`M^2`
:type mv: :class:`Function`
:arg mu_tmp: temporary field
:type mu_tmp: :class:`Function`
:kwarg float relaxation: relaxation coefficient for mixing old and new values
M2 = relaxation*M2_new + (1-relaxation)*M2_old
:kwarg float minval: minimum value for :math:`M^2`
"""
self.mu = mu
self.mv = mv
self.m2 = m2
self.mu_tmp = mu_tmp
self.minval = minval
self.relaxation = relaxation
self.var_solvers = []
for i_comp in range(2):
self.var_solvers.append(VerticalGradSolver(uv[i_comp], mu_tmp))
[docs]
@PETSc.Log.EventDecorator("thetis.ShearFrequencySolver.solve")
def solve(self, init_solve=False):
"""
Computes buoyancy frequency
:kwarg bool init_solve: Set to True if solving for the first time, skips
relaxation
"""
with timed_stage('shear_freq_solv'):
mu_comp = [self.mu, self.mv]
self.m2.assign(0.0)
for i_comp, solver in enumerate(self.var_solvers):
solver.solve()
gamma = self.relaxation if not init_solve else 1.0
mu_comp[i_comp].assign(gamma*self.mu_tmp
+ (1.0 - gamma)*mu_comp[i_comp])
self.m2.interpolate(self.m2 + mu_comp[i_comp]*mu_comp[i_comp])
# crop small/negative values
set_func_min_val(self.m2, self.minval)
[docs]
class BuoyFrequencySolver(object):
r"""
Computes buoyancy frequency squared form the given horizontal
velocity field.
.. math::
N^2 = -\frac{g}{\rho_0}\frac{\partial \rho}{\partial z}
"""
@PETSc.Log.EventDecorator("thetis.BuoyFrequencySolver.__init__")
def __init__(self, rho, n2, n2_tmp, relaxation=1.0, minval=1e-12):
"""
:arg rho: water density field
:type rho: :class:`Function`
:arg n2: :math:`N^2` field
:type n2: :class:`Function`
:arg n2_tmp: temporary field
:type n2_tmp: :class:`Function`
:kwarg float relaxation: relaxation coefficient for mixing old and new
values N2 = relaxation*N2_new + (1-relaxation)*N2_old
:kwarg float minval: minimum value for :math:`N^2`
"""
self._no_op = False
if rho is None:
self._no_op = True
if not self._no_op:
self.n2 = n2
self.n2_tmp = n2_tmp
self.relaxation = relaxation
g = physical_constants['g_grav']
rho0 = physical_constants['rho0']
p = -g/rho0 * rho
solver = VerticalGradSolver(p, self.n2_tmp)
self.var_solver = solver
[docs]
@PETSc.Log.EventDecorator("thetis.BuoyFrequencySolver.solve")
def solve(self, init_solve=False):
"""
Computes buoyancy frequency
:kwarg bool init_solve: Set to True if solving for the first time, skips
relaxation
"""
with timed_stage('buoy_freq_solv'):
if not self._no_op:
self.var_solver.solve()
gamma = self.relaxation if not init_solve else 1.0
self.n2.assign(gamma*self.n2_tmp
+ (1.0 - gamma)*self.n2)
[docs]
class TurbulenceModel(ABC):
"""Base class for all vertical turbulence models"""
[docs]
@abstractmethod
def initialize(self):
"""Initialize all turbulence fields"""
pass
[docs]
@abstractmethod
def preprocess(self, init_solve=False):
"""
Computes all diagnostic variables that depend on the mean flow model
variables.
To be called before updating the turbulence PDEs.
"""
pass
[docs]
@abstractmethod
def postprocess(self):
"""
Updates all diagnostic variables that depend on the turbulence state
variables.
To be called after updating the turbulence PDEs.
"""
pass
[docs]
class GenericLengthScaleModel(TurbulenceModel):
"""
Generic Length Scale turbulence closure model implementation
"""
@PETSc.Log.EventDecorator("thetis.GenericLengthScaleModel.__init__")
def __init__(self, solver, k_field, psi_field, uv_field, rho_field,
l_field, epsilon_field,
eddy_diffusivity, eddy_viscosity,
n2, m2, options=None):
"""
:arg solver: FlowSolver object
:arg k_field: turbulent kinetic energy (TKE) field
:type k_field: :class:`Function`
:arg psi_field: generic length scale field
:type psi_field: :class:`Function`
:arg uv_field: horizontal velocity field
:type uv_field: :class:`Function`
:arg rho_field: water density field
:type rho_field: :class:`Function`
:arg l_field: turbulence length scale field
:type l_field: :class:`Function`
:arg epsilon_field: TKE dissipation rate field
:type epsilon_field: :class:`Function`
:arg eddy_diffusivity: eddy diffusivity field
:type eddy_diffusivity: :class:`Function`
:arg eddy_viscosity: eddy viscosity field
:type eddy_viscosity: :class:`Function`
:arg n2: field for buoyancy frequency squared
:type n2: :class:`Function`
:arg m2: field for vertical shear frequency squared
:type m2: :class:`Function`
:kwarg options: GLS model options
"""
self.solver = solver
# 3d model fields
self.uv = uv_field
self.rho = rho_field
# prognostic fields
self.k = k_field
self.psi = psi_field
# diagnostic fields
# NOTE redundant for k-epsilon model where psi==epsilon
self.epsilon = epsilon_field
self.l = l_field
self.viscosity = eddy_viscosity
self.diffusivity = eddy_diffusivity
self.n2 = n2
self.m2 = m2
self.mu_tmp = Function(self.m2.function_space(),
name='tmp Shear frequency')
self.mu = Function(self.m2.function_space(), name='Shear frequency X')
self.mv = Function(self.m2.function_space(), name='Shear frequency Y')
self.n2_tmp = Function(self.n2.function_space(),
name='tmp buoyancy frequency')
self.n2_pos = Function(self.n2.function_space(),
name='positive buoyancy frequency')
self.n2_neg = Function(self.n2.function_space(),
name='negative buoyancy frequency')
# parameter to mix old and new viscosity values (1 => new only)
self.relaxation = 1.0
if options is not None:
self.options = options
else:
self.options = GLSModelOptions()
o = self.options
stability_function_name = o.stability_function_name
stab_args = {'lim_alpha_shear': True,
'lim_alpha_buoy': True,
'smooth_alpha_buoy_lim': False}
if stability_function_name == 'Kantha-Clayson':
self.stability_func = StabilityFunctionKanthaClayson(**stab_args)
elif stability_function_name == 'Canuto A':
self.stability_func = StabilityFunctionCanutoA(**stab_args)
elif stability_function_name == 'Canuto B':
self.stability_func = StabilityFunctionCanutoB(**stab_args)
elif stability_function_name == 'Cheng':
self.stability_func = StabilityFunctionCheng(**stab_args)
else:
raise Exception('Unknown stability function type: ' + stability_function_name)
if o.compute_cmu0:
o.cmu0 = self.stability_func.compute_cmu0()
if o.compute_kappa:
kappa = self.stability_func.compute_kappa(o.schmidt_nb_psi, o.cmu0, o.n, o.c1, o.c2)
o.kappa = kappa
# update mean flow model value as well
physical_constants['von_karman'].assign(kappa)
elif o.compute_schmidt_nb_psi:
o.schmidt_nb_psi = self.stability_func.compute_sigma_psi(o.kappa, o.cmu0, o.n, o.c1, o.c2)
if o.compute_c3_minus:
o.c3_minus = self.stability_func.compute_c3_minus(o.c1, o.c2, o.ri_st)
if o.compute_psi_min:
o.psi_min = (o.cmu0**(3.0 + o.p/o.n)*o.k_min**(3.0/2.0 + o.m/o.n)*o.eps_min**(-1.0))**o.n
# minimum length scale
if o.compute_len_min:
o.len_min = o.cmu0**3 * o.k_min**1.5 / o.eps_min
if o.compute_galperin_clim:
o.galperin_clim = self.stability_func.compute_length_clim(o.cmu0, o.ri_st)
self.shear_frequency_solver = ShearFrequencySolver(self.uv, self.m2,
self.mu, self.mv,
self.mu_tmp)
if self.rho is not None:
self.buoy_frequency_solver = BuoyFrequencySolver(self.rho, self.n2,
self.n2_tmp)
print_output(self.options)
self._initialized = False
[docs]
@PETSc.Log.EventDecorator("thetis.GenericLengthScaleModel.initialize")
def initialize(self, l_init=0.01):
"""
Initialize turbulent fields.
k is set to k_min value. epsilon and psi are set to a value that
corresponds to l_init length scale.
:arg l_init: initial value for turbulent length scale.
"""
if not self.solver._simulation_continued:
o = self.options
eps_init = o.cmu0**3.0 * o.k_min**(3.0/2.0) / l_init
psi_init = (o.cmu0**(3.0 + o.p/o.n)*o.k_min**(3.0/2.0 + o.m/o.n)*eps_init**(-1.0))**o.n
self.k.assign(self.options.k_min)
self.psi.assign(psi_init)
self.epsilon.assign(eps_init)
self.n2.assign(1e-12)
self.n2_pos.assign(1e-12)
self.n2_neg.assign(0.0)
self._initialized = True
self.preprocess(init_solve=True)
self.postprocess()
[docs]
@PETSc.Log.EventDecorator("thetis.GenericLengthScaleModel.preprocess")
def preprocess(self, init_solve=False):
"""
Computes all diagnostic variables that depend on the mean flow model
variables.
To be called before updating the turbulence PDEs.
"""
if self._initialized is False:
self.initialize()
self.shear_frequency_solver.solve(init_solve=init_solve)
if self.rho is not None:
self.buoy_frequency_solver.solve(init_solve=init_solve)
# split to positive and negative parts
self.n2_pos.assign(1e-12)
self.n2_neg.assign(0.0)
pos_ix = self.n2.dat.data[:] >= 0.0
self.n2_pos.dat.data[pos_ix] = self.n2.dat.data[pos_ix]
self.n2_neg.dat.data[~pos_ix] = self.n2.dat.data[~pos_ix]
[docs]
@PETSc.Log.EventDecorator("thetis.GenericLengthScaleModel.postprocess")
def postprocess(self):
r"""
Updates all diagnostic variables that depend on the turbulence state
variables :math:`k,\psi`.
To be called after updating the turbulence PDEs.
"""
with timed_stage('turb_postproc'):
o = self.options
cmu0 = o.cmu0
p = o.p
n = o.n
m = o.m
# limit k
set_func_min_val(self.k, o.k_min)
k_arr = self.k.dat.data[:]
n2_pos = self.n2_pos.dat.data[:]
n2_pos_eps = 1e-12
galp_clim = o.galperin_clim
if o.limit_psi:
# impose Galperin limit on psi
# psi^(1/n) <= sqrt(0.56)* (cmu0)^(p/n) *k^(m/n+0.5)* n2^(-0.5)
val = (numpy.sqrt(2)*galp_clim * (cmu0)**(p / n) * k_arr**(m / n + 0.5) * (n2_pos + n2_pos_eps)**(-0.5))**n
if n > 0:
# impose max value
numpy.minimum(self.psi.dat.data, val, self.psi.dat.data)
else:
# impose min value
numpy.maximum(self.psi.dat.data, val, self.psi.dat.data)
set_func_min_val(self.psi, o.psi_min)
# udpate epsilon
self.epsilon.interpolate(cmu0**(3.0 + p/n)*self.k**(3.0/2.0 + m/n)*self.psi**(-1.0/n))
if o.limit_eps:
# impose Galperin limit on eps
eps_min = cmu0**3.0/(numpy.sqrt(2)*galp_clim)*numpy.sqrt(n2_pos)*k_arr
numpy.maximum(self.epsilon.dat.data, eps_min, self.epsilon.dat.data)
# impose minimum value
set_func_min_val(self.epsilon, o.eps_min)
# update L
self.l.interpolate(cmu0**3.0 * self.k**(3.0/2.0) / self.epsilon)
if o.limit_len_min:
set_func_min_val(self.l, o.len_min)
if o.limit_len:
# Galperin length scale limitation
len_max = galp_clim*numpy.sqrt(2*k_arr/(n2_pos + n2_pos_eps))
numpy.minimum(self.l.dat.data, len_max, self.l.dat.data)
if self.l.dat.data.max() > 10.0:
warning(' * large L: {:f}'.format(self.l.dat.data.max()))
# update stability functions
s_m, s_h = self.stability_func.evaluate(self.m2.dat.data,
self.n2.dat.data,
self.k.dat.data,
self.epsilon.dat.data)
# update diffusivity/viscosity
b = numpy.sqrt(self.k.dat.data[:])*self.l.dat.data[:]
lam = self.relaxation
new_visc = b*s_m/cmu0**3
new_diff = b*s_h/cmu0**3
self.viscosity.dat.data[:] = lam*new_visc + (1.0 - lam)*self.viscosity.dat.data[:]
self.diffusivity.dat.data[:] = lam*new_diff + (1.0 - lam)*self.diffusivity.dat.data[:]
set_func_min_val(self.viscosity, o.visc_min)
set_func_min_val(self.diffusivity, o.diff_min)
[docs]
def print_debug(self):
"""
Print diagnostic field values for debugging.
"""
fmt = '{:8s} {:10.3e} {:10.3e}'
print_output(' -----')
print_output(fmt.format('k', self.k.dat.data.min(), self.k.dat.data.max()))
print_output(fmt.format('psi', self.psi.dat.data.min(), self.psi.dat.data.max()))
print_output(fmt.format('eps', self.epsilon.dat.data.min(), self.epsilon.dat.data.max()))
print_output(fmt.format('L', self.l.dat.data.min(), self.l.dat.data.max()))
print_output(fmt.format('M2', self.m2.dat.data.min(), self.m2.dat.data.max()))
print_output(fmt.format('N2', self.n2.dat.data.min(), self.n2.dat.data.max()))
print_output(fmt.format('N2+', self.n2_pos.dat.data.min(), self.n2_pos.dat.data.max()))
print_output(fmt.format('N2-', self.n2_neg.dat.data.min(), self.n2_neg.dat.data.max()))
print_output(fmt.format('s_h', s_h.min(), s_h.max()))
print_output(fmt.format('s_m', s_m.min(), s_m.max()))
print_output(fmt.format('nuv', self.viscosity.dat.data.min(), self.viscosity.dat.data.max()))
print_output(fmt.format('muv', self.diffusivity.dat.data.min(), self.diffusivity.dat.data.max()))
[docs]
class TKESourceTerm(TracerTerm):
r"""
Production and destruction terms of the TKE equation :eq:`turb_tke_eq`
.. math::
F_k = P + B - \varepsilon
To ensure positivity we use Patankar-type time discretization: all source
terms are treated explicitly and sink terms are treated implicitly.
To this end the buoyancy production term :math:`B` is split in two:
.. math::
F_k = P + B_{source} + \frac{k^{n+1}}{k^n}(B_{sink} - \varepsilon)
with :math:`B_{source} \ge 0` and :math:`B_{sink} < 0`.
"""
def __init__(self, function_space, gls_model,
bathymetry=None, v_elem_size=None, h_elem_size=None):
"""
:arg function_space: :class:`FunctionSpace` where the solution belongs
:arg gls_model: :class:`.GenericLengthScaleModel` object
:kwarg bathymetry: bathymetry of the domain
:type bathymetry: 3D :class:`Function` or :class:`Constant`
:kwarg v_elem_size: scalar :class:`Function` that defines the vertical
element size
:kwarg h_elem_size: scalar :class:`Function` that defines the horizontal
element size
"""
super(TKESourceTerm, self).__init__(function_space,
bathymetry, v_elem_size, h_elem_size)
self.gls_model = gls_model
[docs]
def residual(self, solution, solution_old, fields, fields_old, bnd_conditions):
# TKE: P + B - eps
# P = viscosity M**2 (production)
# B = - diffusivity N**2 (byoyancy production)
# M**2 = (du/dz)**2 + (dv/dz)**2 (shear frequency)
# N**2 = -g\rho_0 (drho/dz) (buoyancy frequency)
# eps = (cmu0)**(3+p/n)*tke**(3/2+m/n)*psi**(-1/n)
# (tke dissipation rate)
eddy_viscosity = fields_old['viscosity_v']
eddy_diffusivity = fields_old['diffusivity_v']
epsilon = fields_old['epsilon']
shear_freq2 = fields_old['shear_freq2']
buoy_freq2_neg = fields_old['buoy_freq2_neg']
buoy_freq2_pos = fields_old['buoy_freq2_pos']
p = eddy_viscosity * shear_freq2
b_source = - eddy_diffusivity * buoy_freq2_neg
b_sink = - eddy_diffusivity * buoy_freq2_pos
source = p + b_source + (b_sink - epsilon)/solution_old*solution # patankar
f = inner(source, self.test)*dx
return f
[docs]
class PsiSourceTerm(TracerTerm):
r"""
Production and destruction terms of the Psi equation :eq:`turb_psi_eq`
.. math::
F_\psi = \frac{\psi}{k} (c_1 P + c_3 B - c_2 \varepsilon f_{wall})
To ensure positivity we use Patankar-type time discretization: all source
terms are treated explicitly and sink terms are treated implicitly.
To this end the buoyancy production term :math:`c_3 B` is split in two:
.. math::
F_\psi = \frac{\psi^n}{k^n} (c_1 P + B_{source})
+ \frac{\psi^{n+1}}{k^n} (B_{sink} - c_2 \varepsilon f_{wall})
with :math:`B_{source} \ge 0` and :math:`B_{sink} < 0`.
Also implements Neumann boundary condition at top and bottom [7]
.. math::
\left( \frac{\nu}{\sigma_\psi} \frac{\psi}{z} \right)\Big|_{\Gamma_b} =
n_z \frac{\nu}{\sigma_\psi} (c_\mu^0)^p k^m \kappa^n (z_b + z_0)^{n-1}
where :math:`z_b` is the distance from boundary, and :math:`z_0` is the
roughness length.
"""
def __init__(self, function_space, gls_model,
bathymetry=None, v_elem_size=None, h_elem_size=None):
"""
:arg function_space: :class:`FunctionSpace` where the solution belongs
:arg gls_model: :class:`.GenericLengthScaleModel` object
:kwarg bathymetry: bathymetry of the domain
:type bathymetry: 3D :class:`Function` or :class:`Constant`
:kwarg v_elem_size: scalar :class:`Function` that defines the vertical
element size
:kwarg h_elem_size: scalar :class:`Function` that defines the horizontal
element size
"""
super(PsiSourceTerm, self).__init__(function_space,
bathymetry, v_elem_size, h_elem_size)
self.gls_model = gls_model
[docs]
def residual(self, solution, solution_old, fields, fields_old, bnd_conditions):
# psi: psi/k*(c1*P + c3*B - c2*eps*f_wall)
# P = viscosity M**2 (production)
# B = - diffusivity N**2 (byoyancy production)
# M**2 = (du/dz)**2 + (dv/dz)**2 (shear frequency)
# N**2 = -g\rho_0 (drho/dz) (buoyancy frequency)
# eps = (cmu0)**(3+p/n)*tke**(3/2+m/n)*psi**(-1/n)
# (tke dissipation rate)
eddy_viscosity = fields_old['viscosity_v']
eddy_diffusivity = fields_old['diffusivity_v']
diffusivity_v = eddy_viscosity/self.gls_model.options.schmidt_nb_psi
k = fields_old['k']
epsilon = fields_old['epsilon']
shear_freq2 = fields_old['shear_freq2']
buoy_freq2_neg = fields_old['buoy_freq2_neg']
buoy_freq2_pos = fields_old['buoy_freq2_pos']
p = eddy_viscosity * shear_freq2
c1 = self.gls_model.options.c1
c2 = self.gls_model.options.c2
# c3 switch: c3 = c3_minus if n2 > 0 else c3_plus
c3_minus = self.gls_model.options.c3_minus
c3_plus = self.gls_model.options.c3_plus # > 0
assert c3_plus >= 0, 'c3_plus has unexpected sign'
b_shear = c3_plus * -eddy_diffusivity * buoy_freq2_neg
b_buoy = c3_minus * -eddy_diffusivity * buoy_freq2_pos
if c3_minus > 0:
b_source = b_shear
b_sink = b_buoy
else:
b_source = b_shear + b_buoy
b_sink = 0
f_wall = self.gls_model.options.f_wall
source = solution_old/k*(c1*p + b_source) + solution/k*(b_sink - c2*f_wall*epsilon) # patankar
f = inner(source, self.test)*dx
# add bottom/top boundary condition for psi
# (nuv_v/sigma_psi * dpsi/dz)_b = n * nuv_v/sigma_psi * (cmu0)^p * k^m * kappa^n * z_b^(n-1)
# z_b = distance_from_bottom + z_0
cmu0 = self.gls_model.options.cmu0
p = self.gls_model.options.p
m = self.gls_model.options.m
n = self.gls_model.options.n
bfr_roughness = fields_old.get('bottom_roughness')
if bfr_roughness is None:
bfr_roughness = Constant(0)
kappa = physical_constants['von_karman']
if self.v_elem_size is None:
raise Exception('v_elem_size required')
# bottom condition
elem_frac = Constant(0.5)
z_b = elem_frac*self.v_elem_size + bfr_roughness
diff_flux = (n*diffusivity_v*(cmu0)**p
* k**m * kappa**n * z_b**(n - 1.0))
f += diff_flux*self.test*self.normal[2]*ds_bottom
# surface condition
elem_frac = Constant(0.5)
z0_surface = Constant(0.05) # TODO generalize
z_s = elem_frac*self.v_elem_size + z0_surface
diff_flux = -(n*diffusivity_v*(cmu0)**p
* k**m * kappa**n * z_s**(n - 1.0))
f += diff_flux*self.test*self.normal[2]*ds_surf
return f
[docs]
class GLSVerticalDiffusionTerm(VerticalDiffusionTerm):
"""
Vertical diffusion term where the diffusivity is replaced by
viscosity/Schmidt number.
"""
def __init__(self, function_space, schmidt_nb,
bathymetry=None, v_elem_size=None, h_elem_size=None,
sipg_factor=Constant(1.0)):
"""
:arg function_space: :class:`FunctionSpace` where the solution belongs
:arg schmidt_nb: the Schmidt number of TKE or Psi
:kwarg bathymetry: bathymetry of the domain
:type bathymetry: 3D :class:`Function` or :class:`Constant`
:kwarg v_elem_size: scalar :class:`Function` that defines the vertical
element size
:kwarg h_elem_size: scalar :class:`Function` that defines the horizontal
element size
:kwarg sipg_factor: :class: `Constant` or :class: `Function` vertical SIPG penalty scaling factor
"""
super(GLSVerticalDiffusionTerm, self).__init__(function_space,
bathymetry, v_elem_size, h_elem_size,
sipg_factor_vertical=sipg_factor)
self.schmidt_nb = schmidt_nb
[docs]
def residual(self, solution, solution_old, fields, fields_old, bnd_conditions):
d = {'diffusivity_v': fields_old['viscosity_v']/self.schmidt_nb}
f = super(GLSVerticalDiffusionTerm, self).residual(solution, solution_old,
d, d, bnd_conditions)
return f
[docs]
class TKEEquation(Equation):
"""
Turbulent kinetic energy equation :eq:`turb_tke_eq` without advection terms.
Advection of TKE is implemented using the standard tracer equation.
"""
def __init__(self, function_space, gls_model,
bathymetry=None, v_elem_size=None, h_elem_size=None):
"""
:arg function_space: :class:`FunctionSpace` where the solution belongs
:arg gls_model: :class:`.GenericLengthScaleModel` object
:kwarg bathymetry: bathymetry of the domain
:type bathymetry: 3D :class:`Function` or :class:`Constant`
:kwarg v_elem_size: scalar :class:`Function` that defines the vertical
element size
:kwarg h_elem_size: scalar :class:`Function` that defines the horizontal
element size
"""
super(TKEEquation, self).__init__(function_space)
diff = GLSVerticalDiffusionTerm(function_space,
gls_model.options.schmidt_nb_tke,
bathymetry, v_elem_size, h_elem_size)
source = TKESourceTerm(function_space,
gls_model,
bathymetry, v_elem_size, h_elem_size)
self.add_term(source, 'implicit')
self.add_term(diff, 'implicit')
[docs]
class PsiEquation(Equation):
r"""
Generic length scale equation :eq:`turb_psi_eq` without advection terms.
Advection of :math:`\psi` is implemented using the standard tracer equation.
"""
def __init__(self, function_space, gls_model,
bathymetry=None, v_elem_size=None, h_elem_size=None):
"""
:arg function_space: :class:`FunctionSpace` where the solution belongs
:arg gls_model: :class:`.GenericLengthScaleModel` object
:kwarg bathymetry: bathymetry of the domain
:type bathymetry: 3D :class:`Function` or :class:`Constant`
:kwarg v_elem_size: scalar :class:`Function` that defines the vertical
element size
:kwarg h_elem_size: scalar :class:`Function` that defines the horizontal
element size
"""
super(PsiEquation, self).__init__(function_space)
diff = GLSVerticalDiffusionTerm(function_space,
gls_model.options.schmidt_nb_psi,
bathymetry, v_elem_size, h_elem_size)
source = PsiSourceTerm(function_space,
gls_model,
bathymetry, v_elem_size, h_elem_size)
self.add_term(diff, 'implicit')
self.add_term(source, 'implicit')
[docs]
class PacanowskiPhilanderModel(TurbulenceModel):
r"""
Gradient Richardson number based model by Pacanowski and Philander (1981).
Given the gradient Richardson number :math:`Ri` the eddy viscosity and
diffusivity are computed diagnostically as
.. math::
\nu &= \frac{\nu_{max}}{(1 + \alpha Ri)^n} \\
\mu &= \frac{\nu}{1 + \alpha Ri}
where :math:`\nu_{max},\alpha,n` are constant parameters.
In unstably stratified cases where :math:`Ri<0`, value :math:`Ri=0` is used.
Pacanowski and Philander (1981). Parameterization of vertical mixing in
numerical models of tropical oceans. Journal of Physical Oceanography,
11(11):1443-1451.
http://dx.doi.org/10.1175/1520-0485(1981)011%3C1443:POVMIN%3E2.0.CO;2
"""
@PETSc.Log.EventDecorator("thetis.PacanowskiPhilanderModel.__init__")
def __init__(self, solver, uv_field, rho_field,
eddy_diffusivity, eddy_viscosity,
n2, m2, options=None):
"""
:arg solver: FlowSolver object
:arg uv_field: horizontal velocity field
:type uv_field: :class:`Function`
:arg rho_field: water density field
:type rho_field: :class:`Function`
:arg eddy_diffusivity: eddy diffusivity field
:type eddy_diffusivity: :class:`Function`
:arg eddy_viscosity: eddy viscosity field
:type eddy_viscosity: :class:`Function`
:arg n2: field for buoyancy frequency squared
:type n2: :class:`Function`
:arg m2: field for vertical shear frequency squared
:type m2: :class:`Function`
:kwarg options: model options
"""
self.solver = solver
# 3d model fields
self.uv = uv_field
self.rho = rho_field
# diagnostic fields
self.viscosity = eddy_viscosity
self.diffusivity = eddy_diffusivity
self.n2 = n2
self.m2 = m2
self.mu_tmp = Function(self.m2.function_space(),
name='tmp Shear frequency')
self.mu = Function(self.m2.function_space(), name='Shear frequency X')
self.mv = Function(self.m2.function_space(), name='Shear frequency Y')
self.n2_tmp = Function(self.n2.function_space(),
name='tmp buoyancy frequency')
self.n2_pos = Function(self.n2.function_space(),
name='positive buoyancy frequency')
self.options = PacanowskiPhilanderModelOptions()
if options is not None:
self.options.update(options)
self.shear_frequency_solver = ShearFrequencySolver(self.uv, self.m2,
self.mu, self.mv,
self.mu_tmp)
if self.rho is not None:
self.buoy_frequency_solver = BuoyFrequencySolver(self.rho, self.n2,
self.n2_tmp)
self.initialize()
print_output(self.options)
[docs]
@PETSc.Log.EventDecorator("thetis.PacanowskiPhilanderModel.initialize")
def initialize(self):
"""Initializes fields"""
self.n2.assign(1e-12)
self.n2_pos.assign(1e-12)
self.preprocess(init_solve=True)
self.postprocess()
[docs]
@PETSc.Log.EventDecorator("thetis.PacanowskiPhilanderModel.preprocess")
def preprocess(self, init_solve=False):
"""
Computes all diagnostic variables that depend on the mean flow model
variables.
To be called before updating the turbulence PDEs.
"""
self.shear_frequency_solver.solve(init_solve=init_solve)
if self.rho is not None:
self.buoy_frequency_solver.solve(init_solve=init_solve)
# split to positive and negative parts
self.n2_pos.assign(1e-12)
pos_ix = self.n2.dat.data[:] >= 0.0
self.n2_pos.dat.data[pos_ix] = self.n2.dat.data[pos_ix]
[docs]
@PETSc.Log.EventDecorator("thetis.PacanowskiPhilanderModel.postprocess")
def postprocess(self):
"""
Updates all diagnostic variables that depend on the turbulence state
variables.
To be called after evaluating the equations.
"""
ri = self.n2_pos.dat.data[:]/self.m2.dat.data[:]
denom = 1.0 + self.options.alpha*ri
self.viscosity.dat.data[:] = self.options.max_viscosity/denom**self.options.exponent
self.diffusivity.dat.data[:] = self.viscosity.dat.data[:]/denom
