7. Sensitivity analysis

7.1. Introduction

Gradient information may also be used to analyse the sensitivity of the chosen functional to various model parameters - for example the bottom friction:

../../_images/sens_friction.png

This enables the designer to identify which parameters may have a high impact upon the quality of the chosen design (as judged by the choice of functional).

This example shows how to:

  • set up shallow water solver with a turbine farm;

  • define an objective, here the farms power production;

  • compute the sensitivity of this objective with respect to:

    • viscosity;

    • depth;

    • bottom friction.

7.2. Implementation

As with other examples, we begin by defining a steady state shallow water problem, once more this is similar to the Sinusoidal wave in a channel example except that we define steady flow driven by a 0.1 m head difference across the domain. Note that this is facilitated by defining a strong dirichlet boundary condiditon on the walls:

from opentidalfarm import *

# Create a rectangular domain.
domain = FileDomain("mesh/mesh.xml")

# Specify boundary conditions.
bcs = BoundaryConditionSet()
bcs.add_bc("eta", Constant(0.1), facet_id=1)
bcs.add_bc("eta", Constant(0), facet_id=2)
# The no-slip boundary conditions.
bcs.add_bc("u", Constant((0, 0)), facet_id=3, bctype="strong_dirichlet")

Set the shallow water parameters, since we want to extract the spatial variation of the sensitivity of our model parameters, rather than simply defining constant values, we define fields over the domain (in this case we’re simply defining a field of a constant value, using dolfin’s interpolate.

prob_params = SteadySWProblem.default_parameters()
prob_params.domain = domain
prob_params.bcs = bcs
V = FunctionSpace(domain.mesh, "CG", 1)
prob_params.viscosity = interpolate(Constant(3), V)
prob_params.depth = interpolate(Constant(50), V)
prob_params.friction = interpolate(Constant(0.0025), V)

The next step is to specify the array design for which we wish to analyse the sensitivity. For simplicity we will use the starting guess from the Farm layout optimization example; 32 turbines in a regular grid layout. As before we’ll use the default turbine type and define the diameter and friction. In practice, one is likely to want to analyse the sensitivity of the optimised array layout - so one would substitute this grid layout with the optimised one.

turbine = BumpTurbine(diameter=20.0, friction=12.0)
farm = RectangularFarm(domain, site_x_start=160, site_x_end=480,
                       site_y_start=80, site_y_end=240, turbine=turbine)
farm.add_regular_turbine_layout(num_x=8, num_y=4)
prob_params.tidal_farm = farm

Now we can create the shallow water problem

problem = SteadySWProblem(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()
sol_params.dump_period = -1
solver = CoupledSWSolver(problem, sol_params)

We wish to study the effect that various model parameters have on the power. Thus, we select the PowerFunctional

functional = PowerFunctional(problem)

First let’s compute the sensitivity of the power with respect to the turbine positions. So we set the “control” variable to the turbine positions by using TurbineFarmControl. We then intialise the ReducedFunctional

control = TurbineFarmControl(farm)
rf_params = ReducedFunctional.default_parameters()
rf = ReducedFunctional(functional, control, solver, rf_params)
m0 = rf.solver.problem.parameters.tidal_farm.control_array
j = rf.evaluate(m0)
turbine_location_sensitivity = rf.derivative(m0)

print "j for turbine positions: ", j
print "dj w.r.t. turbine positions: ", turbine_location_sensitivity

Next we compute the sensitivity of the power with respect to bottom friction. We redefine the control variable using the class Control into which we pass the parameter of interest

control = Control(prob_params.friction)

Turbine positions are stored in different data structures (numpy arrays) than functions such as bottom friction (dolfin functions), so we need to use a different reduced functional; the FenicsReducedFunctional

rf = FenicsReducedFunctional(functional, control, solver)
j = rf.evaluate()
dj = rf.derivative(project=True)[0]
plot(dj, interactive=True, title="Sensitivity with respect to friction")
file = File("FrictionSensitivity.pvd")
file << dj

Now compute the sensitivity with respect to depth

control = Control(prob_params.depth)
rf = FenicsReducedFunctional(functional, control, solver)
j = rf.evaluate()
dj = rf.derivative(project=True)[0]
print "j with depth = 50 m: ", j
plot(dj, interactive=True, title="Sensitivity with respect to depth at 50m")
file = File("50mDepthSensitivity.pvd")
file << dj
../../_images/sens_depth50.png

Let’s reduce the depth and reevaluate derivative

prob_params.depth.assign(Constant(10))
j = rf.evaluate()
dj = rf.derivative(project=True)[0]
print "j with depth = 10 m: ", j
plot(dj, interactive=True, title="Sensitivity with respect to depth at 10m")
file = File("10mDepthSensitivity.pvd")
file << dj
../../_images/sens_depth10.png

Finally let’s compute the sensitivity with respect to viscosity

control = Control(prob_params.viscosity)
rf = FenicsReducedFunctional(functional, control, solver)
j = rf.evaluate()
dj = rf.derivative(project=True)[0]
plot(dj, title="Sensitivity with respect to viscosity")
file = File("ViscositySensitivity.pvd")
file << dj

interactive() # (For docker users the plots will not be shown, but they are
              #  stored in pvd-files. These files can view with ParaView.)
../../_images/sens_visc.png

7.3. How to run the example

The example code can be found in examples/channel-sensitivity/ in the OpenTidalFarm source tree, and executed as follows:

$ python channel-sensitivity.py