Fix gradient in turbine example

demos/opt_go.py

@@ -183,11 +183,11 @@ def adapt_go(mesh, target=1000.0, alpha=1.0, control=None, **kwargs):
 J = op.J_progress
 dJ = np.array([dj.dat.data[0] for dj in op.dJ_progress]).flatten()
 nc = op.nc_progress
-np.save(f"{demo}/data/go_progress_t_{n}_{method}", t)
-np.save(f"{demo}/data/go_progress_m_{n}_{method}", m)
-np.save(f"{demo}/data/go_progress_J_{n}_{method}", J)
-np.save(f"{demo}/data/go_progress_dJ_{n}_{method}", dJ)
-np.save(f"{demo}/data/go_progress_nc_{n}_{method}", nc)
+np.save(f"{demo}/data/go_progress_t_{target}_{method}", t)
+np.save(f"{demo}/data/go_progress_m_{target}_{method}", m)
+np.save(f"{demo}/data/go_progress_J_{target}_{method}", J)
+np.save(f"{demo}/data/go_progress_dJ_{target}_{method}", dJ)
+np.save(f"{demo}/data/go_progress_nc_{target}_{method}", nc)
 with open(f"{demo}/data/go_{target:.0f}_{method}.log", "w+") as f:
     note = " (FAIL)" if failed else ""
     f.write(f"cpu_time: {cpu_time}{note}\n")

demos/opt_hessian.py

@@ -51,7 +51,7 @@
         "target_inc": 0.1 * target,
         "target_max": target,
         "model_options": {
-            "no_exports": True,
+            "no_exports": not args.debug,
             "outfile": VTKFile(
                 f"{demo}/outputs_hessian/{method}/solution.pvd", adaptive=True
             ),
@@ -131,11 +131,12 @@ def adapt_hessian_based(mesh, target=1000.0, norm_order=1.0, **kwargs):
 J = op.J_progress
 dJ = np.array([dj.dat.data[0] for dj in op.dJ_progress]).flatten()
 nc = op.nc_progress
-np.save(f"{demo}/data/hessian_progress_t_{n}_{method}", t)
-np.save(f"{demo}/data/hessian_progress_m_{n}_{method}", m)
-np.save(f"{demo}/data/hessian_progress_J_{n}_{method}", J)
-np.save(f"{demo}/data/hessian_progress_dJ_{n}_{method}", dJ)
-np.save(f"{demo}/data/hessian_progress_nc_{n}_{method}", nc)
+target = int(target)
+np.save(f"{demo}/data/hessian_progress_t_{target}_{method}", t)
+np.save(f"{demo}/data/hessian_progress_m_{target}_{method}", m)
+np.save(f"{demo}/data/hessian_progress_J_{target}_{method}", J)
+np.save(f"{demo}/data/hessian_progress_dJ_{target}_{method}", dJ)
+np.save(f"{demo}/data/hessian_progress_nc_{target}_{method}", nc)
 with open(f"{demo}/data/hessian_{target:.0f}_{method}.log", "w+") as f:
     note = " (FAIL)" if failed else ""
     f.write(f"cpu_time: {cpu_time}{note}\n")

demos/turbine/plot_results.py

@@ -0,0 +1,63 @@
+import glob
+
+import matplotlib.pyplot as plt
+import numpy as np
+from thetis import create_directory
+
+create_directory("plots")
+
+data_dict = {
+    "uniform": {
+        "label": "Uniform meshing",
+    },
+    "hessian": {
+        "label": "Hessian-based",
+    },
+    # "go": {
+    #     "label": "Goal-oriented",
+    # },
+}
+
+# Load data from files
+variables = ("nc", "J", "t", "m")
+for method, data in data_dict.items():
+    for variable in variables:
+        data[variable] = []
+    for fname in glob.glob(f"data/{method}_*.log"):
+        ext = "_".join(fname.split("_")[1:])[:-4]
+        for variable in variables:
+            try:
+                value = np.load(f"data/{method}_progress_{variable}_{ext}.npy")[-1]
+            except (FileNotFoundError, IndexError):
+                continue
+            if variable == "J":
+                data[variable].append(-value / 1000)
+            else:
+                data[variable].append(value)
+
+metadata = {
+    "J": {"label": "qoi", "name": r"Power output ($\mathrm{kW}$)"},
+    "nc": {"label": "elements", "name": "Number of mesh elements"},
+    "t": {"label": "time", "name": r"CPU time ($\mathrm{s}$)"},
+    "m": {"label": "control", "name": r"Control ($\mathrm{m}$)"},
+}
+
+
+def plot(v1, v2):
+    fig, axes = plt.subplots(figsize=(5, 3))
+    for data in data_dict.values():
+        axes.semilogx(data[v1], data[v2], "x", label=data["label"])
+    axes.set_xlabel(metadata[v1]["name"])
+    axes.set_ylabel(metadata[v2]["name"])
+    axes.grid(True, which="both")
+    axes.legend()
+    plt.tight_layout()
+    fname = f"converged_{metadata[v2]['label']}_vs_{metadata[v1]['label']}"
+    for ext in ("pdf", "jpg"):
+        plt.savefig(f"plots/{fname}.{ext}")
+
+
+plot("nc", "J")
+plot("t", "J")
+plot("nc", "m")
+plot("t", "m")

demos/turbine/setup.py

@@ -3,6 +3,7 @@
 import numpy as np
 import ufl
 from animate.metric import RiemannianMetric
+from firedrake.adjoint import pyadjoint
 from firedrake.assemble import assemble
 from firedrake.constant import Constant
 from firedrake.function import Function
@@ -24,16 +25,18 @@ def initial_mesh(n=4):
 
 def initial_control(mesh):
     R = FunctionSpace(mesh, "R", 0)
-    return Function(R).assign(250.0)
+    return Function(R, val=250.0)
 
 
-def forward_run(mesh, control=None, outfile=None, **model_options):
+def forward_run(mesh, control=None, outfile=None, debug=False, **model_options):
     """
     Solve the shallow water flow-past-a-turbine problem on a given mesh.
 
     Optionally, pass an initial value for the control variable
     (y-coordinate of the centre of the second turbine).
     """
+    if debug:
+        pyadjoint.continue_annotation()
     x, y = ufl.SpatialCoordinate(mesh)
 
     # Setup bathymetry
@@ -72,8 +75,7 @@ def forward_run(mesh, control=None, outfile=None, **model_options):
     # Setup boundary conditions
     P1v_2d = solver_obj.function_spaces.P1v_2d
     u_in = Function(P1v_2d)
-    u_in.dat.data[:, 0] = 5.0
-    u_in.dat.data[:, 1] = 0.0
+    u_in.interpolate(ufl.as_vector([5.0, 0.0]))
     bcs = {
         1: {"uv": u_in},
         2: {"elev": Constant(0.0)},
@@ -87,14 +89,14 @@ def forward_run(mesh, control=None, outfile=None, **model_options):
     R = FunctionSpace(mesh, "R", 0)
     ym = 250.0
     sep = 60.0
-    x1 = Constant(456.0)
-    x2 = Constant(456.0)
-    x3 = Constant(456.0)
-    xc = Constant(744.0)
-    y1 = Constant(ym)
-    y2 = Constant(ym + sep)
-    y3 = Constant(ym - sep)
-    yc = Function(R).assign(control or 250.0)
+    x1 = Function(R, val=456.0)
+    x2 = Function(R, val=456.0)
+    x3 = Function(R, val=456.0)
+    xc = Function(R, val=744.0)
+    y1 = Function(R, val=ym)
+    y2 = Function(R, val=ym + sep)
+    y3 = Function(R, val=ym - sep)
+    yc = Function(R, val=control or 250.0)
     thrust_coefficient = 0.8
     turbine_diameter = 18.0
     turbine_footprint = turbine_diameter**2
@@ -103,14 +105,18 @@ def bump(x0, y0, label):
         r = turbine_diameter / 2
         qx = ((x - x0) / r) ** 2
         qy = ((y - y0) / r) ** 2
-        cond = ufl.And(qx < 1, qy < 1)
-        b = ufl.exp(1 - 1 / (1 - qx)) * ufl.exp(1 - 1 / (1 - qy))
-        cond = ufl.conditional(cond, Constant(1.0) * b, 0)
+        cond = ufl.conditional(
+            ufl.And(qx < 1, qy < 1),
+            ufl.exp(1 - 1 / (1 - qx)) * ufl.exp(1 - 1 / (1 - qy)),
+            0,
+        )
         integral = assemble(cond * ufl.dx)
         assert integral > 0.0, f"Invalid area for {label}"
         return cond / integral
 
-    turbine_density = (
+    P1DG_2d = solver_obj.function_spaces.P1DG_2d
+    turbine_density = Function(P1DG_2d)
+    turbine_density.project(
         bump(x1, y1, "turbine 1")
         + bump(x2, y2, "turbine 2")
         + bump(x3, y3, "turbine 3")
@@ -143,13 +149,13 @@ def bump(x0, y0, label):
     u, eta = ufl.split(solver_obj.fields.solution_2d)
     swiped_area = (ufl.pi * turbine_diameter / 2) ** 2
     area_frac = swiped_area / turbine_footprint
-    coeff = -rho * 0.5 * thrust_coefficient * area_frac * turbine_density
-    J_power = coeff * ufl.dot(u, u) ** 1.5 * ufl.dx
+    coeff = rho * 0.5 * thrust_coefficient * area_frac * turbine_density
+    J_power = -coeff * ufl.dot(u, u) ** 1.5 * ufl.dx
     # NOTE: negative because we want maximum
 
     # Add a regularisation term for constraining the control
     area = assemble(Constant(1.0) * ufl.dx(domain=mesh))
-    alpha = 1.0 / area
+    alpha = Constant(1.0 / area)
     J_reg = (
         alpha
         * ufl.conditional(
@@ -159,6 +165,17 @@ def bump(x0, y0, label):
     )
 
     J = assemble(J_power + J_reg, ad_block_tag="qoi")
+
+    if debug:
+        controls = {"q_2d": solver_obj.fields.solution_2d}
+        for key, control in controls.items():
+            Jhat = pyadjoint.ReducedFunctional(J, pyadjoint.Control(control))
+            h = Function(control)
+            h.assign(0.1)
+            assert pyadjoint.taylor_test(Jhat, control, h) > 1.9
+            print(f"Taylor test for {key} passed")
+        pyadjoint.pause_annotation()
+
     return J, yc
 
 
@@ -202,4 +219,4 @@ def hessian_component(f):
     resolution = 4
     init_mesh = initial_mesh(n=resolution)
     init_control = initial_control(init_mesh)
-    forward_run(init_mesh, init_control, outfile=VTKFile("test.pvd"))
+    forward_run(init_mesh, init_control, outfile=VTKFile("test.pvd"), debug=True)