Fix gradient in turbine example (#38)

Closes #37.

While #36 re-enables
the turbines in the turbine example, there's still something not right
because the optimisation progress is non-monotonic and the control
turbine area occasionally goes to zero.

I added a Taylor test and found that the gradient of the QoI wasn't
being computed correctly with the existing setup. I tracked the problem
down to the `turbine_density` expression, which includes several terms.
By projecting this expression into $\mathbb{P}1_{DG}$ space, the Taylor
test passed.

In addition, this PR overhauls the plotting functionality.

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/Makefile

@@ -25,7 +25,7 @@ go:
 	done
 
 plot:
-	python3 plot_qoi.py
+	python3 plot_results.py
 
 clean:
 	rm *.log

demos/turbine/plot_results.py

@@ -0,0 +1,74 @@
+import glob
+from warnings import warn
+
+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:
+            fname = f"data/{method}_progress_{variable}_{ext}.npy"
+            try:
+                value = np.load(fname)[-1]
+            except (FileNotFoundError, IndexError):
+                print(f"Can't load {fname}")
+                continue
+            if variable == "J":
+                descaled = -value / 10000
+                kw = descaled / 1000
+                if kw <= 0.0:
+                    warn(f"Negative power from {fname}", stacklevel=1)
+                    kw = np.nan
+                data[variable].append(kw)
+            elif variable == "m" and (value <= 0.0 or value >= 500.0):
+                warn(f"Control out of bounds from {fname}", stacklevel=1)
+                data[variable].append(np.nan)
+            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")