Merge branch 'master' into correct-adjoint-scaling

examples/discrete_turbines/channel-optimisation.py

@@ -8,7 +8,7 @@
 import random
 op2.init(log_level=INFO)
 
-output_dir = 'outputs'
+output_dir = 'outputs_optimisation'
 
 if os.getenv('THETIS_REGRESSION_TEST') is not None:
     test_gradient = True  # test gradient using Taylor test (see below)

examples/discrete_turbines/tidal_array.py

@@ -1,14 +1,15 @@
 """
 Basic example of a discrete turbine array placed in the high energy region
-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.
+of flow past a headland. Turbines are based on the SIMEC Atlantis AR1500 and
+AR2000 models with cut-in, rated and cut-out speeds.
 Flow becomes steady after an initial ramp up.
-Two callbacks are used - one for the farm and one for discrete turbines.
+Two callbacks are used - one for the overall farm and one for discrete turbines.
 """
 
 from thetis import *
 from firedrake.output.vtk_output import VTKFile
 from turbine_callback import TidalPowerCallback
+from copy import deepcopy
 
 
 # Set output directory, load mesh, set simulation export and end times
@@ -17,8 +18,8 @@
     os.system('gmsh -2 headland.geo -o headland.msh')
 mesh2d = Mesh('headland.msh')
 site_ID = 2  # mesh PhysID for subdomain where turbines are to be sited
-print_output('Loaded mesh ' + mesh2d.name)
-print_output('Exporting to ' + outputdir)
+print_output(f'Loaded mesh {mesh2d.name}')
+print_output(f'Exporting to {outputdir}')
 
 t_end = 1.5 * 3600
 t_export = 200.0
@@ -41,35 +42,63 @@
 
 # Define the thrust curve of the turbine using a tabulated approach:
 # speeds_AR2000: speeds for corresponding thrust coefficients - thrusts_AR2000
+# powers_AR2000: list of idealised power coefficients of an AR2000 tidal turbine
 # 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]
+powers_AR2000 = [0.0105, 0.032, 0.0385, 0.116, 0.437, 0.437, 0.345, 0.277, 0.226, 0.187, 0.156, 0.132, 0.112, 0.0993,
+                 0.0595, 0.0051, 0.00151, 0.000889, 0.000652, 0.000523, 0.000441, 0.000384, 0.000341, 0.000308]
 thrusts_AR2000 = [0.010531, 0.032281, 0.038951, 0.119951, 0.516484, 0.516484, 0.387856, 0.302601, 0.242037, 0.197252,
                   0.16319, 0.136716, 0.115775, 0.102048, 0.060513, 0.005112, 0.00151, 0.00089, 0.000653, 0.000524,
                   0.000442, 0.000384, 0.000341, 0.000308]
 
+# Set the water density to match the thrust and power curves (defaults to 1000kg/m3)
+physical_constants['rho0'] = 1026.0
+
 # initialise discrete turbine farm characteristics
-farm_options = DiscreteTidalTurbineFarmOptions()
-farm_options.turbine_type = turbine_thrust_def
+farm_options_AR2000 = DiscreteTidalTurbineFarmOptions()
+farm_options_AR2000.turbine_type = turbine_thrust_def
 if turbine_thrust_def == 'table':
-    farm_options.turbine_options.thrust_speeds = speeds_AR2000
-    farm_options.turbine_options.thrust_coefficients = thrusts_AR2000
+    farm_options_AR2000.turbine_options.thrust_speeds = speeds_AR2000
+    farm_options_AR2000.turbine_options.thrust_coefficients = thrusts_AR2000
+    farm_options_AR2000.turbine_options.power_coefficients = powers_AR2000
 else:
-    farm_options.turbine_options.thrust_coefficient = 0.6
+    farm_options_AR2000.turbine_options.thrust_coefficient = 0.6
+    farm_options_AR2000.turbine_options.power_coefficient = 0.55
 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.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)
-farm_options.turbine_density = turbine_density
-farm_options.turbine_coordinates = [[Constant(x), Constant(y)]
-                                    for x in numpy.arange(940, 1061, 60)
-                                    for y in numpy.arange(260, 341, 40)]
+    farm_options_AR2000.turbine_options.C_support = 0.7  # support structure thrust coefficient
+    farm_options_AR2000.turbine_options.A_support = 2.6 * 14.0  # cross-sectional area of support structure
+farm_options_AR2000.turbine_options.diameter = 20
+farm_options_AR2000.upwind_correction = True  # See https://arxiv.org/abs/1506.03611 for more details
+turbine_density = Function(FunctionSpace(mesh2d, "CG", 1), name='turbine_density_AR2000').assign(0.0)
+farm_options_AR2000.turbine_density = turbine_density
+farm_options_AR2000.turbine_coordinates = [[Constant(x), Constant(y)]
+                                           for x in numpy.arange(1000, 1061, 60)
+                                           for y in numpy.arange(260, 341, 40)]
+
+# Now create a second farm with AR1500s
+speeds_AR1500 = speeds_AR2000.copy()
+powers_AR1500 = [0.00953, 0.0291, 0.035, 0.106, 0.405, 0.405, 0.32, 0.257, 0.209, 0.173, 0.145, 0.122, 0.104, 0.0919,
+                 0.0551, 0.00471, 0.00139, 0.000821, 0.000602, 0.000483, 0.000408, 0.000355, 0.000315, 0.000285]
+thrusts_AR1500 = [0.00955, 0.0293, 0.0353, 0.109, 0.468, 0.468, 0.355, 0.278, 0.223, 0.182, 0.15, 0.126, 0.107, 0.0942,
+                  0.0559, 0.00472, 0.00139, 0.000821, 0.000602, 0.000483, 0.000408, 0.000355, 0.000315, 0.000285]
+farm_options_AR1500 = deepcopy(farm_options_AR2000)
+if turbine_thrust_def == 'table':
+    farm_options_AR1500.turbine_options.thrust_speeds = speeds_AR1500
+    farm_options_AR1500.turbine_options.thrust_coefficients = thrusts_AR1500
+    farm_options_AR1500.turbine_options.power_coefficients = powers_AR1500
+else:
+    farm_options_AR1500.turbine_options.thrust_coefficient = 0.6
+    farm_options_AR1500.turbine_options.power_coefficient = 0.55
+farm_options_AR1500.turbine_options.diameter = 18
+turbine_density = Function(FunctionSpace(mesh2d, "CG", 1), name='turbine_density_AR1500').assign(0.0)
+farm_options_AR1500.turbine_density = turbine_density
+farm_options_AR1500.turbine_coordinates = [[Constant(940), Constant(y)] for y in numpy.arange(260, 341, 40)]
+
 
 # --- create solver ---
 solver_obj = solver2d.FlowSolver2d(mesh2d, bathymetry_2d)
@@ -87,7 +116,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_AR1500, farm_options_AR2000]
 
 # Use direct solver instead of default iterative settings
 # (see SemiImplicitSWETimeStepperOptions2d in thetis/options.py)
@@ -109,52 +138,72 @@
 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(f'{options.swe_timestepper_type} solver options:')
 print_output(options.swe_timestepper_options.solver_parameters)
 
 # 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
+# export the turbine density
+turbine_density_function = Function(P1_2d, name="Turbine Density")
+turbine_density_function.project(solver_obj.tidal_farms[0].turbine_density)
+VTKFile(outputdir + '/turbine_density_AR1500.pvd').write(turbine_density_function)
+turbine_density_function.project(solver_obj.tidal_farms[1].turbine_density)
+VTKFile(outputdir + '/turbine_density_AR2000.pvd').write(turbine_density_function)
 
 # 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')
+num_AR1500s = len(farm_options_AR1500.turbine_coordinates)
+turbine_names = [f'AR1500_{i+1}' for i in range(num_AR1500s)]
+cb_turbines_AR1500 = TidalPowerCallback(solver_obj, site_ID, farm_options_AR1500, ['pow'], name='AR1500array',
+                                        turbine_names=turbine_names)
+solver_obj.add_callback(cb_turbines_AR1500, 'timestep')
+turbine_names = [f'AR2000_{i+1}' for i in range(len(farm_options_AR2000.turbine_coordinates))]
+cb_turbines_AR2000 = TidalPowerCallback(solver_obj, site_ID, farm_options_AR2000, ['pow'], name='AR2000array',
+                                        turbine_names=turbine_names)
+solver_obj.add_callback(cb_turbines_AR2000, 'timestep')
 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]) + ')')
+for i, name in enumerate(cb_turbines_AR1500.variable_names):
+    print_output(f'{name}: ' f'({cb_turbines_AR1500.get_turbine_coordinates[i][0]}, '
+                 f'{cb_turbines_AR1500.get_turbine_coordinates[i][1]})')
+for i, name in enumerate(cb_turbines_AR2000.variable_names):
+    print_output(f'{name}: ' f'({cb_turbines_AR2000.get_turbine_coordinates[i][0]}, '
+                 f'{cb_turbines_AR2000.get_turbine_coordinates[i][1]})')
 print_output("")
 
 
 def update_forcings(t_new):
     ramp = tanh(t_new / 2000.)
     tidal_vel.project(Constant(ramp * 3.))
-    power_farm.append(cb_farm.instantaneous_power[0])
-    powers_turbines.append(cb_turbines())
+    power_farm.append(cb_farm.instantaneous_power.copy())
+    powers_turbines.append(numpy.concat((cb_turbines_AR1500(), cb_turbines_AR2000())))
 
 
 # See channel-optimisation example for a completely steady state simulation (no ramp)
 solver_obj.iterate(update_forcings=update_forcings)
 
-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!
+power_farm.append(cb_farm.instantaneous_power)  # add final powers, should be the same as callback hdf5 file!
+power_farm = np.array(power_farm).T
+AR1500farm_energy = np.sum(power_farm[0]) * options.timestep / 3600
+AR2000farm_energy = np.sum(power_farm[1]) * options.timestep / 3600
+farm_energy = AR1500farm_energy + AR2000farm_energy
+powers_turbines.append(np.concat((cb_turbines_AR1500(), cb_turbines_AR2000())))
 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')
+for i, name in enumerate(cb_turbines_AR1500.variable_names):
+    print_output(f'{name}: {turbine_energies[i][0]*10**(-3):.2f}kWh')
+for i, name in enumerate(cb_turbines_AR2000.variable_names):
+    print_output(f'{name}: {turbine_energies[num_AR1500s+i][0]*10**(-3):.2f}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(f"Farm callback total energy recorded: {farm_energy*10**(-3):.2f}kWh")
+print_output(f"Turbines callback total energy recorded: {np.sum(turbine_energies*10**(-3), axis=0)[0]:.2f}kWh")
+print_output(f"Difference: {callback_diff[0]:.5f}%")
 print_output("")

examples/tidalfarm/tidalfarm.py

@@ -102,11 +102,11 @@ def update_forcings(t):
 # - farm_options.turbine_options.diameter (default 16.0)
 farm_options = TidalTurbineFarmOptions()
 farm_options.turbine_density = turbine_density
-# amount of power produced per turbine (kW) on average to "break even" (cost = revenue)
+# amount of power produced per turbine (W) on average to "break even" (cost = revenue)
 # this is used to scale the cost, which is assumed to be linear with the number of turbines,
-# in such a way that the cost is expressed in kW which can be subtracted from the profit
+# in such a way that the cost is expressed in W which can be subtracted from the profit
 # which is calculated as the power extracted by the turbines
-farm_options.break_even_wattage = 200
+farm_options.break_even_wattage = 200000
 options.tidal_turbine_farms[2] = [farm_options]
 
 # we first run the "forward" model with no turbines
@@ -144,7 +144,7 @@ def update_forcings(t):
 # the scaling is based on the maximum cost term
 # also we multiply by -1 so that if we minimize the functional, we maximize profit
 # (maximize is also availble from pyadjoint but currently broken)
-scaling = -1./assemble(max(farm_options.break_even_wattage, 100) * max_density * dx(2, domain=mesh2d))
+scaling = -1./assemble(max(farm_options.break_even_wattage, 100000) * max_density * dx(2, domain=mesh2d))
 scaled_functional = scaling * cb.average_profit[0]
 
 # specifies the control we want to vary in the optimisation

thetis/coupled_timeintegrator_2d.py

@@ -199,11 +199,12 @@ def initialize(self, solution2d):
         self.timesteppers.swe2d.initialize(self.fields.solution_2d)
         if self.nh_options.update_free_surface:
             self.timesteppers.fs2d.initialize(self.fields.elev_2d)
-            self.elev_old.assign(self.fields.elev_2d)
 
     @PETSc.Log.EventDecorator("thetis.NonHydrostaticTimeIntegrator2D.advance")
     def advance(self, t, update_forcings=None):
         """Advances equations for one time step."""
+        if self.nh_options.update_free_surface:
+            self.elev_old.assign(self.fields.elev_2d)
         if self.serial_advancing:
             # --- advance in serial ---
             self.timesteppers.swe2d.advance(t, update_forcings=update_forcings)
@@ -213,7 +214,6 @@ def advance(self, t, update_forcings=None):
             if self.nh_options.update_free_surface:
                 self.fields.elev_2d.assign(self.elev_old)
                 self.timesteppers.fs2d.advance(t, update_forcings=update_forcings)
-                self.elev_old.assign(self.fields.elev_2d)
             # update old solution
             if self.options.swe_timestepper_type == 'SSPIMEX':
                 self.timesteppers.swe2d.erk.solution_old.assign(self.fields.solution_2d)
@@ -234,4 +234,3 @@ def advance(self, t, update_forcings=None):
                     elif last_stage:
                         self.fields.elev_2d.assign(self.elev_old)
                         self.timesteppers.fs2d.advance(t, update_forcings=update_forcings)
-                        self.elev_old.assign(self.fields.elev_2d)

thetis/options.py

@@ -473,16 +473,24 @@ class ConstantTidalTurbineOptions(TidalTurbineOptions):
     name = 'Constant tidal turbine options'
     thrust_coefficient = PositiveFloat(
         0.8, help='Thrust coefficient C_T').tag(config=True)
+    power_coefficient = PositiveFloat(
+        None, allow_none=True, help='Power coefficient C_P, defaults to None and will be calculated using '
+                                    'an empirical formula based on thrust coefficients unless specified (see '
+                                    'ConstantThrustTurbine).').tag(config=True)
 
 
 class TabulatedTidalTurbineOptions(TidalTurbineOptions):
     """Options for tidal turbine with tabulated thrust coefficient"""
     name = 'Tabulated tidal turbine options'
-    thrust_coefficients = List([3.0], help='Table of thrust coefficients')
     thrust_speeds = List(
-        [0.8, 0.8],
-        help="""List of speeds at which thrust_coefficients are applied.
-        First and last entry function as cut-in and cut-out speeds respectively""")
+        [0.9, 1., 3., 5., 5.001],
+        help='List of speeds at which thrust_coefficients are applied.'
+        'First and last entry function as cut-in and cut-out speeds respectively')
+    thrust_coefficients = List([0.01, 0.7, 0.7, 0.1, 0.0001], help='Table of thrust coefficients')
+    power_coefficients = List(
+        None, allow_none=True, help='Table of power coefficients. Defaults to None and will be calculated '
+                                    'using an empirical formula based on thrust coefficients unless specified (see '
+                                    'TabulatedThrustTurbine).').tag(config=True)
 
 
 @attach_paired_options("turbine_type",

thetis/timeintegrator.py

@@ -503,8 +503,8 @@ def __init__(self, equation, solution, fields, dt, options, bnd_conditions, term
 
         fs = self.equation.function_space
         self.solution_old = Function(fs, name='old solution')
-        self.msolution_old = Function(fs, name='dual solution')
-        self.rhs_func = Function(fs, name='rhs linear form')
+        self.msolution_old = Cofunction(fs.dual(), name='dual solution')
+        self.rhs_func = Cofunction(fs.dual(), name='rhs linear form')
 
         # fully explicit evaluation
         self.a = self.equation.mass_term(self.equation.trial)
@@ -639,9 +639,9 @@ def __init__(self, equation, solution, fields, dt, options, bnd_conditions, term
         super(SSPRK22ALE, self).__init__(equation, solution, fields, dt, options)
 
         fs = self.equation.function_space
-        self.mu = Function(fs, name='dual solution')
-        self.mu_old = Function(fs, name='dual solution')
-        self.tendency = Function(fs, name='tendency')
+        self.mu = Cofunction(fs.dual(), name='dual solution')
+        self.mu_old = Cofunction(fs.dual(), name='dual solution')
+        self.tendency = Cofunction(fs.dual(), name='tendency')
 
         # fully explicit evaluation
         self.a = self.equation.mass_term(self.equation.trial)