# thetis package¶

## thetis.assembledschur module¶

class thetis.assembledschur.AssembledSchurPC[source]

Preconditioner for the Schur complement

The preconditioner matrix is assembled by explicitly matrix multiplying $$A10*Minv*A10$$.

Here:

• $$A01$$, $$A10$$ are the assembled sub-blocks of the saddle point
system. The form of this system needs to be supplied in the appctx argument to the solver as appctx['a'].
• $$Minv$$ is the inverse of the mass-matrix which is assembled as
assemble(v*u*dx, inverse=True), i.e. the element-wise inverse, where $$v$$ and $$u$$ are the test and trial of the $$A00$$ block. This gives the exact inverse of the mass matrix for a DG discretisation.

Create a PC context suitable for PETSc.

Matrix free preconditioners should inherit from this class and implement:

apply(pc, X, Y)[source]

Apply the preconditioner to X, putting the result in Y.

Both X and Y are PETSc Vecs, Y is not guaranteed to be zero on entry.

applyTranspose(pc, X, Y)[source]

Apply the transpose of the preconditioner to X, putting the result in Y.

Both X and Y are PETSc Vecs, Y is not guaranteed to be zero on entry.

initialize(pc)[source]

Initialize any state in this preconditioner.

update(pc)[source]

Update any state in this preconditioner.

view(pc, viewer=None)[source]

## thetis.callback module¶

Defines custom callback functions used to compute various metrics at runtime.

class thetis.callback.AccumulatorCallback(scalar_callback, solver_obj, **kwargs)[source]

Callback that evaluates a (scalar) functional involving integrals in both time and space. This callback can also be used to assemble time dependent objective functionals for adjoint simulations. Time integration is achieved using the trapezium rule.

Parameters: scalar_callback – Python function that returns a list of values of an objective functional. solver_obj – Thetis solver object **kwargs – any additional keyword arguments, see DiagnosticCallback
get_val()[source]
message_str(*args)[source]

A string representation.

Parameters: args – If provided, these will be the return value from __call__().
variable_names = ['spatial integral at current timestep']
class thetis.callback.CallbackManager[source]

Stores callbacks in different categories and provides methods for evaluating them.

Create callbacks and register them under 'export' mode

cb1 = VolumeConservation3DCallback(...)
cb2 = TracerMassConservationCallback(...)
cm = CallbackManager()


Evaluate callbacks, calls evaluate() method of all callbacks registered in the given mode.

cm.evaluate('export')

add(callback, mode)[source]

Add a callback under the given mode

Parameters: callback – a DiagnosticCallback object mode (str) – register callback under this mode
evaluate(mode, index=None)[source]

Evaluate all callbacks registered under the given mode

Parameters: mode (str) – evaluate all callbacks under this mode index (int) – if provided, sets the export index. Default behavior is to append to the file or stream.
class thetis.callback.DetectorsCallback(solver_obj, detector_locations, field_names, name, detector_names=None, **kwargs)[source]

Callback that evaluates the specified fields at the specified locations

Parameters: solver_obj – Thetis solver object detector_locations – List of x, y locations in which fields are to be interpolated. field_names – List of fields to be interpolated. name – Unique name for this callback and its associated set of detectors. This determines the name of the output h5 file (prefixed with diagnostic_). detector_names – List of names for each of the detectors. If not provided automatic names of the form detectorNNN are created where NNN is the (0-padded) detector number. kwargs – any additional keyword arguments, see DiagnosticCallback.
message_str(*args)[source]

A string representation.

Parameters: args – If provided, these will be the return value from __call__().
name

The name of the diagnostic

variable_names

Names of all scalar values

class thetis.callback.DiagnosticCallback(solver_obj, array_dim=1, attrs=None, outputdir=None, export_to_hdf5=True, append_to_log=True, include_time=True, hdf5_dtype='d')[source]

Bases: abc.ABC

A base class for all Callback classes

Parameters: solver_obj – Thetis solver object outputdir (str) – Custom directory where hdf5 files will be stored. By default solver’s output directory is used. array_dim (int) – Dimension of the output array. 1 for scalars. attrs (dict) – Additional attributes to be saved in the hdf5 file. export_to_hdf5 (bool) – If True, diagnostics will be stored in hdf5 format append_to_log (bool) – If True, callback output messages will be printed in log include_time (bool) – whether to include time in the hdf5 file hdf5_dtype – Precision to use in hdf5 output: d for double precision (default), and f for single precision
evaluate(index=None)[source]

Evaluates callback and pushes values to log and hdf file (if enabled)

message_str(*args)[source]

A string representation.

Parameters: args – If provided, these will be the return value from __call__().
name

The name of the diagnostic

push_to_hdf5(time, args, index=None)[source]

Append values to HDF5 file.

Parameters: time – time stamp of entry args – the return value from __call__().
push_to_log(time, args)[source]

Push callback status message to log

Parameters: time – time stamp of entry args – the return value from __call__().
set_write_mode(mode)[source]

Define whether to create a new hdf5 file or append to an existing one

Parameters: mode (str) – Either ‘create’ (default) or ‘append’
variable_names

Names of all scalar values

class thetis.callback.DiagnosticHDF5(filename, varnames, array_dim=1, attrs=None, comm=<mpi4py.MPI.Intracomm object>, new_file=True, dtype='d', include_time=True)[source]

Bases: object

A HDF5 file for storing diagnostic time series arrays.

Parameters: filename (str) – Full filename of the HDF5 file. varnames – List of variable names that the diagnostic callback provides array_dim (int) – Dimension of the output array. 1 for scalars. attrs (dict) – Additional attributes to be saved in the hdf5 file. comm – MPI communicator new_file (bool) – Define whether to create a new hdf5 file or append to an existing one (if any)
export(variables, time=None, index=None)[source]

Appends a new entry of (time, variables) to the file.

The HDF5 is updated immediately.

Parameters: time (float) – time stamp of entry variables (tuple of float) – values of entry index (int) – If provided, defines the time index in the file
class thetis.callback.MinMaxConservationCallback(minmax_callback, solver_obj, **kwargs)[source]

Base class for callbacks that check conservation of a minimum/maximum

Parameters: minmax_callback – Python function that takes the solver object as an argument and returns a (min, max) value tuple solver_obj – Thetis solver object kwargs – any additional keyword arguments, see DiagnosticCallback.
message_str(*args)[source]

A string representation.

Parameters: args – If provided, these will be the return value from __call__().
variable_names = ['min_value', 'max_value', 'undershoot', 'overshoot']
class thetis.callback.ScalarConservationCallback(scalar_callback, solver_obj, **kwargs)[source]

Base class for callbacks that check conservation of a scalar quantity

Creates scalar conservation check callback object

Parameters: scalar_callback – Python function that takes the solver object as an argument and returns a scalar quantity of interest solver_obj – Thetis solver object kwargs – any additional keyword arguments, see DiagnosticCallback.
message_str(*args)[source]

A string representation.

Parameters: args – If provided, these will be the return value from __call__().
variable_names = ['integral', 'relative_difference']
class thetis.callback.TracerMassConservation2DCallback(tracer_name, solver_obj, **kwargs)[source]

Checks conservation of depth-averaged tracer mass

Depth-averaged tracer mass is defined as the integral of 2D tracer multiplied by total depth.

Parameters: tracer_name – Name of the tracer. Use canonical field names as in FieldDict. solver_obj – Thetis solver object kwargs – any additional keyword arguments, see DiagnosticCallback.
name = 'tracer mass'
class thetis.callback.TracerMassConservationCallback(tracer_name, solver_obj, **kwargs)[source]

Checks conservation of total tracer mass

Parameters: tracer_name – Name of the tracer. Use canonical field names as in FieldDict. solver_obj – Thetis solver object kwargs – any additional keyword arguments, see DiagnosticCallback.
name = 'tracer mass'
class thetis.callback.TracerOvershootCallBack(tracer_name, solver_obj, **kwargs)[source]

Checks overshoots of the given tracer field.

Parameters: tracer_name – Name of the tracer. Use canonical field names as in FieldDict. solver_obj – Thetis solver object kwargs – any additional keyword arguments, see DiagnosticCallback.
name = 'tracer overshoot'
class thetis.callback.VolumeConservation2DCallback(solver_obj, **kwargs)[source]

Checks conservation of 2D volume (integral of water elevation field)

Parameters: solver_obj – Thetis solver object kwargs – any additional keyword arguments, see DiagnosticCallback.
name = 'volume2d'
class thetis.callback.VolumeConservation3DCallback(solver_obj, **kwargs)[source]

Checks conservation of 3D volume (volume of 3D mesh)

Parameters: solver_obj – Thetis solver object kwargs – any additional keyword arguments, see DiagnosticCallback.
name = 'volume3d'
thetis.callback.timed_region()

Log.Event(type cls, name, klass=None)

thetis.callback.timed_stage()

Log.Stage(type cls, name)

## thetis.configuration module¶

Utility function and extensions to traitlets used for specifying Thetis options

class thetis.configuration.ABCMetaHasTraits(name, bases, classdict)[source]

Bases: abc.ABCMeta, traitlets.traitlets.MetaHasTraits

Combined metaclass of ABCMeta and MetaHasTraits

Finish initializing the HasDescriptors class.

class thetis.configuration.BoundedFloat(default_value=traitlets.Undefined, bounds=None, **kwargs)[source]

Bases: traitlets.traitlets.Float

info()[source]
validate(obj, proposal)[source]
class thetis.configuration.BoundedInteger(default_value=traitlets.Undefined, bounds=None, **kwargs)[source]

Bases: traitlets.traitlets.Int

info()[source]
validate(obj, proposal)[source]
class thetis.configuration.FiredrakeCoefficient(default_value=traitlets.Undefined, allow_none=False, read_only=None, help=None, config=None, **kwargs)[source]

Bases: traitlets.traitlets.TraitType

Declare a traitlet.

If allow_none is True, None is a valid value in addition to any values that are normally valid. The default is up to the subclass. For most trait types, the default value for allow_none is False.

Extra metadata can be associated with the traitlet using the .tag() convenience method or by using the traitlet instance’s .metadata dictionary.

default_value = None
default_value_repr()[source]
info_text = 'a Firedrake Constant or Function'
validate(obj, value)[source]
class thetis.configuration.FiredrakeConstantTraitlet(default_value=traitlets.Undefined, allow_none=False, read_only=None, help=None, config=None, **kwargs)[source]

Bases: traitlets.traitlets.TraitType

Declare a traitlet.

If allow_none is True, None is a valid value in addition to any values that are normally valid. The default is up to the subclass. For most trait types, the default value for allow_none is False.

Extra metadata can be associated with the traitlet using the .tag() convenience method or by using the traitlet instance’s .metadata dictionary.

default_value = None
default_value_repr()[source]
info_text = 'a Firedrake Constant'
validate(obj, value)[source]
class thetis.configuration.FiredrakeScalarExpression(default_value=traitlets.Undefined, allow_none=False, read_only=None, help=None, config=None, **kwargs)[source]

Bases: traitlets.traitlets.TraitType

Declare a traitlet.

If allow_none is True, None is a valid value in addition to any values that are normally valid. The default is up to the subclass. For most trait types, the default value for allow_none is False.

Extra metadata can be associated with the traitlet using the .tag() convenience method or by using the traitlet instance’s .metadata dictionary.

default_value = None
default_value_repr()[source]
info_text = 'a scalar UFL expression'
validate(obj, value)[source]
class thetis.configuration.FiredrakeVectorExpression(default_value=traitlets.Undefined, allow_none=False, read_only=None, help=None, config=None, **kwargs)[source]

Bases: traitlets.traitlets.TraitType

Declare a traitlet.

If allow_none is True, None is a valid value in addition to any values that are normally valid. The default is up to the subclass. For most trait types, the default value for allow_none is False.

Extra metadata can be associated with the traitlet using the .tag() convenience method or by using the traitlet instance’s .metadata dictionary.

default_value = None
default_value_repr()[source]
info_text = 'a vector UFL expression'
validate(obj, value)[source]
class thetis.configuration.FrozenConfigurable(*args, **kwargs)[source]

Bases: thetis.configuration.OptionsBase, traitlets.config.configurable.Configurable

A Configurable class that only allows adding new attributes in the class definition or when self._isfrozen is False.

class thetis.configuration.FrozenHasTraits(*args, **kwargs)[source]

Bases: thetis.configuration.OptionsBase, traitlets.traitlets.HasTraits

class thetis.configuration.NonNegativeFloat(default_value=traitlets.Undefined, allow_none=False, **kwargs)[source]

Bases: traitlets.traitlets.Float

info()[source]
validate(obj, proposal)[source]
class thetis.configuration.NonNegativeInteger(default_value=traitlets.Undefined, allow_none=False, **kwargs)[source]

Bases: traitlets.traitlets.Int

info()[source]
validate(obj, proposal)[source]
class thetis.configuration.OptionsBase[source]

Bases: object

Abstract base class for all options classes

name

Human readable name of the configurable object

update(options)[source]

Assign options from another container

Parameters: options – Either a dictionary of options or another HasTraits object
class thetis.configuration.PETScSolverParameters(trait=None, traits=None, default_value=traitlets.Undefined, **kwargs)[source]

Bases: traitlets.traitlets.Dict

PETSc solver options dictionary

Create a dict trait type from a Python dict.

The default value is created by doing dict(default_value), which creates a copy of the default_value.

trait : TraitType [ optional ]
The specified trait type to check and use to restrict contents of the Container. If unspecified, trait types are not checked.
traits : Dictionary of trait types [ optional ]
A Python dictionary containing the types that are valid for restricting the content of the Dict Container for certain keys.
default_value : SequenceType [ optional ]
The default value for the Dict. Must be dict, tuple, or None, and will be cast to a dict if not None. If trait is specified, the default_value must conform to the constraints it specifies.
default_value = None
info_text = 'a PETSc solver options dictionary'
validate(obj, value)[source]
class thetis.configuration.PairedEnum(values, paired_name, default_value=traitlets.Undefined, **kwargs)[source]

Bases: traitlets.traitlets.Enum

A enum whose value must be in a given sequence.

This enum controls a slaved option, with default values provided here.

Parameters: values – iterable of (value, HasTraits class) pairs. The HasTraits class will be called (with no arguments) to create default values if necessary. paired_name – trait name this enum is paired with. default_value – default value.
info()[source]

Returns a description of the trait.

class thetis.configuration.PositiveFloat(default_value=traitlets.Undefined, allow_none=False, **kwargs)[source]

Bases: traitlets.traitlets.Float

info()[source]
validate(obj, proposal)[source]
class thetis.configuration.PositiveInteger(default_value=traitlets.Undefined, allow_none=False, **kwargs)[source]

Bases: traitlets.traitlets.Int

info()[source]
validate(obj, proposal)[source]
thetis.configuration.attach_paired_options(name, name_trait, value_trait)[source]

Attach paired options to a Configurable object.

Parameters: name – the name of the enum trait name_trait – the enum trait (a PairedEnum) value_trait – the slaved value trait.
thetis.configuration.rst_all_options(cls, nspace=0, prefix=None)[source]

Recursively generate rST for a provided Configurable class.

Parameters: cls – The Configurable class. nspace – Indentation level. prefix – Prefix to use for new traits.

## thetis.coordsys module¶

Generic methods for converting data between different spatial coordinate systems. Uses pyproj library.

class thetis.coordsys.VectorCoordSysRotation(source_sys, target_sys, x, y)[source]

Bases: object

Rotates vectors defined in source_sys coordinates to a different coordinate system.

Parameters: source_sys – pyproj coordinate system where (x, y) are defined in target_sys – target pyproj coordinate system x – x coordinate y – y coordinate
thetis.coordsys.convert_coords(source_sys, target_sys, x, y)[source]

Converts coordinates from source_sys to target_sys

This function extends pyproj.transform method by handling NaNs correctly.

Parameters: source_sys – pyproj coordinate system where (x, y) are defined in target_sys – target pyproj coordinate system x (float or numpy.array_like) – x coordinate y (float or numpy.array_like) – y coordinate
thetis.coordsys.get_vector_rotation_matrix(source_sys, target_sys, x, y, delta=None)[source]

Estimate rotation matrix that converts vectors defined in source_sys to target_sys.

Assume that we have a vector field defined in source_sys: vectors located at (x, y) define the x and y components. We can then rotate the vectors to represent x2 and y2 components of the target_sys by applying a local rotation:

R, theta = get_vector_rotation_matrix(source_sys, target_sys, x, lat)
v_xy = numpy.array([[v_x], [v_y]])
v_new = numpy.matmul(R, v_xy)
v_x2, v_y2 = v_new


## thetis.coupled_timeintegrator module¶

Time integrators for solving coupled 2D-3D system of equations.

class thetis.coupled_timeintegrator.CoupledLeapFrogAM3(solver)[source]

Leap-Frog Adams-Moulton 3 time integrator for coupled 2D-3D problem

This is an ALE time integrator. Implementation follows the SLIM time integrator by Karna et al (2013)

Karna, et al. (2013). A baroclinic discontinuous Galerkin finite element model for coastal flows. Ocean Modelling, 61(0):1-20. http://dx.doi.org/10.1016/j.ocemod.2012.09.009

advance(t, update_forcings=None, update_forcings3d=None)[source]

Advances the equations for one time step

Parameters: t (float) – simulation time update_forcings – Optional user-defined function that takes simulation time and updates time-dependent boundary conditions of the 2D equations. update_forcings3d – Optional user defined function that updates boundary conditions of the 3D equations
integrator_2d
integrator_3d
integrator_vert_3d
class thetis.coupled_timeintegrator.CoupledTimeIntegrator(solver)[source]

Base class of mode-split time integrators that use 2D, 3D and implicit 3D time integrators.

Parameters: solver – FlowSolver object
initialize()[source]

Initial conditions are read from fields dictionary.

integrator_2d

time integrator for 2D equations

integrator_3d

time integrator for explicit 3D equations (momentum, tracers)

integrator_vert_3d

time integrator for implicit 3D equations (vertical diffusion)

set_dt(dt, dt_2d)[source]

Set time step for the coupled time integrator

Parameters: dt (float) – Time step. This is the master (macro) time step used to march the 3D equations. dt_2d (float) – Time step for 2D equations. For consistency dt_2d must be an integer fraction of dt. If 2D solver is implicit set dt_2d equal to dt.
class thetis.coupled_timeintegrator.CoupledTimeIntegratorBase(solver)[source]

Base class for coupled 2D-3D time integrators

Provides common functionality for updating diagnostic fields etc.

Parameters: solver – FlowSolver object
class thetis.coupled_timeintegrator.CoupledTwoStageRK(solver)[source]

Coupled time integrator based on SSPRK(2,2) scheme

This ALE time integration method uses SSPRK(2,2) scheme to advance the 3D equations and a compatible implicit Trapezoid method to advance the 2D equations. Backward Euler scheme is used for vertical diffusion.

advance(t, update_forcings=None, update_forcings3d=None)[source]

Advances the equations for one time step

Parameters: t (float) – simulation time update_forcings – Optional user-defined function that takes simulation time and updates time-dependent boundary conditions of the 2D equations. update_forcings3d – Optional user defined function that updates boundary conditions of the 3D equations
compute_mesh_velocity(istage)[source]

Computes mesh velocity for stage i

Must be called after updating the 2D mode.

Parameters: istage (int) – stage of the Runge-Kutta iteration
integrator_2d
integrator_3d
integrator_vert_3d
store_elevation(istage)[source]

Store current elevation field for computing mesh velocity

Must be called before updating the 2D mode.

Parameters: istage (int) – stage of the Runge-Kutta iteration
thetis.coupled_timeintegrator.timed_region()

Log.Event(type cls, name, klass=None)

thetis.coupled_timeintegrator.timed_stage()

Log.Stage(type cls, name)

## thetis.coupled_timeintegrator_2d module¶

Time integrators for solving coupled shallow water equations with one tracer.

class thetis.coupled_timeintegrator_2d.CoupledCrankEuler2D(solver)[source]
Parameters: solver – FlowSolver object
swe_integrator
tracer_integrator
class thetis.coupled_timeintegrator_2d.CoupledCrankNicolson2D(solver)[source]
Parameters: solver – FlowSolver object
swe_integrator
tracer_integrator
class thetis.coupled_timeintegrator_2d.CoupledTimeIntegrator2D(solver)[source]

Base class of time integrator for coupled shallow water and tracer equations

Parameters: solver – FlowSolver object
advance(t, update_forcings=None)[source]

Advances equations for one time step

Parameters: t (float) – simulation time update_forcings – user-defined function that takes the simulation time and updates any time-dependent boundary conditions
initialize(solution2d)[source]

Initial conditions are read from fields dictionary.

set_dt(dt)[source]

Set time step for the coupled time integrator

Parameters: dt (float) – Time step.
swe_integrator

time integrator for the shallow water equations

tracer_integrator

time integrator for the tracer equation

thetis.coupled_timeintegrator_2d.timed_region()

Log.Event(type cls, name, klass=None)

thetis.coupled_timeintegrator_2d.timed_stage()

Log.Stage(type cls, name)

## thetis.equation module¶

Implements Equation and Term classes.

class thetis.equation.Equation(function_space)[source]

Bases: object

Implements an equation, made out of terms.

Parameters: function_space – the FunctionSpace the solution belongs to
SUPPORTED_LABELS = frozenset({'source', 'explicit', 'nonlinear', 'implicit'})

Valid labels for terms, indicating how they should be treated in the time integrator.

source
The term is a source term, i.e. does not depend on the solution.
explicit
The term should be treated explicitly
implicit
The term should be treated implicitly
nonlinear
The term is nonlinear and should be treated fully implicitly
add_term(term, label)[source]

Adds a term in the equation

Parameters: term – Term object to add_term label (string) – Assign a label to the term. Valid labels are given by SUPPORTED_LABELS.
jacobian(label, solution, solution_old, fields, fields_old, bnd_conditions)[source]

Returns an UFL form of the Jacobian by summing up all the Jacobians of the terms.

Sign convention: all terms are assumed to be on the right hand side of the equation: d(u)/dt = term.

Parameters: label – string defining the type of terms to sum up. Currently one of ‘source’|’explicit’|’implicit’|’nonlinear’. Can be a list of multiple labels, or ‘all’ in which case all defined terms are summed. solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
label_term(term, label)[source]

Assings a label to the given term(s).

Parameters: term – Term object, or a tuple of terms label – string label to assign
mass_term(solution)[source]

Returns default mass matrix term for the solution function space.

Returns: UFL form of the mass term
residual(label, solution, solution_old, fields, fields_old, bnd_conditions)[source]

Returns an UFL form of the residual by summing up all the terms with the desired label.

Sign convention: all terms are assumed to be on the right hand side of the equation: d(u)/dt = term.

Parameters: label – string defining the type of terms to sum up. Currently one of ‘source’|’explicit’|’implicit’|’nonlinear’. Can be a list of multiple labels, or ‘all’ in which case all defined terms are summed. solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
select_terms(label)[source]

Generator function that selects terms by label(s).

label can be a single label (e.g. ‘explicit’), ‘all’ or a tuple of labels.

class thetis.equation.Term(function_space)[source]

Bases: object

Implements a single term of an equation.

Note

Sign convention: all terms are assumed to be on the right hand side of the equation: d(u)/dt = term.

Parameters: function_space – the FunctionSpace the solution belongs to
jacobian(solution, solution_old, fields, fields_old, bnd_conditions)[source]

Returns an UFL form of the Jacobian of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
residual(solution, solution_old, fields, fields_old, bnd_conditions)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
thetis.equation.timed_region()

Log.Event(type cls, name, klass=None)

thetis.equation.timed_stage()

Log.Stage(type cls, name)

## thetis.exporter module¶

Routines for handling file exports.

class thetis.exporter.ExportManager(outputdir, fields_to_export, functions, field_metadata, export_type='vtk', next_export_ix=0, verbose=False)[source]

Bases: object

Helper object for exporting multiple fields simultaneously

from .field_defs import field_metadata
field_dict = {'elev_2d': Function(...), 'uv_3d': Function(...), ...}
e = exporter.ExportManager('mydirectory',
['elev_2d', 'uv_3d', salt_3d'],
field_dict,
export_type='vtk')
e.export()

Parameters: outputdir (string) – directory where files are stored fields_to_export (list of strings) – list of fields to export functions – dict that contains all existing Function s field_metadata – dict of all field metadata. See field_defs export_type (str) – export format, either ‘vtk’ or ‘hdf5’ next_export_ix (int) – index for next export (default 0) verbose (bool) – print debug info to stdout
export()[source]

Export all designated functions to disk

Increments export index by 1.

export_bathymetry(bathymetry_2d)[source]

Special function to export 2D bathymetry data to disk

Bathymetry does not vary in time so this only needs to be called once.

Parameters: bathymetry_2d – 2D bathymetry Function
set_next_export_ix(next_export_ix)[source]

Set export index to all child exporters

class thetis.exporter.ExporterBase(filename, outputdir, next_export_ix=0, verbose=False)[source]

Bases: object

Base class for exporter objects.

Parameters: filename (string) – output file name (without directory) outputdir (string) – directory where file is stored next_export_ix (int) – set the index for next output verbose (bool) – print debug info to stdout
export(function)[source]

Exports given function to disk

set_next_export_ix(next_export_ix)[source]

Sets the index of next export

class thetis.exporter.HDF5Exporter(function_space, outputdir, filename_prefix, next_export_ix=0, verbose=False)[source]

Stores fields in disk in native discretization using HDF5 containers

Create exporter object for given function.

Parameters: function_space (FunctionSpace) – space where the exported functions belong outputdir (string) – directory where outputs will be stored filename_prefix (string) – prefix of output filename. Filename is prefix_nnnnn.h5 where nnnnn is the export number. next_export_ix (int) – index for next export (default 0) verbose (bool) – print debug info to stdout
export(function)[source]

Export function to disk.

Increments export index by 1.

Parameters: function – Function to export
export_as_index(iexport, function)[source]

Export function to disk using the specified export index number

Parameters: iexport (int) – export index >= 0 function – Function to export
gen_filename(iexport)[source]

Generate file name ‘prefix_nnnnn.h5’ for i-th export

Parameters: iexport (int) – export index >= 0
load(iexport, function)[source]

Loads nodal values from disk and assigns to the given function

Parameters: iexport (int) – export index >= 0 function – target Function
class thetis.exporter.VTKExporter(fs_visu, func_name, outputdir, filename, next_export_ix=0, project_output=False, coords_dg=None, verbose=False)[source]

Class that handles Paraview VTK file exports

Parameters: fs_visu – function space where input function will be cast before exporting func_name – name of the function outputdir – output directory filename – name of the pvd file next_export_ix (int) – index for next export (default 0) project_output (bool) – project function to output space instead of interpolating coords_dg (bool) – Discontinuous coordinate field. Needed to avoid allocating new coordinate field in case of discontinuous export functions. verbose (bool) – print debug info to stdout
export(function)[source]

Exports given function to disk

set_next_export_ix(next_export_ix)[source]

Sets the index of next export

thetis.exporter.get_visu_space(fs)[source]

Returns an appropriate VTK visualization space for a function space

Parameters: fs – function space function space for VTK visualization
thetis.exporter.is_2d(fs)[source]

Tests wether a function space is 2D or 3D

thetis.exporter.timed_region()

Log.Event(type cls, name, klass=None)

thetis.exporter.timed_stage()

Log.Stage(type cls, name)

## thetis.field_defs module¶

Definitions and meta data of fields

thetis.field_defs.field_metadata = {'baroc_head_3d': {'filename': 'BaroHead3d', 'name': 'Baroclinic head', 'shortname': 'Baroclinic head', 'unit': 'm'}, 'bathymetry_2d': {'filename': 'bathymetry2d', 'name': 'Bathymetry', 'shortname': 'Bathymetry', 'unit': 'm'}, 'bathymetry_3d': {'filename': 'bathymetry3d', 'name': 'Bathymetry', 'shortname': 'Bathymetry', 'unit': 'm'}, 'bottom_drag_2d': {'filename': 'BottomDrag2d', 'name': 'Bottom drag coefficient', 'shortname': 'Bottom drag coefficient', 'unit': ''}, 'bottom_drag_3d': {'filename': 'BottomDrag3d', 'name': 'Bottom drag coefficient', 'shortname': 'Bottom drag coefficient', 'unit': ''}, 'buoy_freq_3d': {'filename': 'BuoyFreq3d', 'name': 'Buoyancy frequency squared', 'shortname': 'Buoyancy frequency squared', 'unit': 's-2'}, 'coriolis_2d': {'filename': 'coriolis_2d', 'name': 'Coriolis parameter', 'shortname': 'Coriolis parameter', 'unit': 's-1'}, 'coriolis_3d': {'filename': 'coriolis_3d', 'name': 'Coriolis parameter', 'shortname': 'Coriolis parameter', 'unit': 's-1'}, 'density_3d': {'filename': 'Density3d', 'name': 'Water density', 'shortname': 'Density', 'unit': 'kg m-3'}, 'eddy_diff_3d': {'filename': 'EddyDiff3d', 'name': 'Eddy diffusivity', 'shortname': 'Eddy diffusivity', 'unit': 'm2 s-1'}, 'eddy_visc_3d': {'filename': 'EddyVisc3d', 'name': 'Eddy Viscosity', 'shortname': 'Eddy Viscosity', 'unit': 'm2 s-1'}, 'elev_2d': {'filename': 'Elevation2d', 'name': 'Water elevation', 'shortname': 'Elevation', 'unit': 'm'}, 'elev_3d': {'filename': 'Elevation3d', 'name': 'Water elevation', 'shortname': 'Elevation', 'unit': 'm'}, 'elev_cg_2d': {'filename': 'ElevationCG2d', 'name': 'Water elevation CG', 'shortname': 'Elevation', 'unit': 'm'}, 'elev_cg_3d': {'filename': 'ElevationCG3d', 'name': 'Water elevation CG', 'shortname': 'Elevation', 'unit': 'm'}, 'eps_3d': {'filename': 'TurbEps3d', 'name': 'TKE dissipation rate', 'shortname': 'TKE dissipation rate', 'unit': 'm2 s-2'}, 'h_elem_size_2d': {'filename': 'h_elem_size_2d', 'name': 'Element size in horizontal direction', 'shortname': 'Horizontal element size', 'unit': 'm'}, 'h_elem_size_3d': {'filename': 'h_elem_size_3d', 'name': 'Element size in horizontal direction', 'shortname': 'Horizontal element size', 'unit': 'm'}, 'hcc_metric_3d': {'filename': 'HCCMetric3d', 'name': 'HCC mesh quality', 'shortname': 'HCC metric', 'unit': '-'}, 'int_pg_3d': {'filename': 'IntPG3d', 'name': 'Internal pressure gradient', 'shortname': 'Int. Pressure gradient', 'unit': 'm s-2'}, 'len_3d': {'filename': 'TurbLen3d', 'name': 'Turbulent length scale', 'shortname': 'Turbulent length scale', 'unit': 'm'}, 'max_h_diff': {'filename': 'MaxHDiffusivity3d', 'name': 'Maximum stable horizontal diffusivity', 'shortname': 'Maximum horizontal diffusivity', 'unit': 'm2 s-1'}, 'parab_visc_3d': {'filename': 'ParabVisc3d', 'name': 'Parabolic Viscosity', 'shortname': 'Parabolic Viscosity', 'unit': 'm2 s-1'}, 'psi_3d': {'filename': 'TurbPsi3d', 'name': 'Turbulence psi variable', 'shortname': 'Turbulence psi variable', 'unit': ''}, 'salt_3d': {'filename': 'Salinity3d', 'name': 'Water salinity', 'shortname': 'Salinity', 'unit': 'psu'}, 'shear_freq_3d': {'filename': 'ShearFreq3d', 'name': 'Vertical shear frequency squared', 'shortname': 'Vertical shear frequency squared', 'unit': 's-2'}, 'smag_visc_3d': {'filename': 'SmagViscosity3d', 'name': 'Smagorinsky viscosity', 'shortname': 'Smagorinsky viscosity', 'unit': 'm2 s-1'}, 'split_residual_2d': {'filename': 'SplitResidual2d', 'name': 'Momentum eq. residual for mode splitting', 'shortname': 'Momentum residual', 'unit': 'm s-2'}, 'temp_3d': {'filename': 'Temperature3d', 'name': 'Water temperature', 'shortname': 'Temperature', 'unit': 'C'}, 'tke_3d': {'filename': 'TurbKEnergy3d', 'name': 'Turbulent Kinetic Energy', 'shortname': 'Turbulent Kinetic Energy', 'unit': 'm2 s-2'}, 'tracer_2d': {'filename': 'Tracer2d', 'name': 'Depth averaged tracer', 'shortname': 'Tracer', 'unit': ''}, 'uv_2d': {'filename': 'Velocity2d', 'name': 'Depth averaged velocity', 'shortname': 'Depth averaged velocity', 'unit': 'm s-1'}, 'uv_3d': {'filename': 'Velocity3d', 'name': 'Horizontal velocity', 'shortname': 'Horizontal velocity', 'unit': 'm s-1'}, 'uv_bottom_2d': {'filename': 'BottomVelo2d', 'name': 'Bottom velocity', 'shortname': 'Bottom velocity', 'unit': 'm s-1'}, 'uv_bottom_3d': {'filename': 'BottomVelo3d', 'name': 'Bottom velocity', 'shortname': 'Bottom velocity', 'unit': 'm s-1'}, 'uv_dav_2d': {'filename': 'DAVelocity2d', 'name': 'Depth averaged velocity', 'shortname': 'Depth averaged velocity', 'unit': 'm s-1'}, 'uv_dav_3d': {'filename': 'DAVelocity3d', 'name': 'Depth averaged velocity', 'shortname': 'Depth averaged velocity', 'unit': 'm s-1'}, 'uv_mag_3d': {'filename': 'VeloMag3d', 'name': 'Magnitude of horizontal velocity', 'shortname': 'Velocity magnitude', 'unit': 'm s-1'}, 'uv_p1_3d': {'filename': 'VeloCG3d', 'name': 'P1 projection of horizontal velocity', 'shortname': 'P1 Velocity', 'unit': 'm s-1'}, 'v_elem_size_2d': {'filename': 'VElemSize2d', 'name': 'Element size in vertical direction', 'shortname': 'Vertical element size', 'unit': 'm'}, 'v_elem_size_3d': {'filename': 'VElemSize3d', 'name': 'Element size in vertical direction', 'shortname': 'Vertical element size', 'unit': 'm'}, 'w_3d': {'filename': 'VertVelo3d', 'name': 'Vertical velocity', 'shortname': 'Vertical velocity', 'unit': 'm s-1'}, 'w_mesh_3d': {'filename': 'MeshVelo3d', 'name': 'Mesh velocity', 'shortname': 'Mesh velocity', 'unit': 'm s-1'}, 'w_mesh_surf_2d': {'filename': 'SurfMeshVelo3d', 'name': 'Surface mesh velocity', 'shortname': 'Surface mesh velocity', 'unit': 'm s-1'}, 'w_mesh_surf_3d': {'filename': 'SurfMeshVelo3d', 'name': 'Surface mesh velocity', 'shortname': 'Surface mesh velocity', 'unit': 'm s-1'}, 'wind_stress_3d': {'filename': 'wind_stress_3d', 'name': 'Wind stress', 'shortname': 'Wind stress', 'unit': 'Pa'}, 'z_bottom_2d': {'filename': 'ZBottom2d', 'name': 'Bottom cell z coordinates', 'shortname': 'Bottom cell z coordinates', 'unit': 'm'}, 'z_coord_3d': {'filename': 'ZCoord3d', 'name': 'Mesh z coordinates', 'shortname': 'Z coordinates', 'unit': 'm'}, 'z_coord_ref_3d': {'filename': 'ZCoordRef3d', 'name': 'Static mesh z coordinates', 'shortname': 'Z coordinates', 'unit': 'm'}}

Dictionary that contains the meta data of each field.

Required meta data entries are:

• shortname: description used in visualization etc
• unit: SI unit of the field
• filename: filename for output files

The naming convention for field keys is snake_case: field_name_3d

## thetis.implicitexplicit module¶

Implicit-explicit time integrators

class thetis.implicitexplicit.IMEXEuler(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, solver_parameters_dirk={})[source]

Forward-Backward Euler

dirk_class
erk_class
class thetis.implicitexplicit.IMEXGeneric(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, solver_parameters_dirk={})[source]

Generic implementation of Runge-Kutta Implicit-Explicit schemes

Derived classes must define the implicit dirk_class and explicit erk_class Runge-Kutta time integrator classes.

This method solves the linearized equations: All implicit terms are fed to the implicit solver, while all the other terms are fed to the explicit solver. In case of non-linear terms proper linearization must defined in the equation using the two solution functions (solution, solution_old)

advance(t, update_forcings=None)[source]

Advances equations for one time step.

dirk_class

Implicit DIRK class

erk_class

Explicit Runge-Kutta class

get_final_solution()[source]

Evaluates the final solution.

initialize(solution)[source]

set_dt(dt)[source]

Update time step

Parameters: dt (float) – time step
solve_stage(i_stage, t, update_forcings=None)[source]

Solves i-th stage

update_solver()[source]

Create solver objects

class thetis.implicitexplicit.IMEXLPUM2(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, solver_parameters_dirk={})[source]

SSP-IMEX RK scheme (20) in Higureras et al. (2014)

CFL coefficient is 2.0

Higueras et al (2014). Optimized strong stability preserving IMEX Runge-Kutta methods. Journal of Computational and Applied Mathematics 272(2014) 116-140. http://dx.doi.org/10.1016/j.cam.2014.05.011

dirk_class
erk_class
class thetis.implicitexplicit.IMEXLSPUM2(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, solver_parameters_dirk={})[source]

SSP-IMEX RK scheme (17) in Higureras et al. (2014)

CFL coefficient is 2.0

Higueras et al (2014). Optimized strong stability preserving IMEX Runge-Kutta methods. Journal of Computational and Applied Mathematics 272(2014) 116-140. http://dx.doi.org/10.1016/j.cam.2014.05.011

dirk_class
erk_class
class thetis.implicitexplicit.IMEXMidpoint(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, solver_parameters_dirk={})[source]

Implicit-explicit midpoint scheme (1, 2, 2) from Ascher et al. (1997)

Ascher et al. (1997). Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Applied Numerical Mathematics, 25:151-167. http://dx.doi.org/10.1137/0732037

dirk_class
erk_class
thetis.implicitexplicit.timed_region()

Log.Event(type cls, name, klass=None)

thetis.implicitexplicit.timed_stage()

Log.Stage(type cls, name)

## thetis.interpolation module¶

Methods for interpolating data from structured data sets on Thetis fields.

Simple example of an atmospheric pressure interpolator:

def to_latlon(x, y, positive_lon=False):
# Converts mesh (x,y) points to coordinates used in the atm data
lon, lat = coordsys_spcs.spcs2lonlat(x, y)
if positive_lon and lon < 0.0:
lon += 360.
return lat, lon

class WRFInterpolator(object):
# Interpolates WRF atmospheric model data on 2D fields
def __init__(self, function_space, atm_pressure_field, ncfile_pattern, init_date):
self.atm_pressure_field = atm_pressure_field

# object that interpolates forcing data from structured grid on the local mesh
self.grid_interpolator = NetCDFLatLonInterpolator2d(function_space, to_latlon)
# reader object that can read fields from netCDF files, applies spatial interpolation
# object that can find previous/next time stamps in a collection of netCDF files
self.timesearch_obj = NetCDFTimeSearch(ncfile_pattern, init_date, NetCDFTimeParser)
# finally a linear intepolator class that performs linar interpolation in time

def set_fields(self, time):
# Evaluates forcing fields at the given time
pressure = self.interpolator(time)
self.atm_pressure_field.dat.data_with_halos[:] = pressure


Usage:

atm_pressure_2d = Function(solver_obj.function_spaces.P1_2d, name='atm pressure')
wrf_pattern = 'forcings/atm/wrf/wrf_air.2016_*_*.nc'
wrf_atm = WRFInterpolator(
solver_obj.function_spaces.P1_2d,
wind_stress_2d, atm_pressure_2d, wrf_pattern, init_date)
simulation_time = 3600.
wrf_atm.set_fields(simulation_time)

class thetis.interpolation.FileTreeReader[source]

Bases: object

Abstract base class of file tree reader object

class thetis.interpolation.GridInterpolator(grid_xyz, target_xyz, fill_mode=None, fill_value=nan, normalize=False)[source]

Bases: object

A reuseable griddata interpolator object.

Usage:

interpolator = GridInterpolator(source_xyz, target_xyz)
vals = interpolator(source_data)


Example:

x0 = np.linspace(0, 10, 10)
y0 = np.linspace(5, 10, 10)
X, Y = np.meshgrid(x, y)
x = X.ravel(); y = Y.ravel()
data = x + 25.*y
x_target = np.linspace(1, 10, 20)
y_target = np.linspace(5, 10, 20)
interpolator = GridInterpolator(np.vstack((x, y)).T, np.vstack((target_x, target_y)).T)
vals = interpolator(data)

Parameters: grid_xyz – Array of source grid coordinates, shape (npoints, 2) or (npoints, 3) target_xyz – Array of target grid coordinates, shape (n, 2) or (n, 3)
class thetis.interpolation.LinearTimeInterpolator(timesearch_obj, reader)[source]

Bases: object

Interpolates time series in time

User must provide timesearch_obj that finds time stamps from a file tree, and a reader that can read those time stamps into numpy arrays.

Previous/next data sets are cached in memory to avoid hitting disk every time.

class thetis.interpolation.NetCDFLatLonInterpolator2d(function_space, to_latlon)[source]

Interpolates netCDF data on a local 2D unstructured mesh

The intepolator is constructed for a single netCDF file that defines the source grid. Once the interpolator has been constructed, data can be read from any file that uses the same grid.

This routine returns the data in numpy arrays.

Usage:

fs = FunctionSpace(...)
myfunc = Function(fs, ...)
ncinterp2d = NetCDFLatLonInterpolator2d(fs, to_latlon, nc_filename)
val1, val2 = ncinterp2d.interpolate(nc_filename, ['var1', 'var2'], 10)
myfunc.dat.data_with_halos[:] = val1 + val2

Parameters: function_space – target Firedrake FunctionSpace to_latlon – Python function that converts local mesh coordinates to latitude and longitude: ‘lat, lon = to_latlon(x, y)’
interpolate(nc_filename, variable_list, itime)[source]

Interpolates data from a netCDF file onto Firedrake function space.

Parameters: nc_filename (str) – netCDF file to read variable_list – list of netCDF variable names to read itime (int) – time index to read list of numpy.arrays corresponding to variable_list
class thetis.interpolation.NetCDFSpatialInterpolator(grid_interpolator, variable_list)[source]

Wrapper class that provides FileTreeReader API for grid interpolators

class thetis.interpolation.NetCDFTimeParser(filename, time_variable_name='time', allow_gaps=False)[source]

Describes the time stamps stored in a netCDF file.

Construct a new object by scraping data from the given netcdf file.

Parameters: filename (str) – name of the netCDF file to read time_variable_name (str) – name of the time variable in the netCDF file (default: ‘time’) allow_gaps (bool) – if False, an error is raised if time step is not constant.
find_time_stamp(t, previous=False)[source]

Given time t, returns index of the next (previous) time stamp

raises IndexError if t is out of range, i.e. t > self.get_end_time() or t < self.get_start_time()

get_end_time()[source]

Returns the last time stamp in the file/data set

get_start_time()[source]

Returns the first time stamp in the file/data set

scalars = {'days': 86400.0, 'seconds': 1.0}
class thetis.interpolation.NetCDFTimeSearch(file_pattern, init_date, netcdf_class, *args, **kwargs)[source]

Finds a nearest time stamp in a collection of netCDF files.

find(simulation_time, previous=False)[source]

Find file that contains the given simulation time

Parameters: simulation_time (float) – simulation time in seconds previous (bool) – if True finds previous existing time stamp instead of next (default False). (filename, time index, simulation time) of found data
simulation_time_to_datetime(t)[source]
class thetis.interpolation.NetCDFTimeSeriesInterpolator(ncfile_pattern, variable_list, init_date, time_variable_name='time', scalars=None, allow_gaps=False)[source]

Bases: object

Reads and interpolates scalar time series from a sequence of netCDF files.

Parameters: ncfile_pattern (str) – file search pattern, e.g. “mydir/foo_*.nc” variable_list – list if netCDF variable names to read init_date (datetime.datetime) – simulation start time scalars – (optional) list of scalars; scale output variables by a factor.

Note

All the variables must have the same dimensions in the netCDF files. If the shapes differ, create separate interpolator instances.

class thetis.interpolation.NetCDFTimeSeriesReader(variable_list, time_variable_name='time')[source]

A simple netCDF reader that returns a time slice of the given variable.

This class does not interpolate the data in any way. Useful for interpolating time series.

class thetis.interpolation.SpatialInterpolator(function_space, to_latlon)[source]

Bases: object

Abstract base class for spatial interpolators that read data from disk

Parameters: function_space – target Firedrake FunctionSpace to_latlon – Python function that converts local mesh coordinates to latitude and longitude: ‘lat, lon = to_latlon(x, y)’
interpolate(filename, variable_list, itime)[source]

Interpolates data from the given file at given time step

class thetis.interpolation.SpatialInterpolator2d(function_space, to_latlon)[source]

Abstract spatial interpolator class that can interpolate onto a 2D Function

Parameters: function_space – target Firedrake FunctionSpace to_latlon – Python function that converts local mesh coordinates to latitude and longitude: ‘lat, lon = to_latlon(x, y)’
interpolate(filename, variable_list, time)[source]

Calls the interpolator object

class thetis.interpolation.TimeParser[source]

Bases: object

Abstract base class for time definition objects.

Defines the time span that a file (or data set) covers and provides a time index search routine.

find_time_stamp(t, previous=False)[source]

Given time t, returns index of the next (previous) time stamp

raises IndexError if t is out of range, i.e. t > self.get_end_time() or t < self.get_start_time()

get_end_time()[source]

Returns the last time stamp in the file/data set

get_start_time()[source]

Returns the first time stamp in the file/data set

class thetis.interpolation.TimeSearch[source]

Bases: object

Base class for searching nearest time steps in a file tree or database

find(time, previous=False)[source]

Find a next (previous) time stamp from a given time

Parameters: time (float) – input time stamp previous (bool) – if True, look for last time stamp before requested time. Otherwise returns next time stamp. a (filename, time_index, time) tuple

## thetis.limiter module¶

Slope limiters for discontinuous fields

class thetis.limiter.VertexBasedP1DGLimiter(p1dg_space, time_dependent_mesh=True)[source]

Vertex based limiter for P1DG tracer fields, see Kuzmin (2010)

Note

Currently only scalar fields are supported

Kuzmin (2010). A vertex-based hierarchical slope limiter for p-adaptive discontinuous Galerkin methods. Journal of Computational and Applied Mathematics, 233(12):3077-3085. http://dx.doi.org/10.1016/j.cam.2009.05.028

Parameters: p1dg_space – P1DG function space
apply(field)[source]

Applies the limiter on the given field (in place)

Parameters: field – Function to limit
compute_bounds(field)[source]

Re-compute min/max values of all neighbouring centroids

Parameters: field – Function to limit
thetis.limiter.assert_function_space(fs, family, degree)[source]

Checks the family and degree of function space.

Raises AssertionError if function space differs. If the function space lies on an extruded mesh, checks both spaces of the outer product.

Parameters: fs – function space family (string) – name of element family degree (int) – polynomial degree of the function space
thetis.limiter.timed_region()

Log.Event(type cls, name, klass=None)

thetis.limiter.timed_stage()

Log.Stage(type cls, name)

## thetis.log module¶

Loggers for Thetis

Creates two logger instances, one for general model output and one for debug, warning, error etc. messages.

To print to the model output stream, use print_output().

Debug, warning etc. messages are issued with debug(), info(), warning(), error(), critical() methods.

## thetis.momentum_eq module¶

3D momentum equation for hydrostatic Boussinesq flow.

The three dimensional momentum equation reads

(1)$\frac{\partial \textbf{u}}{\partial t} + \nabla_h \cdot (\textbf{u} \textbf{u}) + \frac{\partial \left(w\textbf{u} \right)}{\partial z} + f\textbf{e}_z\wedge\textbf{u} + g\nabla_h \eta + g\nabla_h r = \nabla_h \cdot \left( \nu_h \nabla_h \textbf{u} \right) + \frac{\partial }{\partial z}\left( \nu \frac{\partial \textbf{u}}{\partial z}\right)$

where $$\textbf{u}$$ and $$w$$ denote horizontal and vertical velocity, $$\nabla_h$$ is the horizontal gradient, $$\wedge$$ denotes the cross product, $$g$$ is the gravitational acceleration, $$f$$ is the Coriolis frequency, $$\textbf{e}_z$$ is a vertical unit vector, and $$\nu_h, \nu$$ stand for horizontal and vertical viscosity. Water density is given by $$\rho = \rho'(T, S, p) + \rho_0$$, where $$\rho_0$$ is a constant reference density. Above $$r$$ denotes the baroclinic head

(2)$r = \frac{1}{\rho_0} \int_{z}^\eta \rho' d\zeta.$

The internal pressure gradient is computed as a separate diagnostic field:

(3)$\mathbf{F}_{pg} = g\nabla_h r.$

In the case of purely barotropic problems the $$r$$ and $$\mathbf{F}_{pg}$$ fields are omitted.

When using mode splitting we split the velocity field into a depth average and a deviation, $$\textbf{u} = \bar{\textbf{u}} + \textbf{u}'$$. Following Higdon and de Szoeke (1997) we write an equation for the deviation $$\textbf{u}'$$:

(4)$\frac{\partial \textbf{u}'}{\partial t} = + \nabla_h \cdot (\textbf{u} \textbf{u}) + \frac{\partial \left(w\textbf{u} \right)}{\partial z} + f\textbf{e}_z\wedge\textbf{u}' + g\nabla_h r = \nabla_h \cdot \left( \nu_h \nabla_h \textbf{u} \right) + \frac{\partial }{\partial z}\left( \nu \frac{\partial \textbf{u}}{\partial z}\right)$

In (4) the external pressure gradient $$g\nabla_h \eta$$ vanishes and the Coriolis term only contains the deviation $$\textbf{u}'$$. Advection and diffusion terms are not changed.

Higdon and de Szoeke (1997). Barotropic-Baroclinic Time Splitting for Ocean Circulation Modeling. Journal of Computational Physics, 135(1):30-53. http://dx.doi.org/10.1006/jcph.1997.5733

class thetis.momentum_eq.MomentumEquation(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_nonlinear_equations=True, use_lax_friedrichs=True, use_bottom_friction=False)[source]

Hydrostatic 3D momentum equation (4) for mode split models

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_nonlinear_equations (bool) – If False defines the linear shallow water equations use_bottom_friction (bool) – If True includes bottom friction term
class thetis.momentum_eq.MomentumTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_nonlinear_equations=True, use_lax_friedrichs=True, use_bottom_friction=False)[source]

Generic term for momentum equation that provides commonly used members and mapping for boundary functions.

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_nonlinear_equations (bool) – If False defines the linear shallow water equations use_bottom_friction (bool) – If True includes bottom friction term
class thetis.momentum_eq.HorizontalAdvectionTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_nonlinear_equations=True, use_lax_friedrichs=True, use_bottom_friction=False)[source]

Horizontal advection term, $$\nabla_h \cdot (\textbf{u} \textbf{u})$$

$\int_\Omega \nabla_h \cdot (\textbf{u} \textbf{u}) \cdot \boldsymbol{\psi} dx = - \int_\Omega \nabla_h \boldsymbol{\psi} : (\textbf{u} \textbf{u}) dx + \int_{\mathcal{I}_h \cup \mathcal{I}_v} \textbf{u}^{\text{up}} \cdot \text{jump}(\boldsymbol{\psi} \textbf{n}_h) \cdot \text{avg}(\textbf{u}) dS$

where the right hand side has been integrated by parts; $$:$$ stand for the Frobenius inner product, $$\textbf{n}_h$$ is the horizontal projection of the normal vector, $$\textbf{u}^{\text{up}}$$ is the upwind value, and $$\text{jump}$$ and $$\text{avg}$$ denote the jump and average operators across the interface.

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_nonlinear_equations (bool) – If False defines the linear shallow water equations use_bottom_friction (bool) – If True includes bottom friction term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.momentum_eq.VerticalAdvectionTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_nonlinear_equations=True, use_lax_friedrichs=True, use_bottom_friction=False)[source]

Vertical advection term, $$\partial \left(w\textbf{u} \right)/(\partial z)$$

$\int_\Omega \frac{\partial \left(w\textbf{u} \right)}{\partial z} \cdot \boldsymbol{\psi} dx = - \int_\Omega \left( w \textbf{u} \right) \cdot \frac{\partial \boldsymbol{\psi}}{\partial z} dx + \int_{\mathcal{I}_{h}} \textbf{u}^{\text{up}} \cdot \text{jump}(\boldsymbol{\psi} n_z) \text{avg}(w) dS$
Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_nonlinear_equations (bool) – If False defines the linear shallow water equations use_bottom_friction (bool) – If True includes bottom friction term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.momentum_eq.HorizontalViscosityTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_nonlinear_equations=True, use_lax_friedrichs=True, use_bottom_friction=False)[source]

Horizontal viscosity term, $$- \nabla_h \cdot \left( \nu_h \nabla_h \textbf{u} \right)$$

Using the symmetric interior penalty method the weak form becomes

$\begin{split}\int_\Omega \nabla_h \cdot \left( \nu_h \nabla_h \textbf{u} \right) \cdot \boldsymbol{\psi} dx =& -\int_\Omega \nu_h (\nabla_h \boldsymbol{\psi}) : (\nabla_h \textbf{u})^T dx \\ &+ \int_{\mathcal{I}_h \cup \mathcal{I}_v} \text{jump}(\boldsymbol{\psi} \textbf{n}_h) \cdot \text{avg}( \nu_h \nabla_h \textbf{u}) dS + \int_{\mathcal{I}_h \cup \mathcal{I}_v} \text{jump}(\textbf{u} \textbf{n}_h) \cdot \text{avg}( \nu_h \nabla_h \boldsymbol{\psi}) dS \\ &- \int_{\mathcal{I}_h \cup \mathcal{I}_v} \sigma \text{avg}(\nu_h) \text{jump}(\textbf{u} \textbf{n}_h) \cdot \text{jump}(\boldsymbol{\psi} \textbf{n}_h) dS\end{split}$

where $$\sigma$$ is a penalty parameter, see Epshteyn and Riviere (2007).

Epshteyn and Riviere (2007). Estimation of penalty parameters for symmetric interior penalty Galerkin methods. Journal of Computational and Applied Mathematics, 206(2):843-872. http://dx.doi.org/10.1016/j.cam.2006.08.029

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_nonlinear_equations (bool) – If False defines the linear shallow water equations use_bottom_friction (bool) – If True includes bottom friction term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.momentum_eq.VerticalViscosityTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_nonlinear_equations=True, use_lax_friedrichs=True, use_bottom_friction=False)[source]

Vertical viscosity term, $$- \frac{\partial }{\partial z}\left( \nu \frac{\partial \textbf{u}}{\partial z}\right)$$

Using the symmetric interior penalty method the weak form becomes

$\begin{split}\int_\Omega \frac{\partial }{\partial z}\left( \nu \frac{\partial \textbf{u}}{\partial z}\right) \cdot \boldsymbol{\psi} dx =& -\int_\Omega \nu \frac{\partial \boldsymbol{\psi}}{\partial z} \cdot \frac{\partial \textbf{u}}{\partial z} dx \\ &+ \int_{\mathcal{I}_h} \text{jump}(\boldsymbol{\psi} n_z) \cdot \text{avg}\left(\nu \frac{\partial \textbf{u}}{\partial z}\right) dS + \int_{\mathcal{I}_h} \text{jump}(\textbf{u} n_z) \cdot \text{avg}\left(\nu \frac{\partial \boldsymbol{\psi}}{\partial z}\right) dS \\ &- \int_{\mathcal{I}_h} \sigma \text{avg}(\nu) \text{jump}(\textbf{u} n_z) \cdot \text{jump}(\boldsymbol{\psi} n_z) dS\end{split}$

where $$\sigma$$ is a penalty parameter, see Epshteyn and Riviere (2007).

Epshteyn and Riviere (2007). Estimation of penalty parameters for symmetric interior penalty Galerkin methods. Journal of Computational and Applied Mathematics, 206(2):843-872. http://dx.doi.org/10.1016/j.cam.2006.08.029

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_nonlinear_equations (bool) – If False defines the linear shallow water equations use_bottom_friction (bool) – If True includes bottom friction term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.momentum_eq.PressureGradientTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_nonlinear_equations=True, use_lax_friedrichs=True, use_bottom_friction=False)[source]

Internal pressure gradient term, $$g\nabla_h r$$

where $$r$$ is the baroclinic head (2). Let $$s$$ denote $$r/H$$. We can then write

$F_{IPG} = g\nabla_h((s -\bar{s}) H) + g\nabla_h\left(\frac{1}{H}\right) H^2\bar{s} + g s_{bot}\nabla_h h$

where $$\bar{s},s_{bot}$$ are the depth average and bottom value of $$s$$.

If $$s$$ belongs to a discontinuous function space, the first term is integrated by parts. Its weak form reads

$\int_\Omega g\nabla_h((s -\bar{s}) H) \cdot \boldsymbol{\psi} dx = - \int_\Omega g (s -\bar{s}) H \nabla_h \cdot \boldsymbol{\psi} dx + \int_{\mathcal{I}_h \cup \mathcal{I}_v} g (s -\bar{s}) H \boldsymbol{\psi} \cdot \textbf{n}_h dx$
Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_nonlinear_equations (bool) – If False defines the linear shallow water equations use_bottom_friction (bool) – If True includes bottom friction term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.momentum_eq.CoriolisTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_nonlinear_equations=True, use_lax_friedrichs=True, use_bottom_friction=False)[source]

Coriolis term, $$f\textbf{e}_z\wedge \bar{\textbf{u}}$$

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_nonlinear_equations (bool) – If False defines the linear shallow water equations use_bottom_friction (bool) – If True includes bottom friction term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.momentum_eq.BottomFrictionTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_nonlinear_equations=True, use_lax_friedrichs=True, use_bottom_friction=False)[source]

Quadratic bottom friction term, $$\tau_b = C_D \| \textbf{u}_b \| \textbf{u}_b$$

$\int_{\Gamma_{bot}} \tau_b \cdot \boldsymbol{\psi} dx = \int_{\Gamma_{bot}} C_D \| \textbf{u}_b \| \textbf{u}_b \cdot \boldsymbol{\psi} dx$

where $$\textbf{u}_b$$ is reconstructed velocity in the middle of the bottom element:

$\textbf{u}_b = \textbf{u}\Big|_{\Gamma_{bot}} + \frac{\partial \textbf{u}}{\partial z}\Big|_{\Gamma_{bot}} h_b,$

$$h_b$$ being half of the element height. For implicit solvers we linearize the stress as $$\tau_b = C_D \| \textbf{u}_b^{n} \| \textbf{u}_b^{n+1}$$

The drag is computed from the law-of-the wall

$C_D = \left( \frac{\kappa}{\ln (h_b + z_0)/z_0} \right)^2$

where $$z_0$$ is the bottom roughness length, read from z0_friction field. The user can override the $$C_D$$ value by providing quadratic_drag_coefficient field.

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_nonlinear_equations (bool) – If False defines the linear shallow water equations use_bottom_friction (bool) – If True includes bottom friction term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.momentum_eq.LinearDragTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_nonlinear_equations=True, use_lax_friedrichs=True, use_bottom_friction=False)[source]

Linear drag term, $$\tau_b = D \textbf{u}_b$$

where $$D$$ is the drag coefficient, read from linear_drag_coefficient field.

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_nonlinear_equations (bool) – If False defines the linear shallow water equations use_bottom_friction (bool) – If True includes bottom friction term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.momentum_eq.SourceTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_nonlinear_equations=True, use_lax_friedrichs=True, use_bottom_friction=False)[source]

Generic momentum source term

$F_s = \int_\Omega \sigma \cdot \boldsymbol{\psi} dx$

where $$\sigma$$ is a user defined vector valued Function.

This term also implements the wind stress, $$-\tau_w/(H \rho_0)$$. $$\tau_w$$ is a user-defined wind stress Function wind_stress. The weak form is

$F_w = \int_{\Gamma_s} \frac{1}{\rho_0} \tau_w \cdot \boldsymbol{\psi} dx$

Wind stress is only included if vertical viscosity is provided.

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_nonlinear_equations (bool) – If False defines the linear shallow water equations use_bottom_friction (bool) – If True includes bottom friction term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.momentum_eq.InternalPressureGradientCalculator(fields, options, bnd_functions, solver_parameters=None)[source]

Computes the internal pressure gradient term, $$g\nabla_h r$$

where $$r$$ is the baroclinic head (2).

If $$r$$ belongs to a discontinuous function space, the term is integrated by parts:

$\int_\Omega g \nabla_h r \cdot \boldsymbol{\psi} dx = - \int_\Omega g r \nabla_h \cdot \boldsymbol{\psi} dx + \int_{\mathcal{I}_h \cup \mathcal{I}_v} g \text{avg}(r) \text{jump}(\boldsymbol{\psi} \cdot \textbf{n}_h) dx$

Note

Due to the Term sign convention this term is assembled on the right-hand-side.

Parameters: solver – classFlowSolver object solver_parameters (dict) – PETSc solver options
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
solve()[source]

Computes internal pressure gradient and stores it in int_pg_3d field

## thetis.optimisation module¶

Some classes to help optimisation problems formulated with thetis_adjoint.

In particular this module contains some OptimisationCallbacks that can be used as callbacks of a ReducedFunctional called at various stages during the optimisation process: - eval_cb_pre(controls) and eval_cb_post(functional, controls) called before and after (re)evaluation of the forward model - derivative_cb_pre(controls) and eval_cb_post(functional, derivative, controls) called before and after the gradient computation using the adjoint of the model - hessian_cb_pre(controls) and eval_cb_post(functional, derivative, controls) called before and after the hessian computation OptimisationCallbacks that (can) use controls, functional and derivative information, work out what is provided by the number of arguments: current control values are always in the last argument; if more than 2 arguments are provided, the first is the latest evaluated functional value.

class thetis.optimisation.ControlsExportOptimisationCallback(solver_obj_or_outputdir, **kwargs)[source]

A callback that exports the current control values (assumed to all be Function s)

The control values are assumed to be the last argument in the callback (as for all ReducedFunctional callbacks).

Parameters: solver_obj_or_outputdir – a FlowSolver2d object, used to determine the output directory. Alternatively, the outputdir can be specified with a string as the first argument. kwargs – any further keyword arguments are passed on to UserExportManager
class thetis.optimisation.DeferredExportManager(solver_obj_or_outputdir, **kwargs)[source]

Bases: object

A wrapper around a UserExportManager that is only created on the first export() call.

In addition the functions provided in the export call are copied into a fixed set of functions, where the functions provided in subsequent calls may be different (they need to be in the same function space). This is used in the ControlsExportOptimisationCallback and DerivativesExportOptimisationCallback.

Parameters: solver_obj_or_outputdir – a FlowSolver2d object, used to determine the output directory. Alternatively, the outputdir can be specified with a string as the first argument. kwargs – any further keyword arguments are passed on to UserExportManager
export(functions, suggested_names=None)[source]

Create the UserExportManager (first call only), and call its export() method.

Parameters: functions – a list of Function s that the UserExportManager will be based on. Their values are first copied. The list may contain different functions in subsequent calls, but their function space should remain the same.
class thetis.optimisation.DerivativesExportOptimisationCallback(solver_obj_or_outputdir, **kwargs)[source]

A callback that exports the derivatives calculated by the adjoint.

The derivatives are assumed to be the second argument in the callback. This can therefore be used as a derivative_cb_post callback in a ReducedFunctional

Parameters: solver_obj_or_outputdir – a FlowSolver2d object, used to determine the output directory. Alternatively, the outputdir can be specified with a string as the first argument. kwargs – any further keyword arguments are passed on to UserExportManager
class thetis.optimisation.DiagnosticOptimisationCallback(solver_obj, **kwargs)[source]

An OptimsationCallback similar to DiagnosticCallback that can be used as callback in a ReducedFunctional.

Note that in this case the computing of the values needs to be defined in the compute_values method, not in the __call__ method (as this one is directly called from the ReducedFunctional). In addition, like any DiagnosticCallback, the name and variable_names properties and a message_str method need to be defined.

Parameters: solver_obj – Thetis solver object kwargs – keyword arguments passed to DiagnosticCallback
compute_values(*args)[source]

Compute diagnostic values.

This method is to be implemented in concrete subclasses of a DiagnosticOptimisationCallback. The number of arguments varies depending on which of the 6 [eval|derivative|hessian]_cb_[pre|post] callbacks this is used as. The last argument always contains the current controls. In the “pre” callbacks this is the only argument. In all “post” callbacks the 0th argument is the current functional value. eval_cb_post is given two arguments: functional and controls. derivative_cb_post and hessian_cb_post are given three arguments with args[1] being the derivative/hessian just calculated.

evaluate(*args, index=None)[source]

Evaluates callback and pushes values to log and hdf file (if enabled)

class thetis.optimisation.FunctionalOptimisationCallback(solver_obj, **kwargs)[source]

A simple OptimisationCallback that records the functional value in the log and/or hdf5 file.

Parameters: solver_obj – Thetis solver object kwargs – keyword arguments passed to DiagnosticCallback
compute_values(*args)[source]

Compute diagnostic values.

This method is to be implemented in concrete subclasses of a DiagnosticOptimisationCallback. The number of arguments varies depending on which of the 6 [eval|derivative|hessian]_cb_[pre|post] callbacks this is used as. The last argument always contains the current controls. In the “pre” callbacks this is the only argument. In all “post” callbacks the 0th argument is the current functional value. eval_cb_post is given two arguments: functional and controls. derivative_cb_post and hessian_cb_post are given three arguments with args[1] being the derivative/hessian just calculated.

message_str(functional)[source]

A string representation.

Parameters: args – If provided, these will be the return value from __call__().
name = 'functional'
variable_names = ['functional']
class thetis.optimisation.OptimisationCallbackList[source]

Bases: list

A list of callbacks that can be used as a single callback itself.

Calls all callbacks in order.

class thetis.optimisation.UserExportManager(solver_obj_or_outputdir, functions_to_export, filenames=None, filename_prefix='', shortnames=None, **kwargs)[source]

ExportManager for user provided functions (not necessarily known to Thetis)

In the standard ExportManager all provided functions need to have standard names present in field_metadata. Here, any functions can be provided. If function.name() is in field_metadata, the standard filename and shortname are used. If the function.name() is unknown, both are based on function.name() directly (with an optional additional filename_prefix). Filenames and shortnames can be overruled by the shortnames and filenames arguments.

Parameters: solver_obj_or_outputdir – a FlowSolver2d object, used to determine the output directory. Alternatively, the outputdir can be specified with a string as the first argument. functions_to_export – a list of Function s filenames – a list of strings that specify the filename for each provided function. If not provided, filenames are based on function.name(). filename_prefix – a string prefixed to each filename shortnames – a list of strings with the shortnames used for each provided function. If not provided, shortnames are based on function.name(). kwargs – any further keyword arguments are passed on to ExportManager
class thetis.optimisation.UserExportOptimisationCallback(solver_obj_or_outputdir, functions_to_export, **kwargs)[source]

A UserExportManager that can be used as a ReducedFunctional callback

Any callback arguments (functional value, derivatives, controls) are ignored

Parameters: solver_obj_or_outputdir – a FlowSolver2d object, used to determine the output directory. Alternatively, the outputdir can be specified with a string as the first argument. functions_to_export – a list of Function s kwargs – any further keyword arguments are passed on to UserExportManager

## thetis.options module¶

Thetis options for the 2D and 3D model

All options are type-checked and they are stored in traitlets Configurable objects.

class thetis.options.CommonModelOptions(*args, **kwargs)[source]

Options that are common for both 2d and 3d models

atmospheric_pressure
cfl_2d
cfl_3d
check_volume_conservation_2d

A boolean (True, False) trait.

coriolis_frequency
element_family

An enum whose value must be in a given sequence.

export_diagnostics

A boolean (True, False) trait.

fields_to_export

An instance of a Python list.

fields_to_export_hdf5

An instance of a Python list.

horizontal_diffusivity
horizontal_velocity_scale
horizontal_viscosity
horizontal_viscosity_scale
lax_friedrichs_tracer_scaling_factor
lax_friedrichs_velocity_scaling_factor
linear_drag_coefficient
log_output

A boolean (True, False) trait.

manning_drag_coefficient
momentum_source_2d
name = 'Model options'
no_exports

A boolean (True, False) trait.

output_directory

A trait for unicode strings.

polynomial_degree
quadratic_drag_coefficient
simulation_end_time
simulation_export_time
timestep
tracer_source_2d
use_grad_depth_viscosity_term

A boolean (True, False) trait.

use_grad_div_viscosity_term

A boolean (True, False) trait.

use_lax_friedrichs_tracer

A boolean (True, False) trait.

use_lax_friedrichs_velocity

A boolean (True, False) trait.

use_limiter_for_tracers

A boolean (True, False) trait.

use_nonlinear_equations

A boolean (True, False) trait.

verbose

An int trait.

volume_source_2d
wind_stress
class thetis.options.CrankNicolsonTimestepperOptions2d(*args, **kwargs)[source]

Options for 2d Crank-Nicolson time integrator

implicitness_theta
use_semi_implicit_linearization

A boolean (True, False) trait.

class thetis.options.EquationOfStateOptions(*args, **kwargs)[source]

Base class of equation of state options

name = 'Equation of State'
class thetis.options.ExplicitTimestepperOptions(*args, **kwargs)[source]

Options for explicit time integrator

use_automatic_timestep

A boolean (True, False) trait.

class thetis.options.ExplicitTimestepperOptions2d(*args, **kwargs)[source]

Options for 2d explicit time integrator

solver_parameters

PETSc solver options dictionary

class thetis.options.ExplicitTimestepperOptions3d(*args, **kwargs)[source]

Base class for all 3d time stepper options

solver_parameters_2d_swe

PETSc solver options dictionary

solver_parameters_momentum_explicit

PETSc solver options dictionary

solver_parameters_momentum_implicit

PETSc solver options dictionary

solver_parameters_tracer_explicit

PETSc solver options dictionary

solver_parameters_tracer_implicit

PETSc solver options dictionary

class thetis.options.GLSModelOptions(*args, **kwargs)[source]

Options for Generic Length Scale turbulence model

apply_defaults(closure_name)[source]

Applies default parameters for given closure name

Parameters: closure_name (string) – name of the turbulence closure model

Sets default values for parameters p, m, n, schmidt_nb_tke, schmidt_nb_psi, c1, c2, c3_plus, c3_minus, f_wall, k_min, psi_min

c1

A float trait.

c2

A float trait.

c3_minus

A float trait.

c3_plus

A float trait.

closure_name

An enum whose value must be in a given sequence.

cmu0
compute_c3_minus

A boolean (True, False) trait.

compute_cmu0

A boolean (True, False) trait.

compute_kappa

A boolean (True, False) trait.

compute_len_min

A boolean (True, False) trait.

compute_psi_min

A boolean (True, False) trait.

diff_min
eps_min
f_wall

A float trait.

galperin_lim
k_min
kappa

A float trait.

len_min
limit_eps

A boolean (True, False) trait.

limit_len

A boolean (True, False) trait.

limit_len_min

A boolean (True, False) trait.

limit_psi

A boolean (True, False) trait.

m

A float trait.

n

A float trait.

name = 'Generic Lenght Scale turbulence closure model'
p

A float trait.

psi_min
ri_st

A float trait.

schmidt_nb_psi
schmidt_nb_tke
stability_function_name

An enum whose value must be in a given sequence.

visc_min
class thetis.options.LinearEquationOfStateOptions(*args, **kwargs)[source]

Linear equation of state options

alpha

A float trait.

beta

A float trait.

name = 'Linear Equation of State'
rho_ref
s_ref
th_ref

A float trait.

class thetis.options.ModelOptions2d(*args, **kwargs)[source]

Options for 2D depth-averaged shallow water model

check_tracer_conservation

A boolean (True, False) trait.

check_tracer_overshoot

A boolean (True, False) trait.

name = 'Depth-averaged 2D model'
solve_tracer

A boolean (True, False) trait.

tidal_turbine_farms

An instance of a Python dict.

timestepper_options

A trait whose value must be an instance of a specified class.

The value can also be an instance of a subclass of the specified class.

Subclasses can declare default classes by overriding the klass attribute

timestepper_type

A enum whose value must be in a given sequence.

This enum controls a slaved option, with default values provided here.

Parameters: values – iterable of (value, HasTraits class) pairs. The HasTraits class will be called (with no arguments) to create default values if necessary. paired_name – trait name this enum is paired with. default_value – default value.
tracer_only

A boolean (True, False) trait.

use_wetting_and_drying

A boolean (True, False) trait.

wetting_and_drying_alpha
class thetis.options.ModelOptions3d(*args, **kwargs)[source]

Options for 3D hydrostatic model

check_salinity_conservation

A boolean (True, False) trait.

check_salinity_overshoot

A boolean (True, False) trait.

check_temperature_conservation

A boolean (True, False) trait.

check_temperature_overshoot

A boolean (True, False) trait.

check_volume_conservation_3d

A boolean (True, False) trait.

constant_salinity
constant_temperature
equation_of_state_options

A trait whose value must be an instance of a specified class.

The value can also be an instance of a subclass of the specified class.

Subclasses can declare default classes by overriding the klass attribute

equation_of_state_type

A enum whose value must be in a given sequence.

This enum controls a slaved option, with default values provided here.

Parameters: values – iterable of (value, HasTraits class) pairs. The HasTraits class will be called (with no arguments) to create default values if necessary. paired_name – trait name this enum is paired with. default_value – default value.
horizontal_diffusivity
momentum_source_3d
name = '3D hydrostatic model'
salinity_source_3d
smagorinsky_coefficient
solve_salinity

A boolean (True, False) trait.

solve_temperature

A boolean (True, False) trait.

temperature_source_3d
timestep_2d
timestepper_options

A trait whose value must be an instance of a specified class.

The value can also be an instance of a subclass of the specified class.

Subclasses can declare default classes by overriding the klass attribute

timestepper_type

A enum whose value must be in a given sequence.

This enum controls a slaved option, with default values provided here.

Parameters: values – iterable of (value, HasTraits class) pairs. The HasTraits class will be called (with no arguments) to create default values if necessary. paired_name – trait name this enum is paired with. default_value – default value.
turbulence_model_options

A trait whose value must be an instance of a specified class.

The value can also be an instance of a subclass of the specified class.

Subclasses can declare default classes by overriding the klass attribute

turbulence_model_type

A enum whose value must be in a given sequence.

This enum controls a slaved option, with default values provided here.

Parameters: values – iterable of (value, HasTraits class) pairs. The HasTraits class will be called (with no arguments) to create default values if necessary. paired_name – trait name this enum is paired with. default_value – default value.
use_ale_moving_mesh

A boolean (True, False) trait.

use_baroclinic_formulation

A boolean (True, False) trait.

use_bottom_friction

A boolean (True, False) trait.

use_implicit_vertical_diffusion

A boolean (True, False) trait.

use_limiter_for_velocity

A boolean (True, False) trait.

use_parabolic_viscosity

A boolean (True, False) trait.

use_quadratic_density

A boolean (True, False) trait.

use_quadratic_pressure

A boolean (True, False) trait.

use_smagorinsky_viscosity

A boolean (True, False) trait.

use_smooth_eddy_viscosity

A boolean (True, False) trait.

use_turbulence

A boolean (True, False) trait.

use_turbulence_advection

A boolean (True, False) trait.

vertical_diffusivity
vertical_velocity_scale
vertical_viscosity
class thetis.options.PacanowskiPhilanderModelOptions(*args, **kwargs)[source]

Options for Pacanowski-Philander turbulence model

alpha
exponent
max_viscosity
name = 'Pacanowski-Philander turbulence closure model'
class thetis.options.PressureProjectionTimestepperOptions2d(*args, **kwargs)[source]

Options for 2d pressure-projection time integrator

implicitness_theta
picard_iterations
solver_parameters_momentum

PETSc solver options dictionary

solver_parameters_pressure

PETSc solver options dictionary

use_semi_implicit_linearization

A boolean (True, False) trait.

class thetis.options.SemiImplicitTimestepperOptions2d(*args, **kwargs)[source]

Options for 2d explicit time integrator

solver_parameters

PETSc solver options dictionary

solver_parameters_tracer

PETSc solver options dictionary

class thetis.options.SemiImplicitTimestepperOptions3d(*args, **kwargs)[source]

Class for all 3d time steppers that have a configurable semi-implicit 2D solver

implicitness_theta_2d
class thetis.options.SteadyStateTimestepperOptions2d(*args, **kwargs)[source]

Options for 2d steady state solver

solver_parameters

PETSc solver options dictionary

class thetis.options.TidalTurbineFarmOptions(*args, **kwargs)[source]

Bases: thetis.configuration.FrozenHasTraits, traitlets.traitlets.TraitType

Tidal turbine farm options

break_even_wattage
name = 'Farm options'
turbine_density
turbine_options = <thetis.options.TidalTurbineOptions object>
class thetis.options.TidalTurbineOptions(*args, **kwargs)[source]

Tidal turbine parameters

diameter
thrust_coefficient
class thetis.options.TimeStepperOptions(*args, **kwargs)[source]

Base class for all time stepper options

name = 'Time stepper'
class thetis.options.TurbulenceModelOptions(*args, **kwargs)[source]

Abstract base class for all turbulence model options

name = 'Turbulence closure model'

## thetis.physical_constants module¶

Default values for physical constants and parameters

## thetis.rungekutta module¶

Implements Runge-Kutta time integration methods.

The abstract class AbstractRKScheme defines the Runge-Kutta coefficients, and can be used to implement generic time integrators.

class thetis.rungekutta.AbstractRKScheme[source]

Bases: object

Abstract class for defining Runge-Kutta schemes.

Derived classes must define the Butcher tableau (arrays a, b, c) and the CFL number (cfl_coeff).

Currently only explicit or diagonally implicit schemes are supported.

a

Runge-Kutta matrix $$a_{i,j}$$ of the Butcher tableau

b

weights $$b_{i}$$ of the Butcher tableau

c

nodes $$c_{i}$$ of the Butcher tableau

cfl_coeff

CFL number of the scheme

Value 1.0 corresponds to Forward Euler time step.

class thetis.rungekutta.BackwardEuler(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.BackwardEulerAbstract[source]

Backward Euler method

a = [[1.0]]
b = [1.0]
c = [1.0]
cfl_coeff = inf
class thetis.rungekutta.CrankNicolsonAbstract[source]

Crack-Nicolson scheme

a = [[0.0, 0.0], [0.5, 0.5]]
b = [0.5, 0.5]
c = [0.0, 1.0]
cfl_coeff = inf
class thetis.rungekutta.CrankNicolsonRK(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.DIRK22(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.DIRK22Abstract[source]

2-stage, 2nd order, L-stable Diagonally Implicit Runge Kutta method

This method has the Butcher tableau

$\begin{split}\begin{array}{c|cc} \gamma & \gamma & 0 \\ 1 & 1-\gamma & \gamma \\ \hline & 1/2 & 1/2 \end{array}\end{split}$

with $$\gamma = (2 + \sqrt{2})/2$$.

From DIRK(2,3,2) IMEX scheme in Ascher et al. (1997)

Ascher et al. (1997). Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Applied Numerical Mathematics, 25:151-167. http://dx.doi.org/10.1137/0732037

a = [[1.7071067811865475, 0], [-0.7071067811865475, 1.7071067811865475]]
b = [-0.7071067811865475, 1.7071067811865475]
c = [1.7071067811865475, 1]
cfl_coeff = inf
gamma = 1.7071067811865475
class thetis.rungekutta.DIRK23(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.DIRK23Abstract[source]

2-stage, 3rd order Diagonally Implicit Runge Kutta method

This method has the Butcher tableau

$\begin{split}\begin{array}{c|cc} \gamma & \gamma & 0 \\ 1-\gamma & 1-2\gamma & \gamma \\ \hline & 1/2 & 1/2 \end{array}\end{split}$

with $$\gamma = (3 + \sqrt{3})/6$$.

From DIRK(2,3,3) IMEX scheme in Ascher et al. (1997)

Ascher et al. (1997). Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Applied Numerical Mathematics, 25:151-167. http://dx.doi.org/10.1137/0732037

a = [[0.7886751345948128, 0], [-0.5773502691896255, 0.7886751345948128]]
b = [0.5, 0.5]
c = [0.7886751345948128, 0.21132486540518725]
cfl_coeff = inf
gamma = 0.7886751345948128
class thetis.rungekutta.DIRK33(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.DIRK33Abstract[source]

3-stage, 3rd order, L-stable Diagonally Implicit Runge Kutta method

From DIRK(3,4,3) IMEX scheme in Ascher et al. (1997)

Ascher et al. (1997). Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Applied Numerical Mathematics, 25:151-167. http://dx.doi.org/10.1137/0732037

a = [[0.4358665215, 0, 0], [0.28206673925000003, 0.4358665215, 0], [1.208496649153235, -0.6443631706532353, 0.4358665215]]
b = [1.208496649153235, -0.6443631706532353, 0.4358665215]
b1 = 1.208496649153235
b2 = -0.6443631706532353
c = [0.4358665215, 0.71793326075, 1]
cfl_coeff = inf
gamma = 0.4358665215
class thetis.rungekutta.DIRK43(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.DIRK43Abstract[source]

4-stage, 3rd order, L-stable Diagonally Implicit Runge Kutta method

From DIRK(4,4,3) IMEX scheme in Ascher et al. (1997)

Ascher et al. (1997). Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Applied Numerical Mathematics, 25:151-167. http://dx.doi.org/10.1137/0732037

a = [[0.5, 0, 0, 0], [0.16666666666666666, 0.5, 0, 0], [-0.5, 0.5, 0.5, 0], [1.5, -1.5, 0.5, 0.5]]
b = [1.5, -1.5, 0.5, 0.5]
c = [0.5, 0.6666666666666666, 0.5, 1.0]
cfl_coeff = inf
class thetis.rungekutta.DIRKGeneric(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]

Generic implementation of Diagonally Implicit Runge Kutta schemes.

All derived classes must define the Butcher tableau coefficients a, b, c.

Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
get_final_solution()[source]

Assign final solution to self.solution

initialize(init_cond)[source]

solve_stage(i_stage, t, update_forcings=None)[source]

Solve i-th stage and assign solution to self.solution.

solve_tendency(i_stage, t, update_forcings=None)[source]

Evaluates the tendency of i-th stage.

update_solution(i_stage)[source]

Tendencies must have been evaluated first.

update_solver()[source]

Create solver objects

class thetis.rungekutta.DIRKLPUM2(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.DIRKLPUM2Abstract[source]

DIRKLPUM2, 3-stage, 2nd order, L-stable Diagonally Implicit Runge Kutta method

From IMEX RK scheme (20) in Higureras et al. (2014).

Higueras et al (2014). Optimized strong stability preserving IMEX Runge-Kutta methods. Journal of Computational and Applied Mathematics 272(2014) 116-140. http://dx.doi.org/10.1016/j.cam.2014.05.011

a = [[0.18181818181818182, 0, 0], [0.2662337662337662, 0.18181818181818182, 0], [0.3412042502951594, 0.34710743801652894, 0.18181818181818182]]
b = [0.3333333333333333, 0.3333333333333333, 0.3333333333333333]
c = [0.18181818181818182, 0.44805194805194803, 0.8701298701298701]
cfl_coeff = 4.34
class thetis.rungekutta.DIRKLSPUM2(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.DIRKLSPUM2Abstract[source]

DIRKLSPUM2, 3-stage, 2nd order, L-stable Diagonally Implicit Runge Kutta method

From IMEX RK scheme (17) in Higureras et al. (2014).

Higueras et al (2014). Optimized strong stability preserving IMEX Runge-Kutta methods. Journal of Computational and Applied Mathematics 272(2014) 116-140. http://dx.doi.org/10.1016/j.cam.2014.05.011

a = [[0.18181818181818182, 0, 0], [0.44372294372294374, 0.18181818181818182, 0], [0.44004329004329007, 0.19090909090909092, 0.18181818181818182]]
b = [0.43636363636363634, 0.2, 0.36363636363636365]
c = [0.18181818181818182, 0.6255411255411255, 0.8127705627705628]
cfl_coeff = 4.34
class thetis.rungekutta.ERKEuler(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.ERKGeneric(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]

Generic explicit Runge-Kutta time integrator.

Implements the Butcher form. All terms in the equation are treated explicitly.

Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
get_final_solution(additive=False)[source]

Assign final solution to self.solution

If additive=False, will overwrite solution function, otherwise will add to it.

initialize(solution)[source]

solve_stage(i_stage, t, update_forcings=None)[source]

Solve i-th stage and assign solution to self.solution.

solve_tendency(i_stage, t, update_forcings=None)[source]

Evaluates the tendency of i-th stage

update_solution(i_stage, additive=False)[source]

Computes the solution of the i-th stage

Tendencies must have been evaluated first.

If additive=False, will overwrite solution function, otherwise will add to it.

update_solver()[source]
class thetis.rungekutta.ERKGenericShuOsher(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]

Generic explicit Runge-Kutta time integrator.

Implements the Shu-Osher form.

Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
advance(t, update_forcings=None)[source]

Advances equations for one time step.

initialize(solution)[source]

Initialize the time integrator

Parameters: init_solution – initial solution
solve_stage(i_stage, t, update_forcings=None)[source]

Solve i-th stage and assign solution to self.solution.

update_solver()[source]
class thetis.rungekutta.ERKLPUM2(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.ERKLPUM2Abstract[source]

ERKLPUM2, 3-stage, 2nd order Explicit Runge Kutta method

From IMEX RK scheme (20) in Higureras et al. (2014).

Higueras et al (2014). Optimized strong stability preserving IMEX Runge-Kutta methods. Journal of Computational and Applied Mathematics 272(2014) 116-140. http://dx.doi.org/10.1016/j.cam.2014.05.011

a = [[0, 0, 0], [0.5, 0, 0], [0.5, 0.5, 0]]
b = [0.3333333333333333, 0.3333333333333333, 0.3333333333333333]
c = [0, 0.5, 1.0]
cfl_coeff = 2.0
class thetis.rungekutta.ERKLSPUM2(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.ERKLSPUM2Abstract[source]

ERKLSPUM2, 3-stage, 2nd order Explicit Runge Kutta method

From IMEX RK scheme (17) in Higureras et al. (2014).

Higueras et al (2014). Optimized strong stability preserving IMEX Runge-Kutta methods. Journal of Computational and Applied Mathematics 272(2014) 116-140. http://dx.doi.org/10.1016/j.cam.2014.05.011

a = [[0, 0, 0], [0.8333333333333334, 0, 0], [0.4583333333333333, 0.4583333333333333, 0]]
b = [0.43636363636363634, 0.2, 0.36363636363636365]
c = [0, 0.8333333333333334, 0.9166666666666666]
cfl_coeff = 1.2
class thetis.rungekutta.ERKMidpoint(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.ERKMidpointAbstract[source]
a = [[0.0, 0.0], [0.5, 0.0]]
b = [0.0, 1.0]
c = [0.0, 0.5]
cfl_coeff = 1.0
class thetis.rungekutta.ESDIRKMidpoint(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.ESDIRKMidpointAbstract[source]
a = [[0.0, 0.0], [0.0, 0.5]]
b = [0.0, 1.0]
c = [0.0, 0.5]
cfl_coeff = 1.0
class thetis.rungekutta.ESDIRKTrapezoid(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.ESDIRKTrapezoidAbstract[source]
a = [[0.0, 0.0], [0.5, 0.5]]
b = [0.5, 0.5]
c = [0.0, 1.0]
cfl_coeff = inf
class thetis.rungekutta.ForwardEulerAbstract[source]

Forward Euler method

a = [[0]]
b = [1.0]
c = [0]
cfl_coeff = 1.0
class thetis.rungekutta.ImplicitMidpoint(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.ImplicitMidpointAbstract[source]

Implicit midpoint method, second order.

This method has the Butcher tableau

$\begin{split}\begin{array}{c|c} 0.5 & 0.5 \\ \hline & 1.0 \end{array}\end{split}$
a = [[0.5]]
b = [1.0]
c = [0.5]
cfl_coeff = inf
class thetis.rungekutta.RungeKuttaTimeIntegrator(equation, solution, fields, dt, solver_parameters={})[source]

Abstract base class for all Runge-Kutta time integrators

Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds solver_parameters (dict) – PETSc solver options
advance(t, update_forcings=None)[source]

Advances equations for one time step.

get_final_solution(additive=False)[source]

Evaluates the final solution

solve_stage(i_stage, t, update_forcings=None)[source]

Solves a single stage of step from t to t+dt. All functions that the equation depends on must be at right state corresponding to each sub-step.

class thetis.rungekutta.SSPRK33(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]
Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
class thetis.rungekutta.SSPRK33Abstract[source]

3rd order Strong Stability Preserving Runge-Kutta scheme, SSP(3,3).

This scheme has Butcher tableau

$\begin{split}\begin{array}{c|ccc} 0 & \\ 1 & 1 \\ 1/2 & 1/4 & 1/4 & \\ \hline & 1/6 & 1/6 & 2/3 \end{array}\end{split}$

CFL coefficient is 1.0

a = [[0, 0, 0], [1.0, 0, 0], [0.25, 0.25, 0]]
b = [0.16666666666666666, 0.16666666666666666, 0.6666666666666666]
c = [0, 1.0, 0.5]
cfl_coeff = 1.0
thetis.rungekutta.butcher_to_shuosher_form(a, b)[source]

Converts Butcher tableau to Shu-Osher form.

The Shu-Osher form of a s-stage scheme is defined by two s+1 by s+1 arrays $$\alpha$$ and $$\beta$$:

$\begin{split}u^{0} &= u^n \\ u^{(i)} &= \sum_{j=0}^s \alpha_{i,j} u^{(j)} + \sum_{j=0}^s \beta_{i,j} F(u^{(j)}) \\ u^{n+1} &= u^{(s)}\end{split}$

The Shu-Osher form is not unique. Here we construct the form where beta values are the diagonal entries (for DIRK schemes) or sub-diagonal entries (for explicit schemes) of the concatenated Butcher tableau [$$a$$; $$b$$].

thetis.rungekutta.timed_region()

Log.Event(type cls, name, klass=None)

thetis.rungekutta.timed_stage()

Log.Stage(type cls, name)

## thetis.shallowwater_eq module¶

Depth averaged shallow water equations

### Equations¶

The state variables are water elevation, $$\eta$$, and depth averaged velocity $$\bar{\textbf{u}}$$.

Denoting the total water depth by $$H=\eta + h$$, the non-conservative form of the free surface equation is

(5)$\frac{\partial \eta}{\partial t} + \nabla \cdot (H \bar{\textbf{u}}) = 0$

(6)$\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 ),$

where $$g$$ is the gravitational acceleration, $$f$$ is the Coriolis frequency, $$\wedge$$ is the cross product, $$\textbf{e}_z$$ is a vertical unit vector, $$p_a$$ is the atmospheric pressure at the free surface, and $$\nu_h$$ is viscosity. Water density is given by $$\rho = \rho'(T, S, p) + \rho_0$$, where $$\rho_0$$ is a constant reference density.

Above $$r$$ denotes the baroclinic head

$r = \frac{1}{\rho_0} \int_{z}^\eta \rho' d\zeta.$

In the case of purely barotropic problems the $$r$$ and the internal pressure gradient are omitted.

If the option ModelOptions.use_nonlinear_equations is False, we solve the linear shallow water equations (i.e. wave equation):

(7)$\frac{\partial \eta}{\partial t} + \nabla \cdot (h \bar{\textbf{u}}) = 0$
(8)$\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 ).$

In case of a 3D problem with mode splitting, we use a simplified 2D system that contains nothing but the rotational external gravity waves:

(9)$\frac{\partial \eta}{\partial t} + \nabla \cdot (H \bar{\textbf{u}}) = 0$
(10)$\frac{\partial \bar{\textbf{u}}}{\partial t} + f\textbf{e}_z\wedge \bar{\textbf{u}} + g \nabla \eta = \textbf{G},$

where $$\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 $$\eta$$ and $$\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

swe_bnd_funcs = {}
swe_bnd_funcs[1] = {'elev':myfunc1, 'flux':myfunc2}


where myfunc1 and myfunc2 are Constant or Function objects.

The user can provide $$\eta$$ and/or $$\bar{\textbf{u}}$$ values. Supported combinations are:

• unspecified : impermeable (land) boundary, implies symmetric $$\eta$$ condition and zero normal velocity
• 'elev': elevation only, symmetric velocity (usually unstable)
• 'uv': 2d velocity vector $$\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 FlowSolver2d or FlowSolver objects:

solver_obj = solver2d.FlowSolver2d(...)
...
solver_obj.bnd_functions['shallow_water'] = swe_bnd_funcs


Internally the boundary conditions passed to the Term.residual() method of each term:

adv_term = shallowwater_eq.HorizontalAdvectionTerm(...)


### Wetting and drying¶

If the option 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 $$\tilde{h} = h + f(H)$$, which ensures positive total water depth, with $$f(H)$$ defined by

$f(H) = \frac{1}{2}(\sqrt{H^2 + \alpha^2} - H),$

introducing a wetting-drying parameter $$\alpha$$, with dimensions of length. This results in a modified total water depth $$\tilde{H}=H+f(H)$$.

The value for $$\alpha$$ is specified by the user through the option ModelOptions.wetting_and_drying_alpha, in units of meters. The default value for 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 $$\alpha$$ is suggested by Karna et al. (2011). Defining $$L_x$$ as the horizontal length scale of the mesh elements at the wet-dry front, it can be reasoned that $$\alpha \approx |L_x \nabla h|$$ yields a suitable choice. Smaller $$\alpha$$ leads to a more accurate solution to the shallow water equations in wet regions, but if $$\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, $$H$$, are replaced by $$\tilde{H}$$ (as defined above);
2. An additional displacement term $$\frac{\partial \tilde{h}}{\partial t}$$ is added to the bathymetry in the free surface equation.

The free surface and momentum equations then become:

(11)$\frac{\partial \eta}{\partial t} + \frac{\partial \tilde{h}}{\partial t} + \nabla \cdot (\tilde{H} \bar{\textbf{u}}) = 0,$
(12)$\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 ).$
class thetis.shallowwater_eq.BaseShallowWaterEquation(function_space, bathymetry, options)[source]

Abstract base class for ShallowWaterEquations, ShallowWaterMomentumEquation and FreeSurfaceEquation.

Provides common functionality to compute time steps and add either momentum or continuity terms.

add_continuity_terms(*args)[source]
add_momentum_terms(*args)[source]
residual_uv_eta(label, uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions)[source]
class thetis.shallowwater_eq.ShallowWaterEquations(function_space, bathymetry, options)[source]

2D depth-averaged shallow water equations in non-conservative form.

This defines the full 2D SWE equations (5) - (6).

Parameters: function_space – Mixed function space where the solution belongs bathymetry (Function or Constant) – bathymetry of the domain options – AttrDict object containing all circulation model options
mass_term(solution)[source]

Returns default mass matrix term for the solution function space.

Returns: UFL form of the mass term
residual(label, solution, solution_old, fields, fields_old, bnd_conditions)[source]

Returns an UFL form of the residual by summing up all the terms with the desired label.

Sign convention: all terms are assumed to be on the right hand side of the equation: d(u)/dt = term.

Parameters: label – string defining the type of terms to sum up. Currently one of ‘source’|’explicit’|’implicit’|’nonlinear’. Can be a list of multiple labels, or ‘all’ in which case all defined terms are summed. solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.ModeSplit2DEquations(function_space, bathymetry, options)[source]

2D depth-averaged shallow water equations for mode splitting schemes.

Defines the equations (9) - (10).

Parameters: function_space – Mixed function space where the solution belongs bathymetry (Function or Constant) – bathymetry of the domain options – AttrDict object containing all circulation model options
add_momentum_terms(*args)[source]
residual(label, solution, solution_old, fields, fields_old, bnd_conditions)[source]

Returns an UFL form of the residual by summing up all the terms with the desired label.

Sign convention: all terms are assumed to be on the right hand side of the equation: d(u)/dt = term.

Parameters: label – string defining the type of terms to sum up. Currently one of ‘source’|’explicit’|’implicit’|’nonlinear’. Can be a list of multiple labels, or ‘all’ in which case all defined terms are summed. solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.ShallowWaterMomentumEquation(u_test, u_space, eta_space, bathymetry, options)[source]

2D depth averaged momentum equation (6) in non-conservative form.

Parameters: u_test – test function of the velocity function space u_space – velocity function space eta_space – elevation function space bathymetry (Function or Constant) – bathymetry of the domain options – AttrDict object containing all circulation model options
residual(label, solution, solution_old, fields, fields_old, bnd_conditions)[source]

Returns an UFL form of the residual by summing up all the terms with the desired label.

Sign convention: all terms are assumed to be on the right hand side of the equation: d(u)/dt = term.

Parameters: label – string defining the type of terms to sum up. Currently one of ‘source’|’explicit’|’implicit’|’nonlinear’. Can be a list of multiple labels, or ‘all’ in which case all defined terms are summed. solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.FreeSurfaceEquation(eta_test, eta_space, u_space, bathymetry, options)[source]

2D free surface equation (5) in non-conservative form.

Parameters: eta_test – test function of the elevation function space eta_space – elevation function space u_space – velocity function space function_space – Mixed function space where the solution belongs bathymetry (Function or Constant) – bathymetry of the domain options – AttrDict object containing all circulation model options
mass_term(solution)[source]

Returns default mass matrix term for the solution function space.

Returns: UFL form of the mass term
residual(label, solution, solution_old, fields, fields_old, bnd_conditions)[source]

Returns an UFL form of the residual by summing up all the terms with the desired label.

Sign convention: all terms are assumed to be on the right hand side of the equation: d(u)/dt = term.

Parameters: label – string defining the type of terms to sum up. Currently one of ‘source’|’explicit’|’implicit’|’nonlinear’. Can be a list of multiple labels, or ‘all’ in which case all defined terms are summed. solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.ShallowWaterTerm(space, bathymetry=None, options=None)[source]

Generic term in the shallow water equations that provides commonly used members and mapping for boundary functions.

get_bnd_functions(eta_in, uv_in, bnd_id, bnd_conditions)[source]

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.

get_total_depth(eta)[source]

Returns total water column depth

wd_bathymetry_displacement(eta)[source]

Returns wetting and drying bathymetry displacement as described in: Karna et al., 2011.

class thetis.shallowwater_eq.ShallowWaterMomentumTerm(u_test, u_space, eta_space, bathymetry=None, options=None)[source]

Generic term in the shallow water momentum equation that provides commonly used members and mapping for boundary functions.

class thetis.shallowwater_eq.ShallowWaterContinuityTerm(eta_test, eta_space, u_space, bathymetry=None, options=None)[source]

Generic term in the depth-integrated continuity equation that provides commonly used members and mapping for boundary functions.

class thetis.shallowwater_eq.HUDivTerm(eta_test, eta_space, u_space, bathymetry=None, options=None)[source]

Divergence term, $$\nabla \cdot (H \bar{\textbf{u}})$$

$\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; $$\textbf{n}$$ denotes the unit normal of the element interfaces, and $$\text{jump}$$ and $$\text{avg}$$ denote the jump and average operators across the interface. $$H^*, \bar{\textbf{u}}^*$$ are values at the interface obtained from an approximate Riemann solver.

If $$\bar{\textbf{u}}$$ belongs to a discontinuous function space, the form on the right hand side is used.

residual(uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.ContinuitySourceTerm(eta_test, eta_space, u_space, bathymetry=None, options=None)[source]

Generic source term in the depth-averaged continuity equation

$F_s = \int_\Omega S \phi dx$

where $$S$$ is a user defined scalar Function.

residual(uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.HorizontalAdvectionTerm(u_test, u_space, eta_space, bathymetry=None, options=None)[source]

Advection of momentum term, $$\bar{\textbf{u}} \cdot \nabla\bar{\textbf{u}}$$

The weak form is

$\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; $$\textbf{n}$$ is the unit normal of the element interfaces, and $$\text{jump}$$ and $$\text{avg}$$ denote the jump and average operators across the interface.

residual(uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.HorizontalViscosityTerm(u_test, u_space, eta_space, bathymetry=None, options=None)[source]

Viscosity of momentum term

If option ModelOptions.use_grad_div_viscosity_term is True, we use the symmetric viscous stress $$\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

$\begin{split}\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\end{split}$

where $$\sigma$$ is a penalty parameter, see Epshteyn and Riviere (2007).

If option ModelOptions.use_grad_div_viscosity_term is False, we use viscous stress $$\boldsymbol{\tau}_\nu = \nu_h \nabla \bar{\textbf{u}}$$. In this case the weak form is

$\begin{split}\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\end{split}$

If option ModelOptions.use_grad_depth_viscosity_term is True, we also include the term

$\boldsymbol{\tau}_{\nabla H} = - \frac{\nu_h \nabla(H)}{H} \cdot ( \nabla \bar{\textbf{u}} + (\nabla \bar{\textbf{u}})^T )$

as a source term.

Epshteyn and Riviere (2007). Estimation of penalty parameters for symmetric interior penalty Galerkin methods. Journal of Computational and Applied Mathematics, 206(2):843-872. http://dx.doi.org/10.1016/j.cam.2006.08.029

residual(uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.ExternalPressureGradientTerm(u_test, u_space, eta_space, bathymetry=None, options=None)[source]

External pressure gradient term, $$g \nabla \eta$$

$\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; $$\textbf{n}$$ denotes the unit normal of the element interfaces, $$n^*$$ is value at the interface obtained from an approximate Riemann solver.

If $$\eta$$ belongs to a discontinuous function space, the form on the right hand side is used.

residual(uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.CoriolisTerm(u_test, u_space, eta_space, bathymetry=None, options=None)[source]

Coriolis term, $$f\textbf{e}_z\wedge \bar{\textbf{u}}$$

residual(uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.LinearDragTerm(u_test, u_space, eta_space, bathymetry=None, options=None)[source]

Linear friction term, $$C \bar{\textbf{u}}$$

Here $$C$$ is a user-defined drag coefficient.

residual(uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.QuadraticDragTerm(u_test, u_space, eta_space, bathymetry=None, options=None)[source]

Quadratic Manning bottom friction term $$C_D \| \bar{\textbf{u}} \| \bar{\textbf{u}}$$

where the drag term is computed with the Manning formula

$C_D = g \frac{\mu^2}{H^{1/3}}$

if the Manning coefficient $$\mu$$ is defined (see field manning_drag_coefficient). Otherwise $$C_D$$ is taken as a constant (see field quadratic_drag_coefficient).

residual(uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.BottomDrag3DTerm(u_test, u_space, eta_space, bathymetry=None, options=None)[source]

Bottom drag term consistent with the 3D mode, $$C_D \| \textbf{u}_b \| \textbf{u}_b$$

Here $$\textbf{u}_b$$ is the bottom velocity used in the 3D mode, and $$C_D$$ the corresponding bottom drag. These fields are computed in the 3D model.

residual(uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.MomentumSourceTerm(u_test, u_space, eta_space, bathymetry=None, options=None)[source]

Generic source term in the shallow water momentum equation

$F_s = \int_\Omega \boldsymbol{\tau} \cdot \boldsymbol{\psi} dx$

where $$\boldsymbol{\tau}$$ is a user defined vector valued Function.

residual(uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.WindStressTerm(u_test, u_space, eta_space, bathymetry=None, options=None)[source]

Wind stress term, $$-\tau_w/(H \rho_0)$$

Here $$\tau_w$$ is a user-defined wind stress Function.

residual(uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.shallowwater_eq.AtmosphericPressureTerm(u_test, u_space, eta_space, bathymetry=None, options=None)[source]

Atmospheric pressure term, $$\nabla (p_a / \rho_0)$$

Here $$p_a$$ is a user-defined atmospheric pressure Function.

residual(uv, eta, uv_old, eta_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.

## thetis.solver module¶

Module for three dimensional baroclinic solver

class thetis.solver.FlowSolver(mesh2d, bathymetry_2d, n_layers, options=None, extrude_options=None)[source]

Main object for 3D solver

Example

Create 2D mesh

from thetis import *
mesh2d = RectangleMesh(20, 20, 10e3, 10e3)


Create bathymetry function, set a constant value

fs_p1 = FunctionSpace(mesh2d, 'CG', 1)
bathymetry_2d = Function(fs_p1, name='Bathymetry').assign(10.0)


Create a 3D model with 6 uniform levels, and set some options (see ModelOptions3d)

solver_obj = solver.FlowSolver(mesh2d, bathymetry_2d, n_layers=6)
options = solver_obj.options
options.element_family = 'dg-dg'
options.polynomial_degree = 1
options.timestepper_type = 'SSPRK22'
options.timestepper_options.use_automatic_timestep = False
options.solve_salinity = False
options.solve_temperature = False
options.simulation_export_time = 50.0
options.simulation_end_time = 3600.
options.timestep = 25.0


solver_obj.create_function_spaces()
init_elev = Function(solver_obj.function_spaces.H_2d)
coords = SpatialCoordinate(mesh2d)
init_elev.project(2.0*exp(-((coords[0] - 4e3)**2 + (coords[1] - 4.5e3)**2)/2.2e3**2))
solver_obj.assign_initial_conditions(elev=init_elev)


Run simulation

solver_obj.iterate()


See the manual for more complex examples.

Parameters: mesh2d – Mesh object of the 2D mesh bathymetry_2d (2D Function) – Bathymetry of the domain. Bathymetry stands for the mean water depth (positive downwards). n_layers (int) – Number of layers in the vertical direction. Elements are distributed uniformly over the vertical. options (ModelOptions3d instance) – Model options (optional). Model options can also be changed directly via the options class property.
M_modesplit = None

Mode split ratio (int)

add_callback(callback, eval_interval='export')[source]

Parameters: callback – DiagnosticCallback instance eval_interval (str) – Determines when callback will be evaluated, either ‘export’ or ‘timestep’ for evaluating after each export or time step.
assign_initial_conditions(elev=None, salt=None, temp=None, uv_2d=None, uv_3d=None, tke=None, psi=None)[source]

Assigns initial conditions

Parameters: elev (scalar 2D Function, Constant, or an expression) – Initial condition for water elevation salt (scalar 3D Function, Constant, or an expression) – Initial condition for salinity field temp (scalar 3D Function, Constant, or an expression) – Initial condition for temperature field uv_2d (vector valued 2D Function, Constant, or an expression) – Initial condition for depth averaged velocity uv_3d (vector valued 3D Function, Constant, or an expression) – Initial condition for horizontal velocity tke (scalar 3D Function, Constant, or an expression) – Initial condition for turbulent kinetic energy field psi (scalar 3D Function, Constant, or an expression) – Initial condition for turbulence generic lenght scale field
callbacks = None

CallbackManager object that stores all callbacks

compute_dt_2d(u_scale)[source]

Computes maximum explicit time step from CFL condition.

$\Delta t = \frac{\Delta x}{U}$

Assumes velocity scale $$U = \sqrt{g H} + U_{scale}$$ where $$U_{scale}$$ is estimated advective velocity.

Parameters: u_scale (float or Constant) – User provided maximum advective velocity scale
compute_dt_diffusion(nu_scale)[source]

Computes maximum explicit time step for horizontal diffusion.

$\Delta t = \alpha \frac{(\Delta x)^2}{\nu_{scale}}$

where $$\nu_{scale}$$ is estimated diffusivity scale.

compute_dt_h_advection(u_scale)[source]

Computes maximum explicit time step for horizontal advection

$\Delta t = \frac{\Delta x}{U_{scale}}$

where $$U_{scale}$$ is estimated horizontal advective velocity.

Parameters: u_scale (float or Constant) – User provided maximum horizontal velocity scale
compute_dt_v_advection(w_scale)[source]

Computes maximum explicit time step for vertical advection

$\Delta t = \frac{\Delta z}{W_{scale}}$

where $$W_{scale}$$ is estimated vertical advective velocity.

Parameters: w_scale (float or Constant) – User provided maximum vertical velocity scale
compute_dx_factor()[source]

Computes normalized distance between nodes in the horizontal direction

The factor depends on the finite element space and its polynomial degree. It is used to compute maximal stable time steps.

compute_dz_factor()[source]

Computes a normalized distance between nodes in the vertical direction

The factor depends on the finite element space and its polynomial degree. It is used to compute maximal stable time steps.

compute_mesh_stats()[source]

Computes number of elements, nodes etc and prints to sdtout

create_equations()[source]

Creates all dynamic equations and time integrators

create_fields()[source]

Creates all fields

create_function_spaces()[source]

Creates function spaces

Function spaces are accessible via function_spaces object.

dt = None

Time step

dt_2d = None

Time of the 2D solver

export()[source]

Export all fields to disk

Also evaluates all callbacks set to ‘export’ interval.

export_initial_state = None

Do export initial state. False if continuing a simulation

fields = None

FieldDict that holds all functions needed by the solver object

function_spaces = None

AttrDict that holds all function spaces needed by the solver object

iterate(update_forcings=None, update_forcings3d=None, export_func=None)[source]

Runs the simulation

Iterates over the time loop until time options.simulation_end_time is reached. Exports fields to disk on options.simulation_export_time intervals.

Parameters: update_forcings – User-defined function that takes simulation time as an argument and updates time-dependent boundary conditions of the 2D system (if any). update_forcings_3d – User-defined function that takes simulation time as an argument and updates time-dependent boundary conditions of the 3D equations (if any). export_func – User-defined function (with no arguments) that will be called on every export.
load_state(i_export, outputdir=None, t=None, iteration=None)[source]

Loads simulation state from hdf5 outputs.

This replaces assign_initial_conditions() in model initilization.

This assumes that model setup is kept the same (e.g. time step) and all pronostic state variables are exported in hdf5 format. The required state variables are: elev_2d, uv_2d, uv_3d, salt_3d, temp_3d, tke_3d, psi_3d

Currently hdf5 field import only works for the same number of MPI processes.

Parameters: i_export (int) – export index to load outputdir (string) – (optional) directory where files are read from. By default options.output_directory. t (float) – simulation time. Overrides the time stamp stored in the hdf5 files. iteration (int) – Overrides the iteration count in the hdf5 files.
mesh = None

3D :classMesh

mesh2d = None

2D :classMesh

options = None

Dictionary of all options. A ModelOptions3d object.

print_state(cputime)[source]

Print a summary of the model state on stdout

Parameters: cputime (float) – Measured CPU time
set_time_step()[source]

Sets the model the model time step

If the time integrator supports automatic time step, and ModelOptions3d.timestepper_options.use_automatic_timestep is True, we compute the maximum time step allowed by the CFL condition. Otherwise uses ModelOptions3d.timestep.

Once the time step is determined, will adjust it to be an integer fraction of export interval options.simulation_export_time.

thetis.solver.timed_region()

Log.Event(type cls, name, klass=None)

thetis.solver.timed_stage()

Log.Stage(type cls, name)

## thetis.solver2d module¶

Module for 2D depth averaged solver

class thetis.solver2d.FlowSolver2d(mesh2d, bathymetry_2d, options=None)[source]

Main object for 2D depth averaged solver

Example

Create mesh

from thetis import *
mesh2d = RectangleMesh(20, 20, 10e3, 10e3)


Create bathymetry function, set a constant value

fs_p1 = FunctionSpace(mesh2d, 'CG', 1)
bathymetry_2d = Function(fs_p1, name='Bathymetry').assign(10.0)


Create solver object and set some options

solver_obj = solver2d.FlowSolver2d(mesh2d, bathymetry_2d)
options = solver_obj.options
options.element_family = 'dg-dg'
options.polynomial_degree = 1
options.timestepper_type = 'CrankNicolson'
options.simulation_export_time = 50.0
options.simulation_end_time = 3600.
options.timestep = 25.0


solver_obj.create_function_spaces()
init_elev = Function(solver_obj.function_spaces.H_2d)
coords = SpatialCoordinate(mesh2d)
init_elev.project(exp(-((coords[0] - 4e3)**2 + (coords[1] - 4.5e3)**2)/2.2e3**2))
solver_obj.assign_initial_conditions(elev=init_elev)


Run simulation

solver_obj.iterate()


See the manual for more complex examples.

Parameters: mesh2d – Mesh object of the 2D mesh bathymetry_2d (Function) – Bathymetry of the domain. Bathymetry stands for the mean water depth (positive downwards). options (ModelOptions2d instance) – Model options (optional). Model options can also be changed directly via the options class property.
add_callback(callback, eval_interval='export')[source]

Parameters: callback – DiagnosticCallback instance eval_interval (string) – Determines when callback will be evaluated, either ‘export’ or ‘timestep’ for evaluating after each export or time step.
assign_initial_conditions(elev=None, uv=None, tracer=None)[source]

Assigns initial conditions

Parameters: elev (scalar Function, Constant, or an expression) – Initial condition for water elevation uv (vector valued Function, Constant, or an expression) – Initial condition for depth averaged velocity
callbacks = None

CallbackManager object that stores all callbacks

compute_time_step(u_scale=Constant(FiniteElement('Real', None, 0), 16))[source]

Computes maximum explicit time step from CFL condition.

$\Delta t = \frac{\Delta x}{U}$

Assumes velocity scale $$U = \sqrt{g H} + U_{scale}$$ where $$U_{scale}$$ is estimated advective velocity.

Parameters: u_scale (float or Constant) – User provided maximum advective velocity scale
create_equations()[source]

Creates shallow water equations

create_exporters()[source]

Creates file exporters

create_function_spaces()[source]

Creates function spaces

Function spaces are accessible via function_spaces object.

create_timestepper()[source]

Creates time stepper instance

dt = None

Time step

export()[source]

Export all fields to disk

Also evaluates all callbacks set to ‘export’ interval.

export_initial_state = None

Do export initial state. False if continuing a simulation

fields = None

FieldDict that holds all functions needed by the solver object

function_spaces = None

AttrDict that holds all function spaces needed by the solver object

initialize()[source]

Creates function spaces, equations, time stepper and exporters

iterate(update_forcings=None, export_func=None)[source]

Runs the simulation

Iterates over the time loop until time options.simulation_end_time is reached. Exports fields to disk on options.simulation_export_time intervals.

Parameters: update_forcings – User-defined function that takes simulation time as an argument and updates time-dependent boundary conditions (if any). export_func – User-defined function (with no arguments) that will be called on every export.
load_state(i_export, outputdir=None, t=None, iteration=None)[source]

Loads simulation state from hdf5 outputs.

This replaces assign_initial_conditions() in model initilization.

This assumes that model setup is kept the same (e.g. time step) and all pronostic state variables are exported in hdf5 format. The required state variables are: elev_2d, uv_2d

Currently hdf5 field import only works for the same number of MPI processes.

Parameters: i_export (int) – export index to load outputdir (string) – (optional) directory where files are read from. By default options.output_directory. t (float) – simulation time. Overrides the time stamp stored in the hdf5 files. iteration (int) – Overrides the iteration count in the hdf5 files.
options = None

Dictionary of all options. A ModelOptions2d object.

print_state(cputime)[source]

Print a summary of the model state on stdout

Parameters: cputime (float) – Measured CPU time
set_time_step(alpha=0.05)[source]

Sets the model the model time step

If the time integrator supports automatic time step, and ModelOptions2d.timestepper_options.use_automatic_timestep is True, we compute the maximum time step allowed by the CFL condition. Otherwise uses ModelOptions2d.timestep.

Parameters: alpha (float) – CFL number scaling factor
thetis.solver2d.timed_region()

Log.Event(type cls, name, klass=None)

thetis.solver2d.timed_stage()

Log.Stage(type cls, name)

## thetis.stability_functions module¶

Implements turbulence closure model stability functions.

$\begin{split}S_m &= S_m(\alpha_M, \alpha_N) \\ S_\rho &= S_\rho(\alpha_M, \alpha_N)\end{split}$

where $$\alpha_M, \alpha_N$$ are the normalized shear and buoyancy frequency

$\begin{split}\alpha_M &= \frac{k^2}{\varepsilon^2} M^2 \\ \alpha_N &= \frac{k^2}{\varepsilon^2} N^2\end{split}$

The following stability functions have been implemented

• Canuto A
• Canuto B
• Kantha-Clayson
• Cheng

References:

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

Burchard, H. and Bolding, K. (2001). Comparative Analysis of Four Second-Moment Turbulence Closure Models for the Oceanic Mixed Layer. Journal of Physical Oceanography, 31(8):1943–1968. http://dx.doi.org/10.1175/1520-0485(2001)031

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

class thetis.stability_functions.StabilityFunction(lim_alpha_shear=True, lim_alpha_buoy=True, smooth_alpha_buoy_lim=True, alpha_buoy_crit=-1.2)[source]

Bases: object

Base class for all stability functions

Parameters: lim_alpha_shear (bool) – limit maximum $$\alpha_M$$ values (see Umlauf and Burchard (2005) eq. 44) lim_alpha_buoy (bool) – limit minimum (negative) $$\alpha_N$$ values (see Umlauf and Burchard (2005)) smooth_alpha_buoy_lim (bool) – if $$\alpha_N$$ is limited, apply a smooth limiter (see Burchard and Bolding (2001) eq. 19). Otherwise $$\alpha_N$$ is clipped at minimum value. alpha_buoy_crit (float) – parameter for $$\alpha_N$$ smooth limiter
compute_c3_minus(c1, c2, ri_st)[source]

Compute c3_minus parameter from c1, c2 and stability functions.

c3_minus is solved from equation

$Ri_{st} = \frac{s_m}{s_h} \frac{c2 - c1}{c2 - c3_minus}$

where $$Ri_{st}$$ is the steady state gradient Richardson number. (see Burchard and Bolding, 2001, eq 32)

compute_cmu0()[source]

Computes the paramenter c_mu_0 from stability function parameters

Umlauf and Buchard (2005) eq A.22

compute_kappa(sigma_psi, n, c1, c2)[source]

Computes von Karman constant from the Psi Schmidt number.

n, c1, c2 are GLS model parameters.

from Umlauf and Burchard (2003) eq (14)

eval_funcs(alpha_buoy, alpha_shear)[source]

Evaluate (unlimited) stability functions

from Burchard and Petersen (1999) eqns (30) and (31)

Parameters: alpha_buoy – normalized buoyancy frequency $$\alpha_N$$ alpha_shear – normalized shear frequency $$\alpha_M$$
evaluate(shear2, buoy2, k, eps)[source]

Evaluates stability functions. Applies limiters on alpha_buoy and alpha_shear.

Parameters: shear2 – $$M^2$$ buoy2 – $$N^2$$ k – turbulent kinetic energy eps – TKE dissipation rate
get_alpha_buoy_min()[source]

Compute minimum alpha buoy

from Umlauf and Buchard (2005) table 3

get_alpha_buoy_smooth_min(alpha_buoy)[source]

Compute smoothed alpha_buoy minimum

from Burchard and Petersen (1999) eq (19)

Parameters: alpha_buoy – normalized buoyancy frequency $$\alpha_N$$
get_alpha_shear_max(alpha_buoy, alpha_shear)[source]

Compute maximum alpha shear

from Umlauf and Buchard (2005) eq (44)

Parameters: alpha_buoy – normalized buoyancy frequency $$\alpha_N$$ alpha_shear – normalized shear frequency $$\alpha_M$$
class thetis.stability_functions.StabilityFunctionCanutoA(lim_alpha_shear=True, lim_alpha_buoy=True, smooth_alpha_buoy_lim=True, alpha_buoy_crit=-1.2)[source]

Canuto et al. (2001) version A stability functions

Parameters are from Umlauf and Buchard (2005), Table 1

class thetis.stability_functions.StabilityFunctionCanutoB(lim_alpha_shear=True, lim_alpha_buoy=True, smooth_alpha_buoy_lim=True, alpha_buoy_crit=-1.2)[source]

Canuto et al. (2001) version B stability functions

Parameters are from Umlauf and Buchard (2005), Table 1

class thetis.stability_functions.StabilityFunctionKanthaClayson(lim_alpha_shear=True, lim_alpha_buoy=True, smooth_alpha_buoy_lim=True, alpha_buoy_crit=-1.2)[source]

Kantha and Clayson (1994) quasi-equilibrium stability functions

Parameters are from Umlauf and Buchard (2005), Table 1

class thetis.stability_functions.StabilityFunctionCheng(lim_alpha_shear=True, lim_alpha_buoy=True, smooth_alpha_buoy_lim=True, alpha_buoy_crit=-1.2)[source]

Cheng et al. (2002) quasi-equilibrium stability functions

Parameters are from Umlauf and Buchard (2005), Table 1

thetis.stability_functions.compute_normalized_frequencies(shear2, buoy2, k, eps)[source]

Computes normalized buoyancy and shear frequency squared.

$\begin{split}\alpha_M &= \frac{k^2}{\varepsilon^2} M^2 \\ \alpha_N &= \frac{k^2}{\varepsilon^2} N^2\end{split}$

From Burchard and Bolding (2001).

Parameters: shear2 – $$M^2$$ buoy2 – $$N^2$$ k – turbulent kinetic energy eps – TKE dissipation rate

## thetis.timeintegrator module¶

Generic time integration schemes to advance equations in time.

class thetis.timeintegrator.CrankNicolson(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, theta=0.5, semi_implicit=False)[source]

Standard Crank-Nicolson time integration scheme.

Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options theta (float) – Implicitness parameter, default 0.5 semi_implicit (bool) – If True use a linearized semi-implicit scheme
advance(t, update_forcings=None)[source]

Advances equations for one time step.

cfl_coeff = inf
initialize(solution)[source]

update_solver()[source]

Create solver objects

class thetis.timeintegrator.ForwardEuler(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={})[source]

Standard forward Euler time integration scheme.

Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options
advance(t, update_forcings=None)[source]

Advances equations for one time step.

cfl_coeff = 1.0
initialize(solution)[source]

update_solver()[source]
class thetis.timeintegrator.LeapFrogAM3(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]

Leap-Frog Adams-Moulton 3 ALE time integrator

Defined in (2.27)-(2.30) in [1]; (2.21)-(2.22) in [2]

[1] Shchepetkin and McWilliams (2005). The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean Modelling, 9(4):347-404. http://dx.doi.org/10.1016/j.ocemod.2013.04.010

[2] Shchepetkin and McWilliams (2009). Computational Kernel Algorithms for Fine-Scale, Multiprocess, Longtime Oceanic Simulations, 14:121-183. http://dx.doi.org/10.1016/S1570-8659(08)01202-0

Parameters: equation (Equation object) – equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
advance(t, update_forcings=None)[source]

Advances equations for one time step.

cfl_coeff = 1.5874
correct()[source]

Correction step from $$t_{n}$$ to $$t_{n+1}$$

Let $$M_n$$ denote the mass matrix at time $$t_{n}$$. The correction step is

$M_{n+1} T_{n+1} = M_{n} T_{n} + \Delta t L_{n+1/2}$

This is Euler ALE step: LHS is evaluated in $$\Omega_{n+1}$$, RHS in $$\Omega_n$$.

eval_rhs()[source]
initialize(solution)[source]

predict()[source]

Prediction step from $$t_{n-1/2}$$ to $$t_{n+1/2}$$

Let $$M_n$$ denote the mass matrix at time $$t_{n}$$. The prediction step is

$\begin{split}T_{n-1/2} &= (1/2 - 2\gamma) T_{n-1} + (1/2 + 2 \gamma) T_{n} \\ M_n T_{n+1/2} &= M_n T_{n-1/2} + \Delta t (1 - 2\gamma) M_n L_{n}\end{split}$

This is computed in a fixed mesh: all terms are evaluated in $$\Omega_n$$.

class thetis.timeintegrator.PressureProjectionPicard(equation, equation_mom, solution, fields, dt, bnd_conditions=None, solver_parameters={}, solver_parameters_mom={}, theta=0.5, semi_implicit=False, iterations=2)[source]

Pressure projection scheme with Picard iteration for shallow water equations

Parameters: equation (Equation object) – free surface equation equation_mom (Equation object) – momentum equation solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options solver_parameters_mom (dict) – PETSc solver options for velocity solver theta (float) – Implicitness parameter, default 0.5 semi_implicit (bool) – If True use a linearized semi-implicit scheme iterations (int) – Number of Picard iterations
advance(t, updateForcings=None)[source]

Advances equations for one time step.

cfl_coeff = 1.0
initialize(solution)[source]

update_solver()[source]

Create solver objects

class thetis.timeintegrator.SSPRK22ALE(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={}, terms_to_add='all')[source]

SSPRK(2,2) ALE time integrator for 3D fields

The scheme is

$\begin{split}u^{(1)} &= u^{n} + \Delta t F(u^{n}) \\ u^{n+1} &= u^{n} + \frac{\Delta t}{2}(F(u^{n}) + F(u^{(1)}))\end{split}$

Both stages are implemented as ALE updates from geometry $$\Omega_n$$ to $$\Omega_{(1)}$$, and $$\Omega_{n+1}$$.

Parameters: equation (Equation object) – equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options terms_to_add ('all' or list of 'implicit', 'explicit', 'source'.) – Defines which terms of the equation are to be added to this solver. Default ‘all’ implies [‘implicit’, ‘explicit’, ‘source’].
advance(t, update_forcings=None)[source]

Advances equations for one time step.

cfl_coeff = 1.0
initialize(solution)[source]

prepare_stage(i_stage, t, update_forcings=None)[source]

Preprocess stage i_stage.

This must be called prior to updating mesh geometry.

solve_stage(i_stage)[source]

Solves i-th stage

stage_one_prep()[source]

Preprocess first stage: compute all forms on the old geometry

stage_one_solve()[source]

First stage: solve $$u^{(1)}$$ given previous solution $$u^n$$.

This is a forward Euler ALE step between domains $$\Omega^n$$ and $$\Omega^{(1)}$$:

$\int_{\Omega^{(1)}} u^{(1)} \psi dx = \int_{\Omega^n} u^n \psi dx + \Delta t \int_{\Omega^n} F(u^n) \psi dx$
stage_two_prep()[source]

Preprocess 2nd stage: compute all forms on the old geometry

stage_two_solve()[source]

2nd stage: solve $$u^{n+1}$$ given previous solutions $$u^n, u^{(1)}$$.

This is an ALE step:

$\begin{split}\int_{\Omega^{n+1}} u^{n+1} \psi dx &= \int_{\Omega^n} u^n \psi dx \\ &+ \frac{\Delta t}{2} \int_{\Omega^n} F(u^n) \psi dx \\ &+ \frac{\Delta t}{2} \int_{\Omega^{(1)}} F(u^{(1)}) \psi dx\end{split}$
class thetis.timeintegrator.SteadyState(equation, solution, fields, dt, bnd_conditions=None, solver_parameters={})[source]

Time integrator that solves the steady state equations, leaving out the mass terms

Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds bnd_conditions (dict) – Dictionary of boundary conditions passed to the equation solver_parameters (dict) – PETSc solver options
advance(t, update_forcings=None)[source]

Advances equations for one time step.

cfl_coeff = inf
initialize(solution)[source]

update_solver()[source]

Create solver objects

class thetis.timeintegrator.TimeIntegrator(equation, solution, fields, dt, solver_parameters={})[source]

Base class for all time integrator objects that march a single equation

Parameters: equation (Equation object) – the equation to solve solution – Function where solution will be stored fields (dict of Function or Constant objects) – Dictionary of fields that are passed to the equation dt (float) – time step in seconds solver_parameters (dict) – PETSc solver options
set_dt(dt)[source]

Update time step

class thetis.timeintegrator.TimeIntegratorBase[source]

Bases: object

Abstract class that defines the API for all time integrators

Both TimeIntegrator and CoupledTimeIntegrator inherit from this class.

advance(t, update_forcings=None)[source]

Advances equations for one time step

Parameters: t (float) – simulation time update_forcings – user-defined function that takes the simulation time and updates any time-dependent boundary conditions
initialize(init_solution)[source]

Initialize the time integrator

Parameters: init_solution – initial solution
thetis.timeintegrator.timed_region()

Log.Event(type cls, name, klass=None)

thetis.timeintegrator.timed_stage()

Log.Stage(type cls, name)

## thetis.timezone module¶

Timezone definitions and conversion methods

class thetis.timezone.FixedTimeZone(offset, name)[source]

Bases: pytz._FixedOffset

Class that represents a fixed time zone defined by UTC offset in hours.

arg int offset: timezone UTC offset in hours arg str name: timezone name

tzname(dt)[source]

datetime -> string name of time zone.

thetis.timezone.datetime_to_epoch(t)[source]

Convert python datetime object to epoch time stamp.

thetis.timezone.epoch_to_datetime(t)[source]

Convert python datetime object to epoch time stamp.

## thetis.tracer_eq module¶

3D advection diffusion equation for tracers.

The advection-diffusion equation of tracer $$T$$ in conservative form reads

(13)$\frac{\partial T}{\partial t} + \nabla_h \cdot (\textbf{u} T) + \frac{\partial (w T)}{\partial z} = \nabla_h \cdot (\mu_h \nabla_h T) + \frac{\partial}{\partial z} \Big(\mu \frac{\partial T}{\partial z}\Big)$

where $$\nabla_h$$ denotes horizontal gradient, $$\textbf{u}$$ and $$w$$ are the horizontal and vertical velocities, respectively, and $$\mu_h$$ and $$\mu$$ denote horizontal and vertical diffusivity.

class thetis.tracer_eq.TracerEquation(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_symmetric_surf_bnd=True, use_lax_friedrichs=True)[source]

3D tracer advection-diffusion equation (13) in conservative form

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_symmetric_surf_bnd (bool) – If True, use symmetric surface boundary condition in the horizontal advection term
class thetis.tracer_eq.TracerTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_symmetric_surf_bnd=True, use_lax_friedrichs=True)[source]

Generic tracer term that provides commonly used members and mapping for boundary functions.

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_symmetric_surf_bnd (bool) – If True, use symmetric surface boundary condition in the horizontal advection term
get_bnd_functions(c_in, uv_in, elev_in, bnd_id, bnd_conditions)[source]

Returns external values of tracer and uv for all supported boundary conditions.

Volume flux (flux) and normal velocity (un) are defined positive out of the domain.

Parameters: c_in – Internal value of tracer uv_in – Internal value of horizontal velocity elev_in – Internal value of elevation bnd_id (int) – boundary id bnd_conditions – dict of boundary conditions: {bnd_id: {field: value, ...}, ...}
class thetis.tracer_eq.HorizontalAdvectionTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_symmetric_surf_bnd=True, use_lax_friedrichs=True)[source]

Horizontal advection term $$\nabla_h \cdot (\textbf{u} T)$$

$\int_\Omega \nabla_h \cdot (\textbf{u} T) \phi dx = -\int_\Omega T\textbf{u} \cdot \nabla_h \phi dx + \int_{\mathcal{I}_h\cup\mathcal{I}_v} T^{\text{up}} \text{avg}(\textbf{u}) \cdot \text{jump}(\phi \textbf{n}_h) dS$

where the right hand side has been integrated by parts; $$\mathcal{I}_h,\mathcal{I}_v$$ denote the set of horizontal and vertical facets, $$\textbf{n}_h$$ is the horizontal projection of the unit normal vector, $$T^{\text{up}}$$ is the upwind value, and $$\text{jump}$$ and $$\text{avg}$$ denote the jump and average operators across the interface.

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_symmetric_surf_bnd (bool) – If True, use symmetric surface boundary condition in the horizontal advection term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.tracer_eq.VerticalAdvectionTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_symmetric_surf_bnd=True, use_lax_friedrichs=True)[source]

Vertical advection term $$\partial (w T)/(\partial z)$$

$\int_\Omega \frac{\partial (w T)}{\partial z} \phi dx = - \int_\Omega T w \frac{\partial \phi}{\partial z} dx + \int_{\mathcal{I}_v} T^{\text{up}} \text{avg}(w) \text{jump}(\phi n_z) dS$

where the right hand side has been integrated by parts; $$\mathcal{I}_v$$ denotes the set of vertical facets, $$n_z$$ is the vertical projection of the unit normal vector, $$T^{\text{up}}$$ is the upwind value, and $$\text{jump}$$ and $$\text{avg}$$ denote the jump and average operators across the interface.

In the case of ALE moving mesh we substitute $$w$$ with $$w - w_m$$, $$w_m$$ being the mesh velocity.

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_symmetric_surf_bnd (bool) – If True, use symmetric surface boundary condition in the horizontal advection term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.tracer_eq.HorizontalDiffusionTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_symmetric_surf_bnd=True, use_lax_friedrichs=True)[source]

Horizontal diffusion term $$-\nabla_h \cdot (\mu_h \nabla_h T)$$

Using the symmetric interior penalty method the weak form becomes

$\begin{split}\int_\Omega \nabla_h \cdot (\mu_h \nabla_h T) \phi dx =& -\int_\Omega \mu_h (\nabla_h \phi) \cdot (\nabla_h T) dx \\ &+ \int_{\mathcal{I}_h\cup\mathcal{I}_v} \text{jump}(\phi \textbf{n}_h) \cdot \text{avg}(\mu_h \nabla_h T) dS + \int_{\mathcal{I}_h\cup\mathcal{I}_v} \text{jump}(T \textbf{n}_h) \cdot \text{avg}(\mu_h \nabla \phi) dS \\ &- \int_{\mathcal{I}_h\cup\mathcal{I}_v} \sigma \text{avg}(\mu_h) \text{jump}(T \textbf{n}_h) \cdot \text{jump}(\phi \textbf{n}_h) dS\end{split}$

where $$\sigma$$ is a penalty parameter, see Epshteyn and Riviere (2007).

Epshteyn and Riviere (2007). Estimation of penalty parameters for symmetric interior penalty Galerkin methods. Journal of Computational and Applied Mathematics, 206(2):843-872. http://dx.doi.org/10.1016/j.cam.2006.08.029

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_symmetric_surf_bnd (bool) – If True, use symmetric surface boundary condition in the horizontal advection term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.tracer_eq.VerticalDiffusionTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_symmetric_surf_bnd=True, use_lax_friedrichs=True)[source]

Vertical diffusion term $$-\frac{\partial}{\partial z} \Big(\mu \frac{T}{\partial z}\Big)$$

Using the symmetric interior penalty method the weak form becomes

$\begin{split}\int_\Omega \frac{\partial}{\partial z} \Big(\mu \frac{T}{\partial z}\Big) \phi dx =& -\int_\Omega \mu \frac{\partial T}{\partial z} \frac{\partial \phi}{\partial z} dz \\ &+ \int_{\mathcal{I}_{h}} \text{jump}(\phi n_z) \text{avg}\Big(\mu \frac{\partial T}{\partial z}\Big) dS + \int_{\mathcal{I}_{h}} \text{jump}(T n_z) \text{avg}\Big(\mu \frac{\partial \phi}{\partial z}\Big) dS \\ &- \int_{\mathcal{I}_{h}} \sigma \text{avg}(\mu) \text{jump}(T n_z) \cdot \text{jump}(\phi n_z) dS\end{split}$

where $$\sigma$$ is a penalty parameter, see Epshteyn and Riviere (2007).

Epshteyn and Riviere (2007). Estimation of penalty parameters for symmetric interior penalty Galerkin methods. Journal of Computational and Applied Mathematics, 206(2):843-872. http://dx.doi.org/10.1016/j.cam.2006.08.029

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_symmetric_surf_bnd (bool) – If True, use symmetric surface boundary condition in the horizontal advection term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.tracer_eq.SourceTerm(function_space, bathymetry=None, v_elem_size=None, h_elem_size=None, use_symmetric_surf_bnd=True, use_lax_friedrichs=True)[source]

Generic source term

$F_s = \int_\Omega \sigma \phi dx$

where $$\sigma$$ is a user defined scalar Function.

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size use_symmetric_surf_bnd (bool) – If True, use symmetric surface boundary condition in the horizontal advection term
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.

## thetis.tracer_eq_2d module¶

2D advection diffusion equation for tracers.

The advection-diffusion equation of tracer $$T$$ in non-conservative form reads

(14)$\frac{\partial T}{\partial t} + \nabla_h \cdot (\textbf{u} T) = \nabla_h \cdot (\mu_h \nabla_h T)$

where $$\nabla_h$$ denotes horizontal gradient, $$\textbf{u}$$ are the horizontal velocities, and $$\mu_h$$ denotes horizontal diffusivity.

class thetis.tracer_eq_2d.TracerEquation2D(function_space, bathymetry=None, use_lax_friedrichs=False)[source]

2D tracer advection-diffusion equation (13) in conservative form

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (2D Function or Constant) – bathymetry of the domain use_symmetric_surf_bnd (bool) – If True, use symmetric surface boundary condition in the horizontal advection term
class thetis.tracer_eq_2d.TracerTerm(function_space, bathymetry=None, use_lax_friedrichs=True)[source]

Generic tracer term that provides commonly used members and mapping for boundary functions.

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (2D Function or Constant) – bathymetry of the domain
get_bnd_functions(c_in, uv_in, elev_in, bnd_id, bnd_conditions)[source]

Returns external values of tracer and uv for all supported boundary conditions.

Volume flux (flux) and normal velocity (un) are defined positive out of the domain.

Parameters: c_in – Internal value of tracer uv_in – Internal value of horizontal velocity elev_in – Internal value of elevation bnd_id (int) – boundary id bnd_conditions – dict of boundary conditions: {bnd_id: {field: value, ...}, ...}
class thetis.tracer_eq_2d.HorizontalAdvectionTerm(function_space, bathymetry=None, use_lax_friedrichs=True)[source]

Advection of tracer term, $$\bar{\textbf{u}} \cdot \nabla T$$

The weak form is

$\int_\Omega \bar{\textbf{u}} \cdot \boldsymbol{\psi} \cdot \nabla T dx = - \int_\Omega \nabla_h \cdot (\bar{\textbf{u}} \boldsymbol{\psi}) \cdot T dx + \int_\Gamma \text{avg}(T) \cdot \text{jump}(\boldsymbol{\psi} (\bar{\textbf{u}}\cdot\textbf{n})) dS$

where the right hand side has been integrated by parts; $$\textbf{n}$$ is the unit normal of the element interfaces, and $$\text{jump}$$ and $$\text{avg}$$ denote the jump and average operators across the interface.

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (2D Function or Constant) – bathymetry of the domain
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.tracer_eq_2d.HorizontalDiffusionTerm(function_space, bathymetry=None, use_lax_friedrichs=True)[source]

Horizontal diffusion term $$-\nabla_h \cdot (\mu_h \nabla_h T)$$

Using the symmetric interior penalty method the weak form becomes

$\begin{split}-\int_\Omega \nabla_h \cdot (\mu_h \nabla_h T) \phi dx =& \int_\Omega \mu_h (\nabla_h \phi) \cdot (\nabla_h T) dx \\ &- \int_{\mathcal{I}_h\cup\mathcal{I}_v} \text{jump}(\phi \textbf{n}_h) \cdot \text{avg}(\mu_h \nabla_h T) dS - \int_{\mathcal{I}_h\cup\mathcal{I}_v} \text{jump}(T \textbf{n}_h) \cdot \text{avg}(\mu_h \nabla \phi) dS \\ &+ \int_{\mathcal{I}_h\cup\mathcal{I}_v} \sigma \text{avg}(\mu_h) \text{jump}(T \textbf{n}_h) \cdot \text{jump}(\phi \textbf{n}_h) dS\end{split}$

where $$\sigma$$ is a penalty parameter, see Epshteyn and Riviere (2007).

Epshteyn and Riviere (2007). Estimation of penalty parameters for symmetric interior penalty Galerkin methods. Journal of Computational and Applied Mathematics, 206(2):843-872. http://dx.doi.org/10.1016/j.cam.2006.08.029

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (2D Function or Constant) – bathymetry of the domain
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.tracer_eq_2d.SourceTerm(function_space, bathymetry=None, use_lax_friedrichs=True)[source]

Generic source term

$F_s = \int_\Omega \sigma \phi dx$

where $$\sigma$$ is a user defined scalar Function.

Parameters: function_space – FunctionSpace where the solution belongs bathymetry (2D Function or Constant) – bathymetry of the domain
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.

## thetis.turbines module¶

Classes related to tidal turbine farms in Thetis.

class thetis.turbines.TurbineFarm(farm_options, subdomain_id, u, v, dt)[source]

Bases: object

Evaluates power output, costs and profit of tidal turbine farm

Cost is simply the number of turbines evaluated as the integral of the turbine density.

Power output is calculated as the amount of energy taken out of the flow by the turbine drag term. This is in general an over-estimate of the ‘usefully extractable’ energy. Furthermore no density is included, so that assuming a density of 1000 kg/m^3 and all further quantities in SI, the power is measured in kW.

Profit is calculated as:

Profit = Power - break_even_wattage * Cost

With the above assumptions, break_even_wattage should be specified in kW and can be interpreted as the average power per turbine required to ‘break even’, i.e. Profit=0.

Power output and profit are time-integrated (simple first order) and time-averages are available as average_power and average_profit.

Parameters: farm_options – a TidalTurbineFarmOptions object that define the farm and the turbines used subdomain_id (int) – the farm is restricted to this subdomain u,v – the depth-averaged velocity field dt (float) – used for time-integration.
evaluate_timestep()[source]

Perform time integration and return current power and time-averaged power and profit.

class thetis.turbines.TurbineFunctionalCallback(solver_obj, **kwargs)[source]

DiagnosticCallback that evaluates the performance of each tidal turbine farm.

Parameters: solver_obj – a FlowSolver2d object containing the tidal_turbine_farms kwargs – see DiagnosticCallback
average_power

The sum of the time-averaged power output of all farms

average_profit

The sum of the time-averaged profit output of all farms

farms = None

The sum of the number of turbines in all farms

integrated_power

The sum of the time-integrated power output of all farms

message_str(current_power, average_power, average_profit)[source]

A string representation.

Parameters: args – If provided, these will be the return value from __call__().
name = 'turbine'
variable_names = ['current_power', 'average_power', 'average_profit']
class thetis.turbines.TurbineOptimisationCallback(solver_obj, turbine_functional_callback, **kwargs)[source]

DiagnosticOptimisationCallback that evaluates the performance of each tidal turbine farm during an optimisation.

See the optimisation module for more info about the use of OptimisationCallbacks.

Parameters: solver_obj – a FlowSolver2d object turbine_functional_callback – a TurbineFunctionalCallback used in the forward model
compute_values(*args)[source]

Compute diagnostic values.

This method is to be implemented in concrete subclasses of a DiagnosticOptimisationCallback. The number of arguments varies depending on which of the 6 [eval|derivative|hessian]_cb_[pre|post] callbacks this is used as. The last argument always contains the current controls. In the “pre” callbacks this is the only argument. In all “post” callbacks the 0th argument is the current functional value. eval_cb_post is given two arguments: functional and controls. derivative_cb_post and hessian_cb_post are given three arguments with args[1] being the derivative/hessian just calculated.

message_str(cost, average_power, average_profit)[source]

A string representation.

Parameters: args – If provided, these will be the return value from __call__().
name = 'farm_optimisation'
variable_names = ['cost', 'average_power', 'average_profit']

## thetis.turbulence module¶

### Generic Length Scale Turbulence Closure model¶

This model solves two dynamic equations, for turbulent kinetic energy (TKE, $$k$$) and one for an additional variable, the generic length scale $$\psi$$ [1]:

(15)$\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$
(16)$\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})$

with the production $$P$$ and buoyancy production $$B$$ are

$\begin{split}P &= \nu M^2 \\ B &= -\mu N^2\end{split}$

where $$M$$ and $$N$$ are the shear and buoyancy frequencies

$\begin{split}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}\end{split}$

The generic lenght scale variable is defined as

$\psi = (c_\mu^0)^p k^m l^n$

where $$p, m, n$$ are parameters and $$c_\mu^0$$ is an empirical constant.

The parameters $$c_1,c_2,c_3,f_{wall}$$ depend on the chosen closure. The parameter $$c_3$$ takes two values: $$c_3^-$$ in stably stratified regime, and $$c_3^+$$ in unstably stratified cases.

Turbulent length scale $$l$$, and the TKE dissipation rate $$\varepsilon$$ are obtained diagnostically as

$\begin{split}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}\end{split}$

Finally the vertical eddy viscosity and diffusivity are also computed diagnostically

$\begin{split}\nu &= \sqrt{2k} l S_m \\ \mu &= \sqrt{2k} l S_\rho\end{split}$

Stability functions $$S_m$$ and $$S_\rho$$ are defined in [2] or [3]. Implementation follows [4].

#### Supported closures¶

The GLS model parameters are controlled via the GLSModelOptions class.

The parameters can be accessed from the solver object:

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:

• $$k-\varepsilon$$ model
This closure is obtained with $$p=3, m=3/2, n=-1$$, resulting in $$\psi=\varepsilon$$. To use this option set closure_name = k-epsilon
• $$k-\omega$$ model
This closure is obtained with $$p=-1, m=1/2, n=-1$$, resulting in $$\psi=\omega$$. To use this option set closure_name = k-omega
• GLS model A
This closure is obtained with $$p=2, m=1, n=-2/3$$, resulting in $$\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 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
class thetis.turbulence.BuoyFrequencySolver(rho, n2, n2_tmp, relaxation=1.0, minval=1e-12)[source]

Bases: object

Computes buoyancy frequency squared form the given horizontal velocity field.

$N^2 = -\frac{g}{\rho_0}\frac{\partial \rho}{\partial z}$
Parameters: rho (Function) – water density field n2 (Function) – $$N^2$$ field n2_tmp (Function) – temporary field relaxation (float) – relaxation coefficient for mixing old and new values N2 = relaxation*N2_new + (1-relaxation)*N2_old minval (float) – minimum value for $$N^2$$
solve(init_solve=False)[source]

Computes buoyancy frequency

Parameters: init_solve (bool) – Set to True if solving for the first time, skips relaxation
class thetis.turbulence.GLSVerticalDiffusionTerm(function_space, schmidt_nb, bathymetry=None, v_elem_size=None, h_elem_size=None)[source]

Vertical diffusion term where the diffusivity is replaced by viscosity/Schmidt number.

Parameters: function_space – FunctionSpace where the solution belongs schmidt_nb – the Schmidt number of TKE or Psi bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.turbulence.GenericLengthScaleModel(solver, k_field, psi_field, uv_field, rho_field, l_field, epsilon_field, eddy_diffusivity, eddy_viscosity, n2, m2, options=None)[source]

Generic Length Scale turbulence closure model implementation

Parameters: solver – FlowSolver object k_field (Function) – turbulent kinetic energy (TKE) field psi_field (Function) – generic length scale field uv_field (Function) – horizontal velocity field rho_field (Function) – water density field l_field (Function) – turbulence length scale field epsilon_field (Function) – TKE dissipation rate field eddy_diffusivity (Function) – eddy diffusivity field eddy_viscosity (Function) – eddy viscosity field n2 (Function) – field for buoyancy frequency squared m2 (Function) – field for vertical shear frequency squared options – GLS model options
initialize()[source]

Initializes fields

postprocess()[source]

Updates all diagnostic variables that depend on the turbulence state variables $$k,\psi$$.

To be called after updating the turbulence PDEs.

preprocess(init_solve=False)[source]

Computes all diagnostic variables that depend on the mean flow model variables.

To be called before updating the turbulence PDEs.

class thetis.turbulence.P1Average(p0, p1, p1dg)[source]

Bases: object

Takes a discontinuous field and computes a P1 field by averaging around nodes

Source must be either a P0 or P1DG Function. The averaging operation is both mass conservative and positivity preserving.

Parameters: p0 – P0 function space p1 – P1 function space p1dg – P1DG function space
apply(source, solution)[source]

Averages discontinuous Function source on P1 Function solution

update_volumes()[source]

Computes nodal volume of the P1 and P1DG function function_spaces

This must be called when the mesh geometry is updated

class thetis.turbulence.PacanowskiPhilanderModel(solver, uv_field, rho_field, eddy_diffusivity, eddy_viscosity, n2, m2, options=None)[source]

Gradient Richardson number based model by Pacanowski and Philander (1981).

Given the gradient Richardson number $$Ri$$ the eddy viscosity and diffusivity are computed diagnostically as

$\begin{split}\nu &= \frac{\nu_{max}}{(1 + \alpha Ri)^n} \\ \mu &= \frac{\nu}{1 + \alpha Ri}\end{split}$

where $$\nu_{max},\alpha,n$$ are constant parameters. In unstably stratified cases where $$Ri<0$$, value $$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

Parameters: solver – FlowSolver object uv_field (Function) – horizontal velocity field rho_field (Function) – water density field eddy_diffusivity (Function) – eddy diffusivity field eddy_viscosity (Function) – eddy viscosity field n2 (Function) – field for buoyancy frequency squared m2 (Function) – field for vertical shear frequency squared options – model options
initialize()[source]

Initializes fields

postprocess()[source]

Updates all diagnostic variables that depend on the turbulence state variables.

To be called after evaluating the equations.

preprocess(init_solve=False)[source]

Computes all diagnostic variables that depend on the mean flow model variables.

To be called before updating the turbulence PDEs.

class thetis.turbulence.PsiEquation(function_space, gls_model, bathymetry=None, v_elem_size=None, h_elem_size=None)[source]

Generic length scale equation (16) without advection terms.

Advection of $$\psi$$ is implemented using the standard tracer equation.

Parameters: function_space – FunctionSpace where the solution belongs gls_model – GenericLengthScaleModel object bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size
class thetis.turbulence.PsiSourceTerm(function_space, gls_model, bathymetry=None, v_elem_size=None, h_elem_size=None)[source]

Production and destruction terms of the Psi equation (16)

$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 $$c_3 B$$ is split in two:

$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 $$B_{source} \ge 0$$ and $$B_{sink} < 0$$.

Also implements Neumann boundary condition at top and bottom [7]

$\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 $$z_b$$ is the distance from boundary, and $$z_0$$ is the roughness length.

Parameters: function_space – FunctionSpace where the solution belongs gls_model – GenericLengthScaleModel object bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.turbulence.ShearFrequencySolver(uv, m2, mu, mv, mu_tmp, relaxation=1.0, minval=1e-12)[source]

Bases: object

Computes vertical shear frequency squared form the given horizontal velocity field.

$M^2 = \left(\frac{\partial u}{\partial z}\right)^2 + \left(\frac{\partial v}{\partial z}\right)^2$
Parameters: uv (Function) – horizontal velocity field m2 (Function) – $$M^2$$ field mu (Function) – field for x component of $$M^2$$ mv (Function) – field for y component of $$M^2$$ mu_tmp (Function) – temporary field relaxation (float) – relaxation coefficient for mixing old and new values M2 = relaxation*M2_new + (1-relaxation)*M2_old minval (float) – minimum value for $$M^2$$
solve(init_solve=False)[source]

Computes buoyancy frequency

Parameters: init_solve (bool) – Set to True if solving for the first time, skips relaxation
class thetis.turbulence.SmoothVerticalGradSolver(source, solution)[source]

Bases: object

Computes vertical gradient in a smooth(er) way.

The source is first interpolated on P0 field. The gradient is computed as $$G = (G_{P0} + G_{P1DG})/2$$.

Parameters: source – A Function or expression to differentiate. solution – A Function where the solution will be stored.
solve()[source]

class thetis.turbulence.TKEEquation(function_space, gls_model, bathymetry=None, v_elem_size=None, h_elem_size=None)[source]

Turbulent kinetic energy equation (15) without advection terms.

Advection of TKE is implemented using the standard tracer equation.

Parameters: function_space – FunctionSpace where the solution belongs gls_model – GenericLengthScaleModel object bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size
class thetis.turbulence.TKESourceTerm(function_space, gls_model, bathymetry=None, v_elem_size=None, h_elem_size=None)[source]

Production and destruction terms of the TKE equation (15)

$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 $$B$$ is split in two:

$F_k = P + B_{source} + \frac{k^{n+1}}{k^n}(B_{sink} - \varepsilon)$

with $$B_{source} \ge 0$$ and $$B_{sink} < 0$$.

Parameters: function_space – FunctionSpace where the solution belongs gls_model – GenericLengthScaleModel object bathymetry (3D Function or Constant) – bathymetry of the domain v_elem_size – scalar Function that defines the vertical element size h_elem_size – scalar Function that defines the horizontal element size
residual(solution, solution_old, fields, fields_old, bnd_conditions=None)[source]

Returns an UFL form of the term.

Parameters: solution – solution Function of the corresponding equation solution_old – a time lagged solution Function fields – a dictionary that provides all the remaining fields that the term depends on. The keys of the dictionary should standard field names in field_metadata fields_old – Time lagged dictionary of fields bnd_conditions – A dictionary describing boundary conditions. E.g. {3: {‘elev_2d’: Constant(1.0)}} replaces elev_2d function by a constant on boundary ID 3.
class thetis.turbulence.TurbulenceModel[source]

Bases: object

Base class for all vertical turbulence models

initialize()[source]

Initialize all turbulence fields

postprocess()[source]

Updates all diagnostic variables that depend on the turbulence state variables.

To be called after updating the turbulence PDEs.

preprocess(init_solve=False)[source]

Computes all diagnostic variables that depend on the mean flow model variables.

To be called before updating the turbulence PDEs.

class thetis.turbulence.VerticalGradSolver(source, solution, solver_parameters=None)[source]

Bases: object

Computes vertical gradient in the weak sense.

Parameters: source – A Function or expression to differentiate. solution – A Function where the solution will be stored. Must be in P0 space. solver_parameters (dict) – PETSc solver options
solve()[source]

thetis.turbulence.set_func_max_val(f, maxval)[source]

Sets a minimum value to a Function

thetis.turbulence.set_func_min_val(f, minval)[source]

Sets a minimum value to a Function

thetis.turbulence.timed_region()

Log.Event(type cls, name, klass=None)

thetis.turbulence.timed_stage()

Log.Stage(type cls, name)

## thetis.utility module¶

Utility functions and classes for 3D hydrostatic ocean model

class thetis.utility.ALEMeshUpdater(solver)[source]

Bases: object

Class that handles vertically moving ALE mesh

Mesh geometry is updated to match the elevation field (solver.fields.elev_2d). First the discontinuous elevation field is projected to continuous space, and this field is used to update the mesh coordinates.

This class stores the reference coordinate field and keeps track of the updated mesh coordinates. It also provides a method for computing the mesh velocity from two adjacent elevation fields.

Parameters: solver – FlowSolver object
compute_mesh_velocity_begin()[source]

Stores the current 2D elevation state as the “old” field

compute_mesh_velocity_finalize(c=1.0)[source]

Computes mesh velocity from the elevation difference

Stores the current 2D elevation state as the “new” field, and computes w_mesh using the given time step factor c.

initialize()[source]

Set values for initial mesh (elevation at rest)

update_elem_height()[source]

update_mesh_coordinates()[source]

Updates 3D mesh coordinates to match current elev_2d field

elev_2d is first projected to continous space

class thetis.utility.AttrDict(*args, **kwargs)[source]

Bases: dict

http://stackoverflow.com/questions/4984647/accessing-dict-keys-like-an-attribute-in-python

class thetis.utility.DensitySolver(salinity, temperature, density, eos_class)[source]

Bases: object

Computes density from salinity and temperature using the equation of state.

Water density is defined as

$\rho = \rho'(T, S, p) + \rho_0$

This method computes the density anomaly $$\rho'$$.

Density is computed point-wise assuming that temperature, salinity and density are in the same function space.

Parameters: salinity (Function) – water salinity field temperature (Function) – water temperature field density (Function) – water density field eos_class (EquationOfState) – equation of state that defines water density
solve()[source]

Compute density

class thetis.utility.DensitySolverWeak(salinity, temperature, density, eos_class)[source]

Bases: object

Computes density from salinity and temperature using the equation of state.

Water density is defined as

$\rho = \rho'(T, S, p) + \rho_0$

This method computes the density anomaly $$\rho'$$.

Density is computed in a weak sense by projecting the analytical expression on the density field.

Parameters: salinity (Function) – water salinity field temperature (Function) – water temperature field density (Function) – water density field eos_class (EquationOfState) – equation of state that defines water density
ensure_positive_salinity()[source]

make sure salinity is not negative

some EOS depend on sqrt(salt).

solve()[source]

Compute density

class thetis.utility.ElementContinuity(horizontal, vertical)

Bases: tuple

A named tuple describing the continuity of an element in the horizontal/vertical direction.

The field value is one of “cg”, “hdiv”, or “dg”.

horizontal

Alias for field number 0

vertical

Alias for field number 1

class thetis.utility.EquationOfState[source]

Bases: object

Base class of all equation of state objects

compute_rho(s, th, p, rho0=0.0)[source]

Compute sea water density.

Parameters: s (float or numpy.array) – Salinity expressed on the Practical Salinity Scale 1978 th (float or numpy.array) – Potential temperature in Celsius, referenced to pressure p_r = 0 dbar. p (float or numpy.array) – Pressure in decibars (1 dbar = 1e4 Pa) rho0 (float) – Optional reference density. If provided computes $$\rho' = \rho(S, Th, p) - \rho_0$$ water density float or numpy.array

All pressures are gauge pressures: they are the absolute pressures minus standard atmosperic pressure 10.1325 dbar.

eval(s, th, p, rho0=0.0)[source]

Compute sea water density.

class thetis.utility.ExpandFunctionTo3d(input_2d, output_3d, elem_height=None)[source]

Bases: object

Copy a 2D field to 3D

Copies a field from 2D mesh to 3D mesh, assigning the same value over the vertical dimension. Horizontal function spaces must be the same.

U = FunctionSpace(mesh, 'DG', 1)
U_2d = FunctionSpace(mesh2d, 'DG', 1)
func2d = Function(U_2d)
func3d = Function(U)
ex = ExpandFunctionTo3d(func2d, func3d)
ex.solve()

Parameters: input_2d (Function) – 2D source field output_3d (Function) – 3D target field elem_height – scalar Function in 3D mesh that defines the vertical element size. Needed only in the case of HDiv function spaces.
solve()[source]
class thetis.utility.FieldDict(*args, **kwargs)[source]

AttrDict that checks that all added fields have proper meta data.

Values can be either Function or Constant objects.

class thetis.utility.FrozenClass[source]

Bases: object

A class where creating a new attribute will raise an exception if _isfrozen == True

class thetis.utility.JackettEquationOfState[source]

Equation of State according of Jackett et al. (2006) for computing sea water density.

(17)$\rho = \rho'(T, S, p) + \rho_0$

$$\rho'(T, S, p)$$ is a nonlinear rational function.

Jackett et al. (2006). Algorithms for Density, Potential Temperature, Conservative Temperature, and the Freezing Temperature of Seawater. Journal of Atmospheric and Oceanic Technology, 23(12):1709-1728. http://dx.doi.org/10.1175/JTECH1946.1

a = array([ 9.99840854e+02, 7.34716259e+00, -5.32112318e-02, 3.64924391e-04, 2.58805710e+00, -6.71682828e-03, 1.92032021e-03, 1.17982637e-02, 9.89202193e-08, 4.69966428e-06, -2.58621871e-08, -3.29214140e-12])
b = array([ 1.00000000e+00, 7.28152101e-03, -4.47872655e-05, 3.38510030e-07, 1.36512024e-10, 1.76321267e-03, -8.80665833e-06, -1.88326894e-10, 5.74637767e-06, 1.47162755e-09, 6.71032463e-06, -2.44616980e-17, -9.15344176e-18])
compute_rho(s, th, p, rho0=0.0)[source]

Compute sea water density.

Parameters: s (float or numpy.array) – Salinity expressed on the Practical Salinity Scale 1978 th (float or numpy.array) – Potential temperature in Celsius, referenced to pressure p_r = 0 dbar. p (float or numpy.array) – Pressure in decibars (1 dbar = 1e4 Pa) rho0 (float) – Optional reference density. If provided computes $$\rho' = \rho(S, Th, p) - \rho_0$$ water density float or numpy.array

All pressures are gauge pressures: they are the absolute pressures minus standard atmosperic pressure 10.1325 dbar.

eval(s, th, p, rho0=0.0)[source]

Compute sea water density.

class thetis.utility.LinearEquationOfState(rho_ref, alpha, beta, th_ref, s_ref)[source]

Linear Equation of State for computing sea water density

$\rho = \rho_{ref} - \alpha (T - T_{ref}) + \beta (S - S_{ref})$
Parameters: rho_ref (float) – reference density alpha (float) – thermal expansion coefficient beta (float) – haline contraction coefficient th_ref (float) – reference temperature s_ref (float) – reference salinity
compute_rho(s, th, p, rho0=0.0)[source]

Compute sea water density.

Parameters: s (float or numpy.array) – Salinity expressed on the Practical Salinity Scale 1978 th (float or numpy.array) – Potential temperature in Celsius p (float or numpy.array) – Pressure in decibars (1 dbar = 1e4 Pa) rho0 (float) – Optional reference density. If provided computes $$\rho' = \rho(S, Th, p) - \rho_0$$ water density float or numpy.array

Pressure is ingored in this equation of state.

eval(s, th, p, rho0=0.0)[source]

Compute sea water density.

class thetis.utility.Mesh3DConsistencyCalculator(solver_obj)[source]

Bases: object

Computes a hydrostatic consistency criterion metric on the 3D mesh.

Let $$\Delta x$$ and $$\Delta z$$ denote the horizontal and vertical element sizes. The hydrostatic consistency criterion (HCC) can then be expressed as

$R = \frac{|\nabla h| \Delta x}{\Delta z} < 1$

where $$\nabla h$$ is the bathymetry gradient (or gradient of the internal horizontal facet).

Violating the hydrostatic consistency criterion leads to internal pressure gradient errors. In practice one can violate the $$R < 1$$ condition without jeopardizing numerical stability; typically $$R < 5$$. Mesh consistency can be improved by coarsening the vertical mesh, refining the horizontal mesh, or smoothing the bathymetry.

For a 3D prism, let $$\delta z_t,\delta z_b$$ denote the maximal $$z$$ coordinate difference in the surface and bottom facets, respectively, and $$\Delta z$$ the height of the prism. We can then compute $$R$$ for the two facets as

$\begin{split}R_t &= \frac{\delta z_t}{\Delta z} \\ R_b &= \frac{\delta z_b}{\Delta z}\end{split}$

For a straight prism we have $$R = 0$$, and $$R = 1$$ in the case where the highest bottom node is at the same level as the lowest surface node.

Parameters: solver_obj – FlowSolver object
solve()[source]

Compute the HCC metric

class thetis.utility.ParabolicViscosity(uv_bottom, bottom_drag, bathymetry, nu, solver_parameters={})[source]

Bases: object

Computes parabolic eddy viscosity profile assuming log layer flow

$\nu = \kappa u_{bf} \frac{(-z)(h + z_0 + z)}{h + z_0}$

with

$u_{bf} = \sqrt{C_D} |\mathbf{u}_b|$
Parameters: uv_bottom (3D Function) – bottom velocity bottom_drag (3D Function) – bottom drag field bathymetry (3D Function) – bathymetry field nu (3D Function) – eddy viscosity field solver_parameters (dict) – PETSc solver options
solve()[source]

Computes viscosity and stores it in nu field

class thetis.utility.SmagorinskyViscosity(uv, output, c_s, h_elem_size, max_val, min_val=1e-10, weak_form=True, solver_parameters={})[source]

Bases: object

Computes Smagorinsky subgrid scale horizontal viscosity

This formulation is according to Ilicak et al. (2012) and Griffies and Hallberg (2000).

$\nu = (C_s \Delta x)^2 |S|$

with the deformation rate

$\begin{split}|S| &= \sqrt{D_T^2 + D_S^2} \\ D_T &= \frac{\partial u}{\partial x} - \frac{\partial v}{\partial y} \\ D_S &= \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\end{split}$

$$\Delta x$$ is the horizontal element size and $$C_s$$ is the Smagorinsky coefficient.

To match a certain mesh Reynolds number $$Re_h$$ set $$C_s = 1/\sqrt{Re_h}$$.

Ilicak et al. (2012). Spurious dianeutral mixing and the role of momentum closure. Ocean Modelling, 45-46(0):37-58. http://dx.doi.org/10.1016/j.ocemod.2011.10.003

Griffies and Hallberg (2000). Biharmonic friction with a Smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models. Monthly Weather Review, 128(8):2935-2946. http://dx.doi.org/10.1175/1520-0493(2000)128%3C2935:BFWASL%3E2.0.CO;2

Parameters: uv_3d (3D vector Function) – horizontal velocity output (3D scalar Function) – Smagorinsky viscosity field c_s (float or Constant) – Smagorinsky coefficient h_elem_size (3D scalar Function or Constant) – field that defines the horizontal element size max_val (float) – Maximum allowed viscosity. Viscosity will be clipped at this value. min_val (float) – Minimum allowed viscosity. Viscosity will be clipped at this value. weak_form (bool) – Compute velocity shear by integrating by parts. Necessary for some function spaces (e.g. P0). solver_parameters (dict) – PETSc solver options
solve()[source]

Compute viscosity

class thetis.utility.SubFunctionExtractor(input_3d, output_2d, boundary='top', elem_facet=None, elem_height=None)[source]

Bases: object

Extract a 2D sub-function from a 3D function in an extruded mesh

Given 2D and 3D functions,

U = FunctionSpace(mesh, 'DG', 1)
U_2d = FunctionSpace(mesh2d, 'DG', 1)
func2d = Function(U_2d)
func3d = Function(U)


Get surface value:

ex = SubFunctionExtractor(func3d, func2d,
boundary='top', elem_facet='top')
ex.solve()


Get bottom value:

ex = SubFunctionExtractor(func3d, func2d,
boundary='bottom', elem_facet='bottom')
ex.solve()


Get value at the top of bottom element:

ex = SubFunctionExtractor(func3d, func2d,
boundary='bottom', elem_facet='top')
ex.solve()

Parameters: input_3d (Function) – 3D source field output_2d (Function) – 2D target field boundary (str) – ‘top’|’bottom’ Defines whether to extract from the surface or bottom 3D elements elem_facet (str) – ‘top’|’bottom’|’average’ Defines which facet of the 3D element is extracted. The ‘average’ computes mean of the top and bottom facets of the 3D element. elem_height – scalar Function in 2D mesh that defines the vertical element size. Needed only in the case of HDiv function spaces.
solve()[source]
class thetis.utility.SubdomainProjector(v, v_out, subdomain_id, solver_parameters=None, constant_jacobian=True)[source]

Bases: object

Projector that projects the restriction of an expression to the specified subdomain.

project()[source]

Apply the projection.

class thetis.utility.SumFunction[source]

Bases: object

Helper class to keep track of sum of Coefficients.

Initialize empty sum.

get operation returns Constant(0)

add(coeff)[source]

get_sum()[source]

Returns a sum of all added Coefficients

class thetis.utility.VelocityMagnitudeSolver(solution, u=None, w=None, min_val=1e-06, solver_parameters={})[source]

Bases: object

Computes magnitude of (u[0],u[1],w) and stores it in solution

Parameters: solution (Function) – scalar field for velocity magnitude scalar Function u (Function) – horizontal velocity w (Function) – vertical velocity min_val (float) – minimum value of magnitude. Minimum value of solution will be clipped to this value solver_parameters (dict) – PETSc solver options

If u is None computes magnitude of (0,0,w).

If w is None computes magnitude of (u[0],u[1],0).

solve()[source]

Compute the magnitude

class thetis.utility.VerticalIntegrator(input, output, bottom_to_top=True, bnd_value=Constant(FiniteElement('Real', None, 0), 5), average=False, bathymetry=None, elevation=None, solver_parameters={})[source]

Bases: object

Computes vertical integral (or average) of a field.

Parameters: input – 3D field to integrate output – 3D field where the integral is stored bottom_to_top – Defines the integration direction: If True integration is performed along the z axis, from bottom surface to top surface. bnd_value – Value of the integral at the bottom (top) boundary if bottom_to_top is True (False) average – If True computes the vertical average instead. Requires bathymetry and elevation fields bathymetry – 3D field defining the bathymetry elevation – 3D field defining the free surface elevation solver_parameters (dict) – PETSc solver options
solve()[source]

Computes the integral and stores it in the output field.

class thetis.utility.VerticalVelocitySolver(solution, uv, bathymetry, boundary_funcs={}, solver_parameters={})[source]

Bases: object

Computes vertical velocity diagnostically from the continuity equation

Vertical velocity is obtained from the continuity equation

(18)$\frac{\partial w}{\partial z} = -\nabla_h \cdot \textbf{u}$

and the bottom impermeability condition ($$h$$ denotes the bathymetry)

$\begin{split}\textbf{n}_h \cdot \textbf{u} + w n_z &= 0 \quad \forall \mathbf{x} \in \Gamma_{b} \\ \Leftrightarrow w &= -\nabla_h h \cdot \mathbf{u} \quad \forall \mathbf{x} \in \Gamma_{b}\end{split}$

$$w$$ can be solved with the weak form

$\begin{split}\int_{\Gamma_s} w n_z \varphi dS + \int_{\mathcal{I}_h} \text{avg}(w) \text{jump}(\varphi n_z) dS - \int_{\Omega} w \frac{\partial \varphi}{\partial z} dx = \\ \int_{\Omega} \mathbf{u} \cdot \nabla_h \varphi dx - \int_{\mathcal{I}_h \cup \mathcal{I}_v} \text{avg}(\mathbf{u}) \cdot \text{jump}(\varphi \mathbf{n}_h) dS - \int_{\Gamma_s} \mathbf{u} \cdot \varphi \mathbf{n}_h dS\end{split}$

where the $$\Gamma_b$$ terms vanish due to the bottom impermeability condition.

Parameters: solution – w Function uv – horizontal velocity Function bathymetry – bathymetry Function boundary_funcs (dict) – boundary conditions used in the 3D momentum equation. Provides external values of uv (if any). solver_parameters (dict) – PETSc solver options
solve()[source]

Compute w

thetis.utility.beta_plane_coriolis_function(latitude, out_function, y_offset=0.0)[source]

Interpolates beta plane Coriolis function to a field

Parameters: latitude (float) – latitude in degrees out_function – Function where to interpolate y_offset (float) – offset (y - y_0) used in Beta-plane approximation. A constant in mesh coordinates.
thetis.utility.beta_plane_coriolis_params(latitude)[source]

Computes beta plane parameters $$f_0,\beta$$ based on latitude

Parameters: latitude (float) – latitude in degrees f_0, beta float
thetis.utility.comp_tracer_mass_2d(eta, bath, scalar_func)[source]

Computes total tracer mass in the 2D domain :arg eta: elevation Function :arg bath: bathymetry Function :arg scalar_func: scalar Function to integrate

thetis.utility.comp_tracer_mass_3d(scalar_func)[source]

Computes total tracer mass in the 3D domain

Parameters: scalar_func – scalar Function to integrate
thetis.utility.comp_volume_2d(eta, bath)[source]

Computes volume of the 2D domain as an integral of the elevation field

thetis.utility.comp_volume_3d(mesh)[source]

Computes volume of the 3D domain as an integral

thetis.utility.compute_baroclinic_head(solver)[source]

Computes the baroclinic head $$r$$ from the density field

$r = \frac{1}{\rho_0} \int_{z}^\eta \rho' d\zeta.$
thetis.utility.compute_bottom_drag(h_b, drag)[source]

Computes bottom drag coefficient (Cd) from the law-of-the wall

$C_D = \left( \frac{\kappa}{\ln (h_b + z_0)/z_0} \right)^2$
Parameters: h_b (Function) – the height above bed where the bottom velocity is evaluated in the law-of-the-wall fit drag (Function) – field where C_D is stored
thetis.utility.compute_bottom_friction(solver, uv_3d, uv_bottom_2d, z_bottom_2d, bathymetry_2d, bottom_drag_2d)[source]

Updates bottom friction related fields for the 3D model

Parameters: solver – FlowSolver object uv_3d (3D vector Function) – horizontal velocity uv_bottom_2d (2D vector Function) – 2D bottom velocity field z_bottom_2d (2D scalar Function) – Bottom element z coordinate bathymetry_2d (2D scalar Function) – Bathymetry field bottom_drag_2d (2D scalar Function) – Bottom grad field
thetis.utility.compute_boundary_length(mesh2d)[source]

Computes the length of the boundary segments in given 2d mesh

thetis.utility.compute_elem_height(zcoord, output)[source]

Computes the element height on an extruded mesh.

Parameters: zcoord (Function) – field that contains the z coordinates of the mesh output (Function) – field where element height is stored
thetis.utility.create_directory(path, comm=<mpi4py.MPI.Intracomm object>)[source]

Create a directory on disk

Raises IOError if a file with the same name already exists.

thetis.utility.element_continuity(ufl_element)[source]

Return an ElementContinuity instance with the continuity of a given element.

Parameters: ufl_element – The UFL element to determine the continuity of. A new ElementContinuity instance.
thetis.utility.extrude_mesh_sigma(mesh2d, n_layers, bathymetry_2d, z_stretch_fact=1.0, min_depth=None)[source]

Extrudes a 2d surface mesh with bathymetry data defined in a 2d field.

Generates a uniform terrain following mesh.

Parameters: mesh2d – 2D mesh n_layers – number of vertical layers bathymetry – 2D Function of the bathymetry (the depth of the domain; positive downwards)
thetis.utility.get_facet_mask(function_space, mode='geometric', facet='bottom')[source]

Returns the top/bottom nodes of extruded 3D elements.

Parameters: function_space – Firedrake FunctionSpace object mode (str) – ‘topological’, to retrieve nodes that lie on the facet, or ‘geometric’ for nodes whose basis functions do not vanish on the facet. facet (str) – ‘top’ or ‘bottom’

Note

The definition of top/bottom depends on the direction of the extrusion. Here we assume that the mesh has been extruded upwards (along positive z axis).

thetis.utility.get_horizontal_elem_size_2d(sol2d)[source]

Computes horizontal element size from the 2D mesh

Parameters: sol2d – 2D Function where result is stored
thetis.utility.get_horizontal_elem_size_3d(sol2d, sol3d)[source]

Computes horizontal element size from the 2D mesh, then copies it on a 3D field

Parameters: sol2d – 2D Function for the element size field sol3d – 3D Function for the element size field
thetis.utility.get_zcoord_from_mesh(zcoord, solver_parameters={})[source]

Evaluates z coordinates from the 3D mesh

Parameters: zcoord – scalar Function where coordinates will be stored
thetis.utility.select_and_move_detectors(mesh, detector_locations, detector_names=None, maximum_distance=0.0)[source]

Select those detectors that are within the domain and/or move them to the nearest cell centre within the domain

Parameters: mesh – Defines the domain in which detectors are to be located detector_locations – List of x, y locations detector_names – List of detector names (optional). If provided, a list of selected locations and a list of selected detector names are returned, otherwise only a list of selected locations is returned maximum_distance – Detectors whose initial locations is outside the domain, but for which the nearest cell centre is within the specified distance, are moved to this location. By default a maximum distance of 0.0 is used, i.e no detectors are moved.
thetis.utility.tensor_jump(v, n)[source]

Jump term for vector functions based on the tensor product

$\text{jump}(\mathbf{u}, \mathbf{n}) = (\mathbf{u}^+ \mathbf{n}^+) + (\mathbf{u}^- \mathbf{n}^-)$

This is the discrete equivalent of grad(u) as opposed to the vectorial UFL jump operator ufl.jump() which represents div(u).

thetis.utility.timed_region()

Log.Event(type cls, name, klass=None)

thetis.utility.timed_stage()

Log.Stage(type cls, name)

## Module contents¶

thetis.timed_region()

Log.Event(type cls, name, klass=None)

thetis.timed_stage`()

Log.Stage(type cls, name)