Discrete turbines callback example

examples/discrete_turbines/tidal_array.py

@@ -3,10 +3,13 @@
 of flow past a headland. Turbines are based on the SIMEC Atlantis AR2000
 with cut-in, rated and cut-out speeds of 1m/s, 3.05m/s and 5m/s respectively.
 Flow becomes steady after an initial ramp up.
+Two callbacks are used - one for the farm and one for discrete turbines.
 """
 
 from thetis import *
 from firedrake.output.vtk_output import VTKFile
+from turbine_callback import TidalPowerCallback
+
 
 # Set output directory, load mesh, set simulation export and end times
 outputdir = 'outputs'
@@ -17,7 +20,7 @@
 print_output('Loaded mesh ' + mesh2d.name)
 print_output('Exporting to ' + outputdir)
 
-t_end = 2 * 3600
+t_end = 1.5 * 3600
 t_export = 200.0
 
 if os.getenv('THETIS_REGRESSION_TEST') is not None:
@@ -34,15 +37,15 @@
 
 # Turbine options
 turbine_thrust_def = 'table'  # 'table' or 'constant'
-include_support_structure = True  # choose whether we want to add additional drag due to the support structure
+include_support_structure = True  # add additional drag due to the support structure?
 
 # Define the thrust curve of the turbine using a tabulated approach:
-# thrusts_AR2000 contains the values for the thrust coefficient of an AR2000 tidal turbine at corresponding speeds in
-# speeds_AR2000 which have been determined using a curve fitting technique based on:
-# cut-in speed = 1m/s
-# rated speed = 3.05m/s
-# cut-out speed = 5m/s
-# There is a ramp up and down to cut-in and at cut-out speeds for model stability.
+# speeds_AR2000: speeds for corresponding thrust coefficients - thrusts_AR2000
+# thrusts_AR2000: list of idealised thrust coefficients of an AR2000 tidal turbine using a curve fitting technique with:
+#   * cut-in speed = 1 m/s
+#   * rated speed = 3.05 m/s
+#   * cut-out speed = 5 m/s
+# (ramp up and down to cut-in and at cut-out speeds for model stability)
 speeds_AR2000 = [0., 0.75, 0.85, 0.95, 1., 3.05, 3.3, 3.55, 3.8, 4.05, 4.3, 4.55, 4.8, 5., 5.001, 5.05, 5.25, 5.5, 5.75,
                  6.0, 6.25, 6.5, 6.75, 7.0]
 thrusts_AR2000 = [0.010531, 0.032281, 0.038951, 0.119951, 0.516484, 0.516484, 0.387856, 0.302601, 0.242037, 0.197252,
@@ -59,7 +62,7 @@
     farm_options.turbine_options.thrust_coefficient = 0.6
 if include_support_structure:
     farm_options.turbine_options.C_support = 0.7  # support structure thrust coefficient
-    farm_options.turbine_options.A_support = 2.6*14.0  # cross-sectional area of support structure
+    farm_options.turbine_options.A_support = 2.6 * 14.0  # cross-sectional area of support structure
 farm_options.turbine_options.diameter = 20
 farm_options.upwind_correction = True  # See https://arxiv.org/abs/1506.03611 for more details
 turbine_density = Function(FunctionSpace(mesh2d, "CG", 1), name='turbine_density').assign(0.0)
@@ -101,26 +104,57 @@
 solver_obj.bnd_functions['shallow_water'] = {right_tag: {'un': tidal_vel},
                                              left_tag: {'elev': tidal_elev}}
 
-# initial conditions, piecewise linear function
+# Initial conditions, piecewise linear function
 elev_init = Function(P1_2d)
 elev_init.assign(0.0)
 solver_obj.assign_initial_conditions(elev=elev_init, uv=(as_vector((1e-3, 0.0))))
 
 print_output(str(options.swe_timestepper_type) + ' solver options:')
 print_output(options.swe_timestepper_options.solver_parameters)
 
-# Operation of tidal turbine farm through a callback
-cb_turbines = turbines.TurbineFunctionalCallback(solver_obj)
+# Operation of tidal turbine farm through a callback (density assumed = 1000kg/m3)
+# 1. In-built farm callback
+cb_farm = turbines.TurbineFunctionalCallback(solver_obj)
+solver_obj.add_callback(cb_farm, 'timestep')
+power_farm = []  # create empty list to append instantaneous powers to
+
+# 2. Tracking of individual turbines through the callback defined in turbine_callback.py
+#    This can be used even if the turbines are not included in the simulation.
+#    This method becomes slow when the domain becomes very large (e.g. 100k+ elements).
+turbine_names = ['TTG' + str(i+1) for i in range(len(farm_options.turbine_coordinates))]
+cb_turbines = TidalPowerCallback(solver_obj, site_ID, farm_options, ['pow'], name='array', turbine_names=turbine_names)
 solver_obj.add_callback(cb_turbines, 'timestep')
-powers = []  # create empty list to append instantaneous powers to
+powers_turbines = []
+
+print_output("\nMonitoring the following turbines:")
+for i in range(len(cb_turbines.variable_names)):
+    print_output(cb_turbines.variable_names[i] + ': (' + str(cb_turbines.get_turbine_coordinates[i][0]) + ', '
+                 + str(cb_turbines.get_turbine_coordinates[i][1]) + ')')
+print_output("")
 
 
 def update_forcings(t_new):
     ramp = tanh(t_new / 2000.)
     tidal_vel.project(Constant(ramp * 3.))
-    powers.append(cb_turbines.instantaneous_power[0])
+    power_farm.append(cb_farm.instantaneous_power[0])
+    powers_turbines.append(cb_turbines())
 
 
 # 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!
+
+power_farm.append(cb_farm.instantaneous_power[0])  # add final power, should be the same as callback hdf5 file!
+farm_energy = sum(power_farm) * options.timestep / 3600
+powers_turbines.append(cb_turbines())  # add final power, should be the same as callback hdf5 file!
+turbine_energies = np.sum(np.array(powers_turbines), axis=0) * options.timestep / 3600
+callback_diff = 100 * (farm_energy - np.sum(turbine_energies, axis=0)) / farm_energy
+
+print_output("\nTurbine energies:")
+for i in range(len(cb_turbines.variable_names)):
+    print_output(cb_turbines.variable_names[i] + ': ' + "{:.2f}".format(turbine_energies[i][0]) + 'kWh')
+
+print_output("Check callbacks sum to the same energy")
+print_output("Farm callback total energy recorded: " + "{:.2f}".format(farm_energy) + 'kWh')
+print_output("Turbines callback total energy recorded: " + "{:.2f}".format(np.sum(turbine_energies, axis=0)[0]) + 'kWh')
+print_output("Difference: " + "{:.5f}".format(callback_diff[0]) + "%")
+print_output("")

examples/discrete_turbines/turbine_callback.py

@@ -0,0 +1,138 @@
+from thetis import *
+from thetis.turbines import DiscreteTidalTurbineFarm
+
+
+class TidalPowerCallback(DiagnosticCallback):
+    """
+    DERIVED Callback to evaluate tidal stream power at the specified locations
+    Based on Thetis' standard `DetectorsCallback`
+    """
+    def __init__(self, solver_obj,
+                 idx,
+                 farm_options,
+                 field_names,
+                 name,
+                 turbine_names=None,
+                 **kwargs):
+        """
+        :arg solver_obj: Thetis solver object
+        :arg idx: Index (int), physical ID of farm
+        :arg farm_options: Farm configuration (farm_option object)
+        :arg field_names: List of fields to be interpolated, e.g. `['pow']`
+        :arg name: Unique name for this callback and its associated set of
+            locations. This determines the name of the output HDF5 file
+            (prefixed with `diagnostic_`).
+        :arg turbine_names (optional): List of turbine names (otherwise
+            named loca0, loca1, ..., locaN)
+        :arg kwargs: any additional keyword arguments, see
+            :class:`.DiagnosticCallback`.
+        """
+        # printing all location output to log is probably not a useful default:
+        kwargs.setdefault('append_to_log', False)
+        self.field_dims = []
+        for field_name in field_names:
+            if field_name != 'pow':
+                self.field_dims.append(solver_obj.fields[field_name].function_space().value_size)
+        attrs = {
+            # use null-padded ascii strings, dtype='U' not supported in hdf5,
+            # see http://docs.h5py.org/en/latest/strings.html
+            'field_names': numpy.array(field_names, dtype='S'),
+            'field_dims': self.field_dims,
+        }
+        super().__init__(solver_obj, array_dim=sum(self.field_dims), attrs=attrs, **kwargs)
+
+        locations = farm_options.turbine_coordinates
+        if turbine_names is None:
+            turbine_names = ['loca{:0{fill}d}'.format(i, fill=len(str(len(locations))))
+                             for i in range(len(locations))]
+        self.loc_names = turbine_names
+        self._variable_names = self.loc_names
+        self.locations = locations
+
+        # similar to solver2d.py
+        p = solver_obj.function_spaces.U_2d.ufl_element().degree()
+        quad_degree = 2*p + 1
+        fdx = dx(idx, degree=quad_degree)
+
+        self.farm = DiscreteTidalTurbineFarm(solver_obj.mesh2d, fdx, farm_options)
+        self.field_names = field_names
+        self._name = name
+
+        # Disassemble density field
+        xyvec = SpatialCoordinate(self.solver_obj.mesh2d)
+        loc_dens = []
+        radius = farm_options.turbine_options.diameter / 2
+        for (xo, yo) in locations:
+            dens = self.farm.turbine_density
+            dens = conditional(And(lt(abs(xyvec[0]-xo), radius), lt(abs(xyvec[1]-yo), radius)), dens, 0)
+            loc_dens.append(dens)
+
+        self.loc_dens = loc_dens
+
+    @property
+    def name(self):
+        return self._name
+
+    @property
+    def variable_names(self):
+        return self.loc_names
+
+    @property
+    def get_turbine_coordinates(self):
+        """
+        Returns a list of turbine locations as x, y coordinates instead
+        of Firedrake constants.
+        """
+        turbine_coordinates = []
+
+        for loc in self.locations:
+            x_val, y_val = loc
+            turbine_coordinates.append([x_val.values()[0], y_val.values()[0]])
+
+        return numpy.array(turbine_coordinates)
+
+    def _values_per_field(self, values):
+        """
+        Given all values evaluated in a detector location, return the values per field
+        """
+        i = 0
+        result = []
+        for dim in self.field_dims:
+            result.append(values[i:i+dim])
+            i += dim
+        return result
+
+    def message_str(self, *args):
+        return '\n'.join(
+            'In {}: '.format(name) + ', '.join(
+                '{}={}'.format(field_name, field_val)
+                for field_name, field_val in zip(self.field_names, self._values_per_field(values)))
+            for name, values in zip(self.loc_names, args))
+
+    def _evaluate_field(self, field_name):
+        if field_name == 'pow':  # intended use
+            _uv, _eta = split(self.solver_obj.fields.solution_2d)
+            _depth = self.solver_obj.fields.bathymetry_2d
+            farm = self.farm
+            self.list_powers = [0] * len(self.locations)
+            for j in range(len(self.locations)):
+                p1 = assemble(farm.turbine.power(_uv, _depth) * self.loc_dens[j] * farm.dx)
+                self.list_powers[j] = p1
+
+            return self.list_powers
+        else:
+            return self.solver_obj.fields[field_name](self.locations)
+
+    def __call__(self):
+        """
+        Evaluate all current fields in all detector locations
+
+        Returns a nturbines x ndims array, where ndims is the sum of the
+        dimension of the fields.
+        """
+        nturbines = len(self.locations)
+        field_vals = []
+        for field_name in self.field_names:
+            field_vals.append(numpy.reshape(self._evaluate_field(field_name), (nturbines, -1)))
+
+        return numpy.hstack(field_vals)