2. Tidal simulation in the Orkney island¶
2.1. Introduction¶
This example demonstrates how OpenTidalFarm can be used for simulating the tides in a realistic domain.
We will be simulating the tides in the Pentland Firth, Scotland for 12.5 hours, starting at 14:40 am on the 18.9.2001. The flow result at the end of the simulation looks like:
This example requires some large data files, that must be downloaded separately by calling in the source code directory:
git submodule init
git submodule update
2.2. Implementation¶
We begin with importing the OpenTidalFarm module.
from opentidalfarm import *
We also need the datetime module for the tidal forcing.
import datetime
Next we define the UTM zone and band, as we we will need it multiple times later on.
utm_zone = 30
utm_band = 'V'
Next we create shallow water problem and attach the domain and boundary conditions
prob_params = SWProblem.default_parameters()
We load the mesh in UTM coordinates, and boundary ids from file
domain = FileDomain("../data/meshes/orkney/orkney_utm.xml")
prob_params.domain = domain
The mesh and boundary ids can be visualised with
#plot(domain.facet_ids, interactive=True)
Next we specify boundary conditions. We apply tidal boundary forcing, by using the TidalForcing class.
eta_expr = TidalForcing(grid_file_name='../data/netcdf/gridES2008.nc',
data_file_name='../data/netcdf/hf.ES2008.nc',
ranges=((-4.0,0.0), (58.0,61.0)),
utm_zone=utm_zone,
utm_band=utm_band,
initial_time=datetime.datetime(2001, 9, 18, 10, 40),
constituents=['Q1', 'O1', 'P1', 'K1', 'N2', 'M2', 'S2', 'K2'])
bcs = BoundaryConditionSet()
bcs.add_bc("eta", eta_expr, facet_id=1)
bcs.add_bc("eta", eta_expr, facet_id=2)
The free-slip boundary conditions are a special case. The boundary condition type weak_dirichlet enforces the boundary value only in the normal direction of the boundary. Hence, a zero weak Dirichlet boundary condition gives us free-slip, while a zero strong_dirichlet boundary condition would give us no-slip.
bcs.add_bc("u", Constant((0, 0)), facet_id=3, bctype="strong_dirichlet")
prob_params.bcs = bcs
Next we load the bathymetry from the NetCDF file.
bathy_expr = BathymetryDepthExpression('../data/netcdf/bathymetry.nc',
utm_zone=utm_zone, utm_band=utm_band, domain=domain.mesh)
prob_params.depth = bathy_expr
The bathymetry can be visualised with
#plot(bathy_expr, mesh=domain.mesh, title="Bathymetry", interactive=True)
Equation settings
For stability reasons, we want to increase the viscosity at the inflow and outflow boundary conditions. For that, we read in a precomputed function (generated by `compute_distance`) and
V = FunctionSpace(domain.mesh, "CG", 1)
dist = Function(V)
File("dist.xml") >> dist
With that we can define an expression that evaluates to a nu_inside value inside the domain and a nu_outside value near the in/outflow boundary.
class ViscosityExpression(Expression):
def __init__(self, dist_function, dist_threshold, nu_inside, nu_boundary):
self.dist_function = dist_function
self.nu_inside = nu_inside
self.nu_boundary = nu_boundary
self.dist_threshold = dist_threshold
def eval(self, value, x):
if self.dist_function(x) > self.dist_threshold:
value[0] = self.nu_inside
else:
value[0] = self.nu_boundary
Finally, we interpolate this expression to a piecewise discontinuous, constant function and attach it as the viscosity value to the shallow water problem.
W = FunctionSpace(domain.mesh, "DG", 0)
nu = ViscosityExpression(dist, dist_threshold=1000, nu_inside=10., nu_boundary=1e3)
nu_func = interpolate(nu, W)
prob_params.viscosity = nu_func
The other parameters are set as usual.
prob_params.friction = Constant(0.0025)
# Temporal settings
prob_params.start_time = Constant(0)
prob_params.finish_time = Constant(12.5*60*60)
prob_params.dt = Constant(5*60)
# The initial condition consists of three components: u_x, u_y and eta.
# Note that we set the velocity components to a small positive number, as some
# components of the Jacobian of the quadratic friction term is
# non-differentiable.
prob_params.initial_condition = Constant((DOLFIN_EPS, DOLFIN_EPS, 1))
# Now we can create the shallow water problem
problem = SWProblem(prob_params)
# Next we create a shallow water solver. Here we choose to solve the shallow
# water equations in its fully coupled form:
sol_params = CoupledSWSolver.default_parameters()
solver = CoupledSWSolver(problem, sol_params)
Now we are ready to solve
f_state = File("results/state.pvd")
timer = Timer('')
for s in solver.solve():
t = float(s["time"])
log(INFO, "Computed solution at time %f in %f s." % (t, timer.stop()))
f_state << (s["state"], t)
timer.start()
The code for this example can be found in examples/tidal-simulation/ in the OpenTidalFarm source tree, and executed as follows:
$ mpirun -n 4 python orkney-coupled.py
where 4 should be replaced by the number of CPU cores available.