diff --git a/examples/discrete_turbines/tidal_array.py b/examples/discrete_turbines/tidal_array.py index b0d7b2e91..77a7f807d 100644 --- a/examples/discrete_turbines/tidal_array.py +++ b/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,7 +104,7 @@ 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)))) @@ -109,18 +112,49 @@ 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("") diff --git a/examples/discrete_turbines/turbine_callback.py b/examples/discrete_turbines/turbine_callback.py new file mode 100644 index 000000000..a0ab5a297 --- /dev/null +++ b/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)