2.3. Solvers

A set of solvers for OpenTidalFarm

2.3.1. Coupled shallow water solver

class opentidalfarm.solvers.coupled_sw_solver.CoupledSWSolver(problem, solver_params)

The coupled solver solves the shallow water equations as a fully coupled system with a \(\theta\)-timestepping discretization.

For a opentidalfarm.problems.steady_sw.SteadySWProblem, it solves \(u\) and \(\eta\) such that:

\[\begin{split}u \cdot \nabla u - \nu \Delta u + g \nabla \eta + \frac{c_b + c_t}{H} \| u\| u & = f_u, \\ \nabla \cdot \left(H u\right) & = 0.\end{split}\]

For a opentidalfarm.problems.multi_steady_sw.MultiSteadySWProblem, it solves \(u^{n+1}\) and \(\eta^{n+1}\) for each timelevel (including the initial time) such that:

\[\begin{split}u^{n+1} \cdot \nabla u^{n+1} - \nu \Delta u^{n+1} + g \nabla \eta^{n+1} + \frac{c_b + c_t}{H^{n+1}} \| u^{n+1}\| u^{n+1} & = f_u^{n+1}, \\ \nabla \cdot \left(H^{n+1} u^{n+1}\right) & = 0.\end{split}\]

For a opentidalfarm.problems.sw.SWProblem, and given an initial condition \(u^{0}\) and \(\eta^{0}\) it solves \(u^{n+1}\) and \(\eta^{n+1}\) for each timelevel such that:

\[\begin{split}\frac{1}{\Delta t} \left(u^{n+1} - u^{n}\right) + u^{n+\theta} \cdot \nabla u^{n+\theta} - \nu \Delta u^{n+\theta} + g \nabla \eta^{n+\theta} + \frac{c_b + c_t}{H^{n+\theta}} \| u^{n+\theta}\| u^{n+\theta} & = f_u^{n+\theta}, \\ \frac{1}{\Delta t} \left(\eta^{n+1} - \eta^{n}\right) + \nabla \cdot \left(H^{n+\theta} u^{n+\theta}\right) & = 0.\end{split}\]

with \(\theta \in [0, 1]\) and \(u^{n+\theta} := \theta u^{n+1} + (1-\theta) u^n\) and \(\eta^{n+\theta} := \theta \eta^{n+1} + (1-\theta) \eta^n\).

Furthermore:

• \(u\) is the velocity,

• \(\eta\) is the free-surface displacement,

• \(H=\eta + h\) is the total water depth where \(h\) is the water depth at rest,

• \(f_u\) is the velocity forcing term,

• \(c_b\) is the (quadratic) natural bottom friction coefficient,

• \(c_t\) is the (quadratic) friction coefficient due to the turbine farm,

• \(\nu\) is the viscosity coefficient,

• \(g\) is the gravitational constant,

• \(\Delta t\) is the time step.

Some terms, such as the advection or the diffusion term, may be deactivated in the problem specification.

The resulting equations are solved with Newton-iterations and the linear problems solved as specified in the dolfin_solver setting in CoupledSWSolverParameters.

static default_parameters()

Return the default parameters for the CoupledSWSolver.

solve(annotate=True)

Returns an iterator for solving the shallow water equations.

class opentidalfarm.solvers.coupled_sw_solver.CoupledSWSolverParameters

A set of parameters for a CoupledSWSolver.

Following parameters are available:

Variables:

• dolfin_solver –

  The dictionary with parameters for the dolfin Newton solver. A list of valid entries can be printed with:

  info(NonlinearVariationalSolver.default_parameters(), True)
  

  By default, the MUMPS direct solver is used for the linear system. If not available, the default solver and preconditioner of FEniCS is used.

• dump_period – Specifies how often the solution should be dumped to disk. Use a negative value to disable it. Default 1.

• cache_forward_state – If True, the shallow water solutions are stored for every timestep and are used as initial guesses for the next solve. If False, the solution of the previous timestep is used as an initial guess. Default: True

• print_individual_turbine_power – Print out the turbine power for each turbine. Default: False

• quadrature_degree – The quadrature degree for the matrix assembly. Default: -1 (automatic)

• cpp_flags – A list of cpp compiler options for the code generation. Default: [“-O3”, “-ffast-math”, “-march=native”]

• revolve_parameters – The adjoint checkpointing settings as a set of the form (strategy, snaps_on_disk, snaps_in_ram, verbose). Default: None

• output_dir – The base directory in which to store the file ouputs. Default: os.curdir

• output_turbine_power – Output the power generation of the individual turbines. Default: False

• output_j – Output the evaluation of the choosen functional (e.g power) for each iteration and store it in a numpy textfile. (This is the same j that is printed each iteration, when the FEniCS log level is INFO or above.) Default: False

• output_temporal_breakdown_of_j – Output the j at each timestep and at iteration. The output is stored in a numpy textfile. (This is the same temporal breakdown which is shown when the FEniCS log level is INFO or above.) Default: False

• output_abs_u_at_turbine_positions – Output the absolute value of the velocity at each turbine position. Default: False

• output_control_array – Output a numpy textfile containing the control array from each optimisation and search iteration. Default: False

• callback – A callback function that is executed for every time-level. The callback function must take a single parameter which contains the dictionary with the solution variables.

2.3.2. Pressure-correction shallow water solver

class opentidalfarm.solvers.ipcs_sw_solver.IPCSSWSolver(problem, parameters)

This incremental pressure correction scheme (IPCS) is an operator splitting scheme that follows the idea of Goda [1] and Simo [2]. This scheme preserves the exact same stability properties as Navier-Stokes and hence does not introduce additional dissipation in the flow (FIXME: needs verification).

The idea is to replace the unknown free-surface with an approximation. This is chosen as the free-surface solution from the previous solution.

The time discretization is done using a \(\theta\)-scheme, the convection, friction and divergence are handled semi-implicitly. Thus, we have a discretized version of the shallow water equations as

\[\begin{split}\frac{1}{\Delta t}\left( u^{n+1}-u^{n} \right)-\nabla\cdot\nu\nabla u^{n+\theta}+u^*\cdot\nabla u^{n+\theta}+g\nabla \eta^{n+\theta} + \frac{c_b + c_t}{H^n} \| u^n\| u^{n+\theta} &= f_u^{n+\theta}, \\ \frac{1}{\Delta t}\left( \eta^{n+1}-\eta^{n} \right) + \nabla \cdot \left( H^{n} u^{n+\theta} \right) &= 0,\end{split}\]

where \(\square^{n+\theta} = \theta{\square}^{n+1}+(1-\theta)\square^n, \theta \in [0, 1]\) and \(u^* = \frac{3}{2}u^n - \frac{1}{2}u^{n-1}\).

This convection term is unconditionally stable, and with \(\theta=0.5\), this equation is second order in time and space [2] (FIXME: Needs verification).

For the operator splitting, we use the free-surface solution from the previous timestep as an estimation, giving an equation for a tentative velocity, \(\tilde{u}^{n+1}\):

\[\frac{1}{\Delta t}\left( \tilde{u}^{n+1}-u^{n} \right) - \nabla\cdot\nu\nabla \tilde{u}^{n+\theta} + u^*\cdot\nabla \tilde u^{n+\theta}+g\nabla \eta^{n} + \frac{c_b + c_t}{H^n} \| u^n\| \tilde u^{n+\theta} = f_u^{n+\theta}.\]

This tenative velocity does not satisfy the divergence equation, and thus we define a velocity correction \(u^c=u^{n+1}-\tilde{u}^{n+1}\). Substracting the second equation from the first, we see that

\[\begin{split}\frac{1}{\Delta t}u^c - \theta \nabla\cdot\nu\nabla u^c + \theta u^*\cdot\nabla u^{c} + g\theta \nabla\left( \eta^{n+1} - \eta^n\right) + \theta \frac{c_b + c_t}{H^n} \| u^n\| u^{c} &= 0, \\ \frac{1}{\Delta t}\left( \eta^{n+1}-\eta^{n} \right) + \theta \nabla \cdot \left( H^{n} u^c \right) &= -\nabla \cdot \left( H^{n} \tilde{u}^{n+\theta} \right).\end{split}\]

The operator splitting is a first order approximation, \(O(\Delta t)\), so we can, without reducing the order of the approximation simplify the above to

\[\begin{split}\frac{1}{\Delta t}u^c + g\theta \nabla\left( \eta^{n+1} - \eta^n\right) &= 0, \\ \frac{1}{\Delta t}\left( \eta^{n+1}-\eta^{n} \right) + \theta \nabla \cdot \left( H^{n} u^c \right) &= -\nabla \cdot \left( H^{n} \tilde{u}^{n+\theta} \right),\end{split}\]

which is reducible to the problem:

\[\eta^{n+1}-\eta^{n} - g \Delta t^2 \theta^2 \nabla \cdot \left( H^{n+1} \nabla\left( \eta^{n+1} - \eta^n\right) \right) = -\Delta t \nabla \cdot \left( H^{n} \tilde{u}^{n+\theta} \right).\]

The corrected velocity is then easily calculated from

\[u^{n+1} = \tilde{u}^{n+1} - \Delta tg\theta\nabla\left(\eta^{n+1}-\eta^n\right)\]

The scheme can be summarized in the following steps:

1. Replace the pressure with a known approximation and solve for a tenative velocity \(\tilde u^{n+1}\).

2. Solve a free-surface correction equation for the free-surface, \(\eta^{n+1}\)

3. Use the corrected pressure to find the velocity correction and calculate \(u^{n+1}\)

4. Update t, and repeat.

Remarks:

• This solver only works with transient problems, that is with opentidalfarm.problems.sw.SWProblem.

• This solver supports large eddy simulation (LES). The LES model is implemented via the opentidalfarm.solvers.les.LES class.

[1]

Goda, Katuhiko. A multistep technique with implicit difference schemes for calculating two-or three-dimensional cavity flows. Journal of Computational Physics 30.1 (1979): 76-95.

[2] (1,2)

Simo, J. C., and F. Armero. Unconditional stability and long-term behavior of transient algorithms for the incompressible Navier-Stokes and Euler equations. Computer Methods in Applied Mechanics and Engineering 111.1 (1994): 111-154.

static default_parameters()

Returns a dictionary with the default parameters.

solve(annotate=True)

Solve the shallow water equations

class opentidalfarm.solvers.ipcs_sw_solver.IPCSSWSolverParameters

A set of parameters for a IPCSSWSolver.

Performance parameters:

Variables:

• dolfin_solver –

  The dictionary with parameters for the dolfin Newton solver. A list of valid entries can be printed with:

  info(NonlinearVariationalSolver.default_parameters(), True)
  

  By default, the MUMPS direct solver is used for the linear system. If not availabe, the default solver and preconditioner of FEniCS is used.

• quadrature_degree – The quadrature degree for the matrix assembly. Default: -1 (auto)

• cpp_flags – A list of cpp compiler options for the code generation. Default: [“-O3”, “-ffast-math”, “-march=native”]

Large eddy simulation parameters:

Variables:

• les_model – De-/Activates the LES model. Default: True

• les_model_parameters – A dictionary with parameters for the LES model. Default: {‘smagorinsky_coefficient’: 1e-2}