r"""
Depth averaged shallow water equations
---------
Equations
---------
The state variables are water elevation, :math:`\eta`, and depth averaged
velocity :math:`\bar{\textbf{u}}`.
Denoting the total water depth by :math:`H=\eta + h`, the non-conservative form of
the free surface equation is
.. math::
\frac{\partial \eta}{\partial t} + \nabla \cdot (H \bar{\textbf{u}}) = 0
:label: swe_freesurf
The non-conservative momentum equation reads
.. math::
\frac{\partial \bar{\textbf{u}}}{\partial t} +
\bar{\textbf{u}} \cdot \nabla\bar{\textbf{u}} +
f\textbf{e}_z\wedge \bar{\textbf{u}} +
g \nabla \eta +
\nabla \left(\frac{p_a}{\rho_0} \right) +
g \frac{1}{H}\int_{-h}^\eta \nabla r dz
= \nabla \cdot \big( \nu_h ( \nabla \bar{\textbf{u}} + (\nabla \bar{\textbf{u}})^T )\big) +
\frac{\nu_h \nabla(H)}{H} \cdot ( \nabla \bar{\textbf{u}} + (\nabla \bar{\textbf{u}})^T ),
:label: swe_momentum
where :math:`g` is the gravitational acceleration, :math:`f` is the Coriolis
frequency, :math:`\wedge` is the cross product, :math:`\textbf{e}_z` is a vertical unit vector,
:math:`p_a` is the atmospheric pressure at the free surface, and :math:`\nu_h`
is viscosity. Water density is given by :math:`\rho = \rho'(T, S, p) + \rho_0`,
where :math:`\rho_0` is a constant reference density.
Above :math:`r` denotes the baroclinic head
.. math::
r = \frac{1}{\rho_0} \int_{z}^\eta \rho' d\zeta.
In the case of purely barotropic problems the :math:`r` and the internal pressure
gradient are omitted.
If the option :attr:`.ModelOptions.use_nonlinear_equations` is ``False``, we solve the linear shallow water
equations (i.e. wave equation):
.. math::
\frac{\partial \eta}{\partial t} + \nabla \cdot (h \bar{\textbf{u}}) = 0
:label: swe_freesurf_linear
.. math::
\frac{\partial \bar{\textbf{u}}}{\partial t} +
f\textbf{e}_z\wedge \bar{\textbf{u}} +
g \nabla \eta
= \nabla \cdot \big( \nu_h ( \nabla \bar{\textbf{u}} + (\nabla \bar{\textbf{u}})^T )\big) +
\frac{\nu_h \nabla(H)}{H} \cdot ( \nabla \bar{\textbf{u}} + (\nabla \bar{\textbf{u}})^T ).
:label: swe_momentum_linear
In case of a 3D problem with mode splitting, we use a simplified 2D
system that contains nothing but the rotational external gravity waves:
.. math::
\frac{\partial \eta}{\partial t} + \nabla \cdot (H \bar{\textbf{u}}) = 0
:label: swe_freesurf_modesplit
.. math::
\frac{\partial \bar{\textbf{u}}}{\partial t} +
f\textbf{e}_z\wedge \bar{\textbf{u}} +
g \nabla \eta
= \textbf{G},
:label: swe_momentum_modesplit
where :math:`\textbf{G}` is a source term used to couple the 2D and 3D momentum
equations.
-------------------
Boundary Conditions
-------------------
All boundary conditions are imposed weakly by providing external values for
:math:`\eta` and :math:`\bar{\textbf{u}}`.
Boundary conditions are set with a dictionary that defines all prescribed
variables at each open boundary.
For example, to assign elevation and volume flux on boundary ``1`` we set
.. code-block:: python
swe_bnd_funcs = {}
swe_bnd_funcs[1] = {'elev':myfunc1, 'flux':myfunc2}
where ``myfunc1`` and ``myfunc2`` are :class:`Constant` or :class:`Function`
objects.
The user can provide :math:`\eta` and/or :math:`\bar{\textbf{u}}` values.
Supported combinations are:
- *unspecified* : impermeable (land) boundary, implies symmetric :math:`\eta` condition and zero normal velocity
- ``'elev'``: elevation only, symmetric velocity (usually unstable)
- ``'uv'``: 2d velocity vector :math:`\bar{\textbf{u}}=(u, v)` (in mesh coordinates), symmetric elevation
- ``'un'``: normal velocity (scalar, positive out of domain), symmetric elevation
- ``'flux'``: normal volume flux (scalar, positive out of domain), symmetric elevation
- ``'elev'`` and ``'uv'``: water elevation and 2d velocity vector
- ``'elev'`` and ``'un'``: water elevation and normal velocity
- ``'elev'`` and ``'flux'``: water elevation and normal flux
The boundary conditions are assigned to the :class:`.FlowSolver2d` or
:class:`.FlowSolver` objects:
.. code-block:: python
solver_obj = solver2d.FlowSolver2d(...)
...
solver_obj.bnd_functions['shallow_water'] = swe_bnd_funcs
Internally the boundary conditions passed to the :meth:`.Term.residual` method
of each term:
.. code-block:: python
adv_term = shallowwater_eq.HorizontalAdvectionTerm(...)
adv_form = adv_term.residual(..., bnd_conditions=swe_bnd_funcs)
------------------
Wetting and drying
------------------
If the option :attr:`.ModelOptions.use_wetting_and_drying` is ``True``, then wetting and
drying is included through the formulation of Karna et al. (2011).
The method introduces a modified bathymetry :math:`\tilde{h} = h + f(H)`, which ensures
positive total water depth, with :math:`f(H)` defined by
.. math::
f(H) = \frac{1}{2}(\sqrt{H^2 + \alpha^2} - H),
introducing a wetting-drying parameter :math:`\alpha`, with dimensions of length. This
results in a modified total water depth :math:`\tilde{H}=H+f(H)`.
The value for :math:`\alpha` is specified by the user through the
option :attr:`.ModelOptions.wetting_and_drying_alpha`, in units of meters. The default value
for :attr:`.ModelOptions.wetting_and_drying_alpha` is 0.5, but the appropriate value is problem
specific and should be set by the user.
An approximate method for selecting a suitable value for :math:`\alpha` is suggested
by Karna et al. (2011). Defining :math:`L_x` as the horizontal length scale of the
mesh elements at the wet-dry front, it can be reasoned that :math:`\alpha \approx |L_x
\nabla h|` yields a suitable choice. Smaller :math:`\alpha` leads to a more accurate
solution to the shallow water equations in wet regions, but if :math:`\alpha` is too
small the simulation will become unstable.
When wetting and drying is turned on, two things occur:
1. All instances of the height, :math:`H`, are replaced by :math:`\tilde{H}` (as defined above);
2. An additional displacement term :math:`\frac{\partial \tilde{h}}{\partial t}` is added to the bathymetry in the free surface equation.
The free surface and momentum equations then become:
.. math::
\frac{\partial \eta}{\partial t} + \frac{\partial \tilde{h}}{\partial t} + \nabla \cdot (\tilde{H} \bar{\textbf{u}}) = 0,
:label: swe_freesurf_wd
.. math::
\frac{\partial \bar{\textbf{u}}}{\partial t} +
\bar{\textbf{u}} \cdot \nabla\bar{\textbf{u}} +
f\textbf{e}_z\wedge \bar{\textbf{u}} +
g \nabla \eta +
g \frac{1}{\tilde{H}}\int_{-h}^\eta \nabla r dz
= \nabla \cdot \big( \nu_h ( \nabla \bar{\textbf{u}} + (\nabla \bar{\textbf{u}})^T )\big) +
\frac{\nu_h \nabla(\tilde{H})}{\tilde{H}} \cdot ( \nabla \bar{\textbf{u}} + (\nabla \bar{\textbf{u}})^T ).
:label: swe_momentum_wd
"""
from .utility import *
from .equation import Term, Equation
__all__ = [
'BaseShallowWaterEquation',
'ShallowWaterEquations',
'ModeSplit2DEquations',
'ShallowWaterMomentumEquation',
'FreeSurfaceEquation',
'ShallowWaterTerm',
'ShallowWaterMomentumTerm',
'ShallowWaterContinuityTerm',
'HUDivTerm',
'ContinuitySourceTerm',
'HorizontalAdvectionTerm',
'HorizontalViscosityTerm',
'ExternalPressureGradientTerm',
'CoriolisTerm',
'LinearDragTerm',
'QuadraticDragTerm',
'BottomDrag3DTerm',
'MomentumSourceTerm',
'WindStressTerm',
'AtmosphericPressureTerm',
]
g_grav = physical_constants['g_grav']
rho_0 = physical_constants['rho0']
[docs]
class ShallowWaterTerm(Term):
"""
Generic term in the shallow water equations that provides commonly used
members and mapping for boundary functions.
"""
def __init__(self, space,
depth, options=None):
super(ShallowWaterTerm, self).__init__(space)
self.depth = depth
self.options = options
# mesh dependent variables
self.cellsize = CellSize(self.mesh)
assert self.mesh.cell_dimension() == 2
self.on_the_sphere = self.mesh.geometric_dimension() == 3
# define measures with a reasonable quadrature degree
p = self.function_space.ufl_element().degree()
self.quad_degree = 2*p + 1
self.dx = dx(degree=self.quad_degree,
domain=self.function_space.ufl_domain())
self.dS = dS(degree=self.quad_degree,
domain=self.function_space.ufl_domain())
[docs]
def get_bnd_functions(self, eta_in, uv_in, bnd_id, bnd_conditions):
"""
Returns external values of elev and uv for all supported
boundary conditions.
Volume flux (flux) and normal velocity (un) are defined positive out of
the domain.
"""
bnd_len = self.boundary_len[bnd_id]
funcs = bnd_conditions[bnd_id]
eta_ext = uv_ext = None
if 'elev' in funcs and 'uv' in funcs:
eta_ext = funcs['elev']
uv_ext = funcs['uv']
elif 'elev' in funcs and 'un' in funcs:
eta_ext = funcs['elev']
uv_ext = funcs['un']*self.normal
elif 'elev' in funcs and 'flux' in funcs:
eta_ext = funcs['elev']
h_ext = self.depth.get_total_depth(eta_ext)
area = h_ext*bnd_len # NOTE using external data only
uv_ext = funcs['flux']/area*self.normal
elif 'elev' in funcs:
eta_ext = funcs['elev']
uv_ext = uv_in # assume symmetry
elif 'uv' in funcs:
eta_ext = eta_in # assume symmetry
uv_ext = funcs['uv']
elif 'un' in funcs:
eta_ext = eta_in # assume symmetry
uv_ext = funcs['un']*self.normal
elif 'flux' in funcs:
eta_ext = eta_in # assume symmetry
h_ext = self.depth.get_total_depth(eta_ext)
area = h_ext*bnd_len # NOTE using internal elevation
uv_ext = funcs['flux']/area*self.normal
if eta_ext is None or uv_ext is None:
raise Exception('Unsupported bnd type, one of "elev", "uv", "un", '
'or "flux" must be defined on boundary '
f'{bnd_marker}: {funcs}')
return eta_ext, uv_ext
[docs]
def impose_dynamic_bnd(self, bnd_funcs, bnd_id):
"""
Check if the prognostic variables have been specified on the boundary.
If any prognostic value has been specified, dynamic boundary terms
will be evaluated. If not, a closed boundary is assumed (which implies
e.g. that volume flux boundary terms are omitted).
:arg bnd_funcs: None or dictionary of boundary coefficients.
:arg bnd_id: Boundary marker id
:returns: True if elev or uv is defined, otherwise False.
"""
open_tags = ['elev', 'uv', 'un', 'flux']
all_tags = open_tags + ['drag']
if bnd_funcs is None:
return False
for k in bnd_funcs.keys():
if k not in all_tags:
raise Exception(f'Invalid boundary tag "{k}" '
f'specified on boundary {bnd_id}')
if k in open_tags:
return True
return False
[docs]
class ShallowWaterMomentumTerm(ShallowWaterTerm):
"""
Generic term in the shallow water momentum equation that provides commonly used
members and mapping for boundary functions.
"""
def __init__(self, u_test, u_space, eta_space,
depth, options=None):
super(ShallowWaterMomentumTerm, self).__init__(u_space, depth, options)
self.options = options
self.u_test = u_test
self.u_space = u_space
self.eta_space = eta_space
self.u_continuity = element_continuity(self.u_space.ufl_element()).horizontal
self.eta_is_dg = element_continuity(self.eta_space.ufl_element()).horizontal == 'dg'
[docs]
class ShallowWaterContinuityTerm(ShallowWaterTerm):
"""
Generic term in the depth-integrated continuity equation that provides commonly used
members and mapping for boundary functions.
"""
def __init__(self, eta_test, eta_space, u_space,
depth, options=None):
super(ShallowWaterContinuityTerm, self).__init__(eta_space, depth, options)
self.eta_test = eta_test
self.eta_space = eta_space
self.u_space = u_space
self.u_continuity = element_continuity(self.u_space.ufl_element()).horizontal
self.eta_is_dg = element_continuity(self.eta_space.ufl_element()).horizontal == 'dg'
[docs]
class ExternalPressureGradientTerm(ShallowWaterMomentumTerm):
r"""
External pressure gradient term, :math:`g \nabla \eta`
The weak form reads
.. math::
\int_\Omega g \nabla \eta \cdot \boldsymbol{\psi} dx
= \int_\Gamma g \eta^* \text{jump}(\boldsymbol{\psi} \cdot \textbf{n}) dS
- \int_\Omega g \eta \nabla \cdot \boldsymbol{\psi} dx
where the right hand side has been integrated by parts; :math:`\textbf{n}`
denotes the unit normal of the element interfaces, :math:`n^*` is value at
the interface obtained from an approximate Riemann solver.
If :math:`\eta` belongs to a discontinuous function space, the form on the
right hand side is used.
"""
[docs]
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
total_h = self.depth.get_total_depth(eta_old)
head = eta
grad_eta_by_parts = self.eta_is_dg
if grad_eta_by_parts:
f = -g_grav*head*nabla_div(self.u_test)*self.dx
if uv is not None:
head_star = avg(head) + sqrt(avg(total_h)/g_grav)*jump(uv, self.normal)
else:
head_star = avg(head)
f += g_grav*head_star*jump(self.u_test, self.normal)*self.dS
for bnd_marker in self.boundary_markers:
funcs = bnd_conditions.get(bnd_marker)
ds_bnd = ds(int(bnd_marker), degree=self.quad_degree)
if self.impose_dynamic_bnd(funcs, bnd_marker):
eta_ext, uv_ext = self.get_bnd_functions(head, uv, bnd_marker, bnd_conditions)
# Compute linear riemann solution with eta, eta_ext, uv, uv_ext
un_jump = inner(uv - uv_ext, self.normal)
eta_rie = 0.5*(head + eta_ext) + sqrt(total_h/g_grav)*un_jump
f += g_grav*eta_rie*dot(self.u_test, self.normal)*ds_bnd
else:
# assume land boundary
# impermeability implies external un=0
un_jump = inner(uv, self.normal)
head_rie = head + sqrt(total_h/g_grav)*un_jump
f += g_grav*head_rie*dot(self.u_test, self.normal)*ds_bnd
else:
f = g_grav*inner(grad(head), self.u_test) * self.dx
for bnd_marker in self.boundary_markers:
funcs = bnd_conditions.get(bnd_marker)
ds_bnd = ds(int(bnd_marker), degree=self.quad_degree)
if self.impose_dynamic_bnd(funcs, bnd_marker):
eta_ext, uv_ext = self.get_bnd_functions(head, uv, bnd_marker, bnd_conditions)
# Compute linear riemann solution with eta, eta_ext, uv, uv_ext
un_jump = inner(uv - uv_ext, self.normal)
eta_rie = 0.5*(head + eta_ext) + sqrt(total_h/g_grav)*un_jump
f += g_grav*(eta_rie-head)*dot(self.u_test, self.normal)*ds_bnd
return -f
[docs]
class HUDivTerm(ShallowWaterContinuityTerm):
r"""
Divergence term, :math:`\nabla \cdot (H \bar{\textbf{u}})`
The weak form reads
.. math::
\int_\Omega \nabla \cdot (H \bar{\textbf{u}}) \phi dx
= \int_\Gamma (H^* \bar{\textbf{u}}^*) \cdot \text{jump}(\phi \textbf{n}) dS
- \int_\Omega H (\bar{\textbf{u}}\cdot\nabla \phi) dx
where the right hand side has been integrated by parts; :math:`\textbf{n}`
denotes the unit normal of the element interfaces, and :math:`\text{jump}`
and :math:`\text{avg}` denote the jump and average operators across the
interface. :math:`H^*, \bar{\textbf{u}}^*` are values at the interface
obtained from an approximate Riemann solver.
If :math:`\bar{\textbf{u}}` belongs to a discontinuous function space,
the form on the right hand side is used.
"""
[docs]
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
total_h = self.depth.get_total_depth(eta_old)
hu_by_parts = self.u_continuity in ['dg', 'hdiv']
if hu_by_parts:
f = -inner(grad(self.eta_test), total_h*uv)*self.dx
if self.eta_is_dg:
h = avg(total_h)
uv_rie = avg(uv) + sqrt(g_grav/h)*jump(eta, self.normal)
hu_star = h*uv_rie
f += inner(jump(self.eta_test, self.normal), hu_star)*self.dS
for bnd_marker in self.boundary_markers:
funcs = bnd_conditions.get(bnd_marker)
ds_bnd = ds(int(bnd_marker), degree=self.quad_degree)
if self.impose_dynamic_bnd(funcs, bnd_marker):
eta_ext, uv_ext = self.get_bnd_functions(eta, uv, bnd_marker, bnd_conditions)
eta_ext_old, uv_ext_old = self.get_bnd_functions(eta_old, uv_old, bnd_marker, bnd_conditions)
# Compute linear riemann solution with eta, eta_ext, uv, uv_ext
total_h_ext = self.depth.get_total_depth(eta_ext_old)
h_av = 0.5*(total_h + total_h_ext)
eta_jump = eta - eta_ext
un_rie = 0.5*inner(uv + uv_ext, self.normal) + sqrt(g_grav/h_av)*eta_jump
un_jump = inner(uv_old - uv_ext_old, self.normal)
eta_rie = 0.5*(eta_old + eta_ext_old) + sqrt(h_av/g_grav)*un_jump
h_rie = self.depth.get_total_depth(eta_rie)
f += h_rie*un_rie*self.eta_test*ds_bnd
else:
f = div(total_h*uv)*self.eta_test*self.dx
for bnd_marker in self.boundary_markers:
funcs = bnd_conditions.get(bnd_marker)
ds_bnd = ds(int(bnd_marker), degree=self.quad_degree)
if not self.impose_dynamic_bnd(funcs, bnd_marker) or 'un' in funcs:
f += -total_h*dot(uv, self.normal)*self.eta_test*ds_bnd
return -f
[docs]
class HorizontalAdvectionTerm(ShallowWaterMomentumTerm):
r"""
Advection of momentum term, :math:`\bar{\textbf{u}} \cdot \nabla\bar{\textbf{u}}`
The weak form is
.. math::
\int_\Omega \bar{\textbf{u}} \cdot \nabla\bar{\textbf{u}} \cdot \boldsymbol{\psi} dx
= - \int_\Omega \nabla_h \cdot (\bar{\textbf{u}} \boldsymbol{\psi}) \cdot \bar{\textbf{u}} dx
+ \int_\Gamma \text{avg}(\bar{\textbf{u}}) \cdot \text{jump}(\boldsymbol{\psi}
(\bar{\textbf{u}}\cdot\textbf{n})) dS
where the right hand side has been integrated by parts;
:math:`\textbf{n}` is the unit normal of
the element interfaces, and :math:`\text{jump}` and :math:`\text{avg}` denote the
jump and average operators across the interface.
"""
[docs]
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
if not self.options.use_nonlinear_equations:
return 0
horiz_advection_by_parts = True
if horiz_advection_by_parts:
f = -inner(div(outer(self.u_test, uv_old)), uv)*self.dx
if self.u_continuity in ['dg', 'hdiv']:
un_av = dot(avg(uv_old), self.normal('-'))
# NOTE mean flux
uv_avg = avg(uv)
f += inner(uv_avg, jump(outer(self.u_test, uv_old), self.normal))*self.dS
# Lax-Friedrichs stabilization
if self.options.use_lax_friedrichs_velocity:
uv_lax_friedrichs = fields_old.get('lax_friedrichs_velocity_scaling_factor')
gamma = 0.5*abs(un_av)*uv_lax_friedrichs
f += gamma*dot(jump(self.u_test), jump(uv))*self.dS
for bnd_marker in self.boundary_markers:
funcs = bnd_conditions.get(bnd_marker)
ds_bnd = ds(int(bnd_marker), degree=self.quad_degree)
if not self.impose_dynamic_bnd(funcs, bnd_marker):
# impose impermeability with mirror velocity
n = self.normal
uv_ext = uv - 2*dot(uv, n)*n
gamma = 0.5*abs(dot(uv_old, n))*uv_lax_friedrichs
f += gamma*dot(self.u_test, uv-uv_ext)*ds_bnd
for bnd_marker in self.boundary_markers:
funcs = bnd_conditions.get(bnd_marker)
ds_bnd = ds(int(bnd_marker), degree=self.quad_degree)
if self.impose_dynamic_bnd(funcs, bnd_marker):
eta_ext, uv_ext = self.get_bnd_functions(eta, uv, bnd_marker, bnd_conditions)
eta_ext_old, uv_ext_old = self.get_bnd_functions(eta_old, uv_old, bnd_marker, bnd_conditions)
# Compute linear riemann solution with eta, eta_ext, uv, uv_ext
eta_jump = eta_old - eta_ext_old
total_h = self.depth.get_total_depth(eta_old)
un_rie = 0.5*inner(uv_old + uv_ext_old, self.normal) + sqrt(g_grav/total_h)*eta_jump
uv_av = 0.5*(uv_ext + uv)
f += un_rie*inner(uv_av, self.u_test)*ds_bnd
return -f
[docs]
class HorizontalViscosityTerm(ShallowWaterMomentumTerm):
r"""
Viscosity of momentum term
If option :attr:`.ModelOptions.use_grad_div_viscosity_term` is ``True``, we
use the symmetric viscous stress :math:`\boldsymbol{\tau}_\nu = \nu_h ( \nabla \bar{\textbf{u}} + (\nabla \bar{\textbf{u}})^T )`.
Using the symmetric interior penalty method the weak form then reads
.. math::
\int_\Omega -\nabla \cdot \boldsymbol{\tau}_\nu \cdot \boldsymbol{\psi} dx
=& \int_\Omega (\nabla \boldsymbol{\psi}) : \boldsymbol{\tau}_\nu dx \\
&- \int_\Gamma \text{jump}(\boldsymbol{\psi} \textbf{n}) \cdot \text{avg}(\boldsymbol{\tau}_\nu) dS
- \int_\Gamma \text{avg}(\nu_h)\big(\text{jump}(\bar{\textbf{u}} \textbf{n}) + \text{jump}(\bar{\textbf{u}} \textbf{n})^T\big) \cdot \text{avg}(\nabla \boldsymbol{\psi}) dS \\
&+ \int_\Gamma \sigma \text{avg}(\nu_h) \big(\text{jump}(\bar{\textbf{u}} \textbf{n}) + \text{jump}(\bar{\textbf{u}} \textbf{n})^T\big) \cdot \text{jump}(\boldsymbol{\psi} \textbf{n}) dS
where :math:`\sigma` is a penalty parameter, see Hillewaert (2013).
If option :attr:`.ModelOptions.use_grad_div_viscosity_term` is ``False``,
we use viscous stress :math:`\boldsymbol{\tau}_\nu = \nu_h \nabla \bar{\textbf{u}}`.
In this case the weak form is
.. math::
\int_\Omega -\nabla \cdot \boldsymbol{\tau}_\nu \cdot \boldsymbol{\psi} dx
=& \int_\Omega (\nabla \boldsymbol{\psi}) : \boldsymbol{\tau}_\nu dx \\
&- \int_\Gamma \text{jump}(\boldsymbol{\psi} \textbf{n}) \cdot \text{avg}(\boldsymbol{\tau}_\nu) dS
- \int_\Gamma \text{avg}(\nu_h)\text{jump}(\bar{\textbf{u}} \textbf{n}) \cdot \text{avg}(\nabla \boldsymbol{\psi}) dS \\
&+ \int_\Gamma \sigma \text{avg}(\nu_h) \text{jump}(\bar{\textbf{u}} \textbf{n}) \cdot \text{jump}(\boldsymbol{\psi} \textbf{n}) dS
If option :attr:`.ModelOptions.use_grad_depth_viscosity_term` is ``True``, we also include
the term
.. math::
\boldsymbol{\tau}_{\nabla H} = - \frac{\nu_h \nabla(H)}{H} \cdot ( \nabla \bar{\textbf{u}} + (\nabla \bar{\textbf{u}})^T )
as a source term.
Hillewaert, Koen (2013). Development of the discontinuous Galerkin method
for high-resolution, large scale CFD and acoustics in industrial
geometries. PhD Thesis. Université catholique de Louvain.
https://dial.uclouvain.be/pr/boreal/object/boreal:128254/
"""
[docs]
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
total_h = self.depth.get_total_depth(eta_old)
nu = fields_old.get('viscosity_h')
sipg_factor = self.options.sipg_factor
if nu is None:
return 0
n = self.normal
if self.options.use_grad_div_viscosity_term:
stress = nu*2.*sym(grad(uv))
stress_jump = avg(nu)*2.*sym(tensor_jump(uv, n))
else:
stress = nu*grad(uv)
stress_jump = avg(nu)*tensor_jump(uv, n)
f = inner(grad(self.u_test), stress)*self.dx
if self.u_continuity in ['dg', 'hdiv']:
cell = self.mesh.ufl_cell()
p = self.function_space.ufl_element().degree()
cp = (p + 1) * (p + 2) / 2 if cell == triangle else (p + 1)**2
l_normal = CellVolume(self.mesh) / FacetArea(self.mesh)
# by default the factor is multiplied by 2 to ensure convergence
sigma = sipg_factor * cp / l_normal
sp = sigma('+')
sm = sigma('-')
sigma_max = conditional(sp > sm, sp, sm)
f += (
+ sigma_max*inner(tensor_jump(self.u_test, n), stress_jump)*self.dS
- inner(avg(grad(self.u_test)), stress_jump)*self.dS
- inner(tensor_jump(self.u_test, n), avg(stress))*self.dS
)
# Dirichlet bcs only for DG
for bnd_marker in self.boundary_markers:
funcs = bnd_conditions.get(bnd_marker)
ds_bnd = ds(int(bnd_marker), degree=self.quad_degree)
if self.impose_dynamic_bnd(funcs, bnd_marker):
if 'un' in funcs:
delta_uv = (dot(uv, n) - funcs['un'])*n
else:
eta_ext, uv_ext = self.get_bnd_functions(eta, uv, bnd_marker, bnd_conditions)
if uv_ext is uv:
continue
delta_uv = uv - uv_ext
if self.options.use_grad_div_viscosity_term:
stress_jump = nu*2.*sym(outer(delta_uv, n))
else:
stress_jump = nu*outer(delta_uv, n)
f += (
sigma*inner(outer(self.u_test, n), stress_jump)*ds_bnd
- inner(grad(self.u_test), stress_jump)*ds_bnd
- inner(outer(self.u_test, n), stress)*ds_bnd
)
if self.options.use_grad_depth_viscosity_term:
f += -dot(self.u_test, dot(grad(total_h)/total_h, stress))*self.dx
return -f
[docs]
class CoriolisTerm(ShallowWaterMomentumTerm):
r"""
Coriolis term, :math:`f\textbf{e}_z\wedge \bar{\textbf{u}}`
"""
[docs]
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
coriolis = fields_old.get('coriolis')
f = 0
if coriolis is not None:
if self.on_the_sphere:
outward_normals = CellNormal(self.mesh)
u = cross(outward_normals, uv)
f = coriolis*inner(self.u_test, u)*self.dx
else:
ez_x_uv = -uv[1]*self.u_test[0] + uv[0]*self.u_test[1]
f += coriolis*ez_x_uv*self.dx
return -f
[docs]
class WindStressTerm(ShallowWaterMomentumTerm):
r"""
Wind stress term, :math:`-\tau_w/(H \rho_0)`
Here :math:`\tau_w` is a user-defined wind stress :class:`Function`.
"""
[docs]
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
wind_stress = fields_old.get('wind_stress')
total_h = self.depth.get_total_depth(eta_old)
f = 0
if wind_stress is not None:
f += dot(wind_stress, self.u_test)/total_h/rho_0*self.dx
return f
[docs]
class AtmosphericPressureTerm(ShallowWaterMomentumTerm):
r"""
Atmospheric pressure term, :math:`\nabla (p_a / \rho_0)`
Here :math:`p_a` is a user-defined atmospheric pressure :class:`Function`.
"""
[docs]
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
atmospheric_pressure = fields_old.get('atmospheric_pressure')
f = 0
if atmospheric_pressure is not None:
f += dot(grad(atmospheric_pressure), self.u_test)/rho_0*self.dx
return -f
[docs]
class QuadraticDragTerm(ShallowWaterMomentumTerm):
r"""
Quadratic Manning bottom friction term
:math:`C_D \| \bar{\textbf{u}} \| \bar{\textbf{u}}`
where the drag coefficient is computed with the Manning formula
.. math::
C_D = g \frac{\mu^2}{H^{1/3}}
if the Manning coefficient :math:`\mu` is defined (see field :attr:`manning_drag_coefficient`).
Otherwise :math:`C_D` is taken as a constant (see field :attr:`quadratic_drag_coefficient`).
"""
[docs]
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
total_h = self.depth.get_total_depth(eta_old)
manning_drag_coefficient = fields_old.get('manning_drag_coefficient')
nikuradse_bed_roughness = fields_old.get('nikuradse_bed_roughness')
C_D = fields_old.get('quadratic_drag_coefficient')
f = 0
if manning_drag_coefficient is not None:
if C_D is not None:
raise Exception('Cannot set both dimensionless and Manning drag parameter')
C_D = g_grav * manning_drag_coefficient**2 / total_h**(1./3.)
if nikuradse_bed_roughness is not None:
if manning_drag_coefficient is not None:
raise Exception('Cannot set both Nikuradse drag and Manning drag parameter')
if C_D is not None:
raise Exception('Cannot set both dimensionless and Nikuradse drag parameter')
kappa = physical_constants['von_karman']
C_D = conditional(total_h > nikuradse_bed_roughness, 2*(kappa**2)/(ln(11.036*total_h/nikuradse_bed_roughness)**2), Constant(0.0))
if C_D is not None:
f += C_D * sqrt(dot(uv_old, uv_old) + self.options.norm_smoother**2) * inner(self.u_test, uv) / total_h * self.dx
return -f
class BoundaryDragTerm(ShallowWaterMomentumTerm):
r"""
Quadratic friction term on the boundary
:math:`C_D \| \bar{\textbf{u}}_t \| \bar{\textbf{u}}_t`
where :math:`\bar{\textbf{u}}_t` denotes the tangential velocity component
and the drag coefficient :math:`C_D` is user-defined.
"""
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
f = 0
for bnd_marker in self.boundary_markers:
funcs = bnd_conditions.get(bnd_marker)
ds_bnd = ds(int(bnd_marker), degree=self.quad_degree)
if funcs is not None and 'drag' in funcs:
C_D = funcs['drag']
# compute tangential velocity
n = self.normal
ut = uv - dot(uv, n) * n
ut_old = uv_old - dot(uv_old, n) * n
ut_mag = sqrt(dot(ut_old, ut_old))
f += C_D * ut_mag * inner(self.u_test, ut) * ds_bnd
return -f
[docs]
class LinearDragTerm(ShallowWaterMomentumTerm):
r"""
Linear friction term, :math:`C \bar{\textbf{u}}`
Here :math:`C` is a user-defined drag coefficient.
"""
[docs]
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
linear_drag_coefficient = fields_old.get('linear_drag_coefficient')
f = 0
if linear_drag_coefficient is not None:
bottom_fri = linear_drag_coefficient*inner(self.u_test, uv)*self.dx
f += bottom_fri
return -f
[docs]
class BottomDrag3DTerm(ShallowWaterMomentumTerm):
r"""
Bottom drag term consistent with the 3D mode,
:math:`C_D \| \textbf{u}_b \| \textbf{u}_b`
Here :math:`\textbf{u}_b` is the bottom velocity used in the 3D mode, and
:math:`C_D` the corresponding bottom drag.
These fields are computed in the 3D model.
"""
[docs]
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
total_h = self.depth.get_total_depth(eta_old)
bottom_drag = fields_old.get('bottom_drag')
uv_bottom = fields_old.get('uv_bottom')
f = 0
if bottom_drag is not None and uv_bottom is not None:
uvb_mag = sqrt(uv_bottom[0]**2 + uv_bottom[1]**2)
stress = bottom_drag*uvb_mag*uv_bottom/total_h
bot_friction = dot(stress, self.u_test)*self.dx
f += bot_friction
return -f
class TurbineDragTerm(ShallowWaterMomentumTerm):
r"""
Turbine drag parameterisation implemented through quadratic drag term
:math:`c_t \| \bar{\textbf{u}} \| \bar{\textbf{u}}`
where the turbine drag :math:`c_t` is related to the turbine thrust coefficient
:math:`C_T`, the turbine diameter :math:`A_T`, and the turbine density :math:`d`
(n/o turbines per unit area), by:
.. math::
c_t = (C_T A_T d)/2
"""
def __init__(self, u_test, u_space, eta_space,
depth, options=None, tidal_farms=None):
super().__init__(u_test, u_space, eta_space, depth, options=options)
self.tidal_farms = tidal_farms
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
total_h = self.depth.get_total_depth(eta_old)
f = 0
for farm in self.tidal_farms:
density = farm.turbine_density
c_t = farm.friction_coefficient(uv_old, total_h)
unorm = sqrt(dot(uv_old, uv_old))
f += c_t * density * unorm * inner(self.u_test, uv) / total_h * farm.dx
return -f
[docs]
class MomentumSourceTerm(ShallowWaterMomentumTerm):
r"""
Generic source term in the shallow water momentum equation
The weak form reads
.. math::
F_s = \int_\Omega \boldsymbol{\tau} \cdot \boldsymbol{\psi} dx
where :math:`\boldsymbol{\tau}` is a user defined vector valued :class:`Function`.
"""
[docs]
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
f = 0
momentum_source = fields_old.get('momentum_source')
if momentum_source is not None:
f += inner(momentum_source, self.u_test)*self.dx
return f
[docs]
class ContinuitySourceTerm(ShallowWaterContinuityTerm):
r"""
Generic source term in the depth-averaged continuity equation
The weak form reads
.. math::
F_s = \int_\Omega S \phi dx
where :math:`S` is a user defined scalar :class:`Function`.
"""
[docs]
def residual(self, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
f = 0
volume_source = fields_old.get('volume_source')
if volume_source is not None:
f += inner(volume_source, self.eta_test)*self.dx
return f
class BathymetryDisplacementMassTerm(ShallowWaterContinuityTerm):
r"""
Bathmetry mass displacement term, :math:`\partial \eta / \partial t + \partial \tilde{h} / \partial t`
The weak form reads
.. math::
\int_\Omega ( \partial \eta / \partial t + \partial \tilde{h} / \partial t ) \phi dx
= \int_\Omega (\partial \tilde{H} / \partial t) \phi dx
"""
def residual(self, solution):
if isinstance(solution, list):
uv, eta = solution
else:
uv, eta = split(solution)
f = inner(self.depth.wd_bathymetry_displacement(eta), self.eta_test)*self.dx
return -f
[docs]
class BaseShallowWaterEquation(Equation):
"""
Abstract base class for ShallowWaterEquations, ShallowWaterMomentumEquation
and FreeSurfaceEquation.
Provides common functionality to compute time steps and add either momentum
or continuity terms.
"""
def __init__(self, function_space,
depth, options):
super(BaseShallowWaterEquation, self).__init__(function_space)
self.depth = depth
self.options = options
[docs]
def add_momentum_terms(self, *args, tidal_farms=None):
self.add_term(ExternalPressureGradientTerm(*args), 'implicit')
self.add_term(HorizontalAdvectionTerm(*args), 'implicit')
self.add_term(HorizontalViscosityTerm(*args), 'explicit')
self.add_term(CoriolisTerm(*args), 'implicit')
self.add_term(WindStressTerm(*args), 'source')
self.add_term(AtmosphericPressureTerm(*args), 'source')
self.add_term(QuadraticDragTerm(*args), 'implicit')
self.add_term(LinearDragTerm(*args), 'implicit')
self.add_term(BoundaryDragTerm(*args), 'implicit')
self.add_term(BottomDrag3DTerm(*args), 'source')
self.add_term(MomentumSourceTerm(*args), 'source')
if tidal_farms:
self.add_term(TurbineDragTerm(*args, tidal_farms), 'implicit')
[docs]
def add_continuity_terms(self, *args):
self.add_term(HUDivTerm(*args), 'implicit')
self.add_term(ContinuitySourceTerm(*args), 'source')
[docs]
def residual_uv_eta(self, label, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions):
f = 0
for term in self.select_terms(label):
f += term.residual(uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions)
return f
[docs]
class ShallowWaterEquations(BaseShallowWaterEquation):
"""
2D depth-averaged shallow water equations in non-conservative form.
This defines the full 2D SWE equations :eq:`swe_freesurf` -
:eq:`swe_momentum`.
"""
def __init__(self, function_space, depth, options, tidal_farms=None):
"""
:arg function_space: Mixed function space where the solution belongs
:arg depth: :class: `DepthExpression` containing depth info
:arg options: :class:`.AttrDict` object containing all circulation model options
"""
super(ShallowWaterEquations, self).__init__(function_space, depth, options)
u_test, eta_test = TestFunctions(function_space)
u_space, eta_space = function_space.subfunctions
self.add_momentum_terms(u_test, u_space, eta_space, depth, options, tidal_farms=tidal_farms)
self.add_continuity_terms(eta_test, eta_space, u_space, depth, options)
self.bathymetry_displacement_mass_term = BathymetryDisplacementMassTerm(
eta_test, eta_space, u_space, depth, options)
[docs]
def mass_term(self, solution):
f = super(ShallowWaterEquations, self).mass_term(solution)
f += -self.bathymetry_displacement_mass_term.residual(solution)
return f
[docs]
def residual(self, label, solution, solution_old, fields, fields_old, bnd_conditions):
if isinstance(solution, list):
uv, eta = solution
else:
uv, eta = split(solution)
uv_old, eta_old = split(solution_old)
return self.residual_uv_eta(label, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions)
[docs]
class ModeSplit2DEquations(BaseShallowWaterEquation):
r"""
2D depth-averaged shallow water equations for mode splitting schemes.
Defines the equations :eq:`swe_freesurf_modesplit` -
:eq:`swe_momentum_modesplit`.
"""
def __init__(self, function_space, depth, options):
"""
:arg function_space: Mixed function space where the solution belongs
:arg depth: :class: `DepthExpression` containing depth info
:arg options: :class:`.AttrDict` object containing all circulation model options
"""
# TODO remove include_grad_* options as viscosity operator is omitted
super(ModeSplit2DEquations, self).__init__(function_space, depth, options)
u_test, eta_test = TestFunctions(function_space)
u_space, eta_space = function_space.subfunctions
self.add_momentum_terms(u_test, u_space, eta_space, depth, options)
self.add_continuity_terms(eta_test, eta_space, u_space, depth, options)
[docs]
def add_momentum_terms(self, *args):
self.add_term(ExternalPressureGradientTerm(*args), 'implicit')
self.add_term(CoriolisTerm(*args), 'explicit')
self.add_term(MomentumSourceTerm(*args), 'source')
self.add_term(AtmosphericPressureTerm(*args), 'source')
[docs]
def residual(self, label, solution, solution_old, fields, fields_old, bnd_conditions):
if isinstance(solution, list):
uv, eta = solution
else:
uv, eta = split(solution)
uv_old, eta_old = split(solution_old)
return self.residual_uv_eta(label, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions)
[docs]
class FreeSurfaceEquation(BaseShallowWaterEquation):
"""
2D free surface equation :eq:`swe_freesurf` in non-conservative form.
"""
def __init__(self, eta_test, eta_space, u_space, depth, options):
"""
:arg eta_test: test function of the elevation function space
:arg eta_space: elevation function space
:arg u_space: velocity function space
:arg function_space: Mixed function space where the solution belongs
:arg depth: :class: `DepthExpression` containing depth info
:arg options: :class:`.AttrDict` object containing all circulation model options
"""
super(FreeSurfaceEquation, self).__init__(eta_space, depth, options)
self.add_continuity_terms(eta_test, eta_space, u_space, depth, options)
self.bathymetry_displacement_mass_term = BathymetryDisplacementMassTerm(
eta_test, eta_space, u_space, depth, options)
[docs]
def mass_term(self, solution):
f = super().mass_term(solution)
f += -self.bathymetry_displacement_mass_term.residual([0, solution]) # expects [uv, eta] but uv is not used
return f
[docs]
def residual(self, label, solution, solution_old, fields, fields_old, bnd_conditions):
uv = fields['uv']
uv_old = fields_old['uv']
eta = solution
eta_old = solution_old
return self.residual_uv_eta(label, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions)
[docs]
class ShallowWaterMomentumEquation(BaseShallowWaterEquation):
"""
2D depth averaged momentum equation :eq:`swe_momentum` in non-conservative
form.
"""
def __init__(self, u_test, u_space, eta_space, depth, options, tidal_farms=None):
"""
:arg u_test: test function of the velocity function space
:arg u_space: velocity function space
:arg eta_space: elevation function space
:arg depth: :class: `DepthExpression` containing depth info
:arg options: :class:`.AttrDict` object containing all circulation model options
"""
super(ShallowWaterMomentumEquation, self).__init__(u_space, depth, options)
self.add_momentum_terms(u_test, u_space, eta_space, depth, options, tidal_farms=tidal_farms)
[docs]
def residual(self, label, solution, solution_old, fields, fields_old, bnd_conditions):
uv = solution
uv_old = solution_old
eta = fields['eta']
eta_old = fields_old['eta']
return self.residual_uv_eta(label, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions)
