Allow multiple tidal energy farms per subdomain (#358)

- Adds list of farm options for each subdomain - Updates corresponding examples accordingly - Farm options now need to be in a list, even if there is only one entry (warning message included for this) - Added instantaneous_power to the farm power callback for convenience --------- Co-authored-by: Connor Jordan <[email protected]> Co-authored-by: Stephan Kramer <[email protected]>
thetisproject · Feb 28, 2024 · 9572dc7 · 9572dc7
1 parent 70aa071
commit 9572dc7
Show file tree

Hide file tree

Showing 6 changed files with 32 additions and 17 deletions.
diff --git a/examples/discrete_turbines/channel-optimisation.py b/examples/discrete_turbines/channel-optimisation.py
@@ -106,7 +106,7 @@
 # everywhere, you can use the key "everywhere" instead of 2 below.
 # Regardless of this, the area where the turbines are places need to have sufficiently high
 # resolution to resolve the Gaussian bump functions that represent the turbines
-options.discrete_tidal_turbine_farms[2] = farm_options
+options.discrete_tidal_turbine_farms[2] = [farm_options]
 
 # set initial condition (initial velocity should not be exactly 0 to avoid failures in the Newton solve)
 solver_obj.assign_initial_conditions(uv=(as_vector((1e-3, 0.0))))

diff --git a/examples/discrete_turbines/tidal_array.py b/examples/discrete_turbines/tidal_array.py
@@ -83,7 +83,7 @@
 options.wetting_and_drying_alpha = Constant(0.5)
 if not hasattr(options.swe_timestepper_options, 'use_automatic_timestep'):
     options.timestep = 50.0
-options.discrete_tidal_turbine_farms[site_ID] = farm_options
+options.discrete_tidal_turbine_farms[site_ID] = [farm_options]
 
 # Use direct solver instead of default iterative settings
 # (see SemiImplicitSWETimeStepperOptions2d in thetis/options.py)
@@ -111,14 +111,15 @@
 # Operation of tidal turbine farm through a callback
 cb_turbines = turbines.TurbineFunctionalCallback(solver_obj)
 solver_obj.add_callback(cb_turbines, 'timestep')
-powers = []  # create empty list to append to if we would like to plot power over time
+powers = []  # create empty list to append instantaneous powers to
 
 
 def update_forcings(t_new):
     ramp = tanh(t_new / 2000.)
     tidal_vel.project(Constant(ramp * 3.))
-    powers.append(cb_turbines.integrated_power[0] / 100 - sum(powers))
+    powers.append(cb_turbines.instantaneous_power[0])
 
 
 # See channel-optimisation example for a completely steady state simulation (no ramp)
 solver_obj.iterate(update_forcings=update_forcings)
+powers.append(cb_turbines.instantaneous_power[0])  # add final power, should be the same as callback hdf5 file!
diff --git a/examples/tidalfarm/tidalfarm.py b/examples/tidalfarm/tidalfarm.py
@@ -107,7 +107,7 @@ def update_forcings(t):
 # in such a way that the cost is expressed in kW which can be subtracted from the profit
 # which is calculated as the power extracted by the turbines
 farm_options.break_even_wattage = 200
-options.tidal_turbine_farms[2] = farm_options
+options.tidal_turbine_farms[2] = [farm_options]
 
 # we first run the "forward" model with no turbines
 turbine_density.assign(0.0)

diff --git a/thetis/options.py b/thetis/options.py
@@ -881,12 +881,13 @@ class ModelOptions2d(CommonModelOptions):
 
         Note this is only relevant if `use_automatic_wetting_and_drying_alpha` is set to ``True``.
         """).tag(config=True)
-    tidal_turbine_farms = Dict(trait=TidalTurbineFarmOptions(),
-                               default_value={}, help='Dictionary mapping subdomain ids to the options of the corresponding farm')
+    tidal_turbine_farms = Dict(trait=List(TidalTurbineFarmOptions()), default_value={},
+                               help='Dictionary mapping subdomain ids to a list of TidalTurbineFarmOptions instances '
+                                    'corresponding to one or more farms.')
 
-    discrete_tidal_turbine_farms = Dict(trait=DiscreteTidalTurbineFarmOptions(),
-                                        default_value={},
-                                        help='Dictionary mapping subdomain ids to the options of the corresponding farm')
+    discrete_tidal_turbine_farms = Dict(trait=List(DiscreteTidalTurbineFarmOptions()), default_value={},
+                                        help='Dictionary mapping subdomain ids to a list of DiscreteTidalTurbineFarmOptions '
+                                             'instances corresponding to one or more farms.')
 
     check_tracer_conservation = Bool(
         False, help="""

diff --git a/thetis/solver2d.py b/thetis/solver2d.py
@@ -463,13 +463,24 @@ def create_equations(self):
             self.tidal_farms = []
             p = self.function_spaces.U_2d.ufl_element().degree()
             quad_degree = 2*p + 1
-            for subdomain, farm_options in self.options.tidal_turbine_farms.items():
-                fdx = dx(subdomain, degree=quad_degree)
-                self.tidal_farms.append(TidalTurbineFarm(farm_options.turbine_density,
-                                                         fdx, farm_options))
-            for subdomain, farm_options in self.options.discrete_tidal_turbine_farms.items():
-                fdx = dx(subdomain, degree=farm_options.quadrature_degree)
-                self.tidal_farms.append(DiscreteTidalTurbineFarm(self.mesh2d, fdx, farm_options))
+            for subdomain, farm_options_list in self.options.tidal_turbine_farms.items():
+                if not isinstance(farm_options_list, list):
+                    error_msg = "Farm options must be entered as a list e.g. " \
+                                "solver2d.FlowSolver2d(mesh2d, bathymetry_2d).options.discrete_tidal_turbine_farms[site_ID] = " \
+                                "[farm_options]"
+                    raise TypeError(error_msg)
+                for farm_options in farm_options_list:
+                    fdx = dx(subdomain, degree=quad_degree)
+                    self.tidal_farms.append(TidalTurbineFarm(farm_options.turbine_density, fdx, farm_options))
+            for subdomain, farm_options_list in self.options.discrete_tidal_turbine_farms.items():
+                if not isinstance(farm_options_list, list):
+                    error_msg = "Farm options must be entered as a list e.g. " \
+                                "solver2d.FlowSolver2d(mesh2d, bathymetry_2d).options.discrete_tidal_turbine_farms[site_ID] = " \
+                                "[farm_options]"
+                    raise TypeError(error_msg)
+                for farm_options in farm_options_list:
+                    fdx = dx(subdomain, degree=farm_options.quadrature_degree)
+                    self.tidal_farms.append(DiscreteTidalTurbineFarm(self.mesh2d, fdx, farm_options))
         else:
             self.tidal_farms = None
         # Shallow water equations for hydrodynamic modelling

diff --git a/thetis/turbines.py b/thetis/turbines.py
@@ -183,6 +183,7 @@ def __init__(self, solver_obj, **kwargs):
             print_output('Number of turbines = {}'.format(sum(self.cost)))
         self.break_even_wattage = [farm.break_even_wattage for farm in self.farms]
 
+        self.instantaneous_power = [0] * nfarms
         # time-integrated quantities:
         self.integrated_power = [0] * nfarms
         self.average_power = [0] * nfarms
@@ -196,6 +197,7 @@ def _evaluate_timestep(self):
         for i, farm in enumerate(self.farms):
             power = farm.power_output(self.uv, self.depth)
             current_power.append(power)
+            self.instantaneous_power[i] = power
             self.integrated_power[i] += power * self.dt
             self.average_power[i] = self.integrated_power[i] / self.time_period
             self.average_profit[i] = self.average_power[i] - self.break_even_wattage[i] * self.cost[i]