diff --git a/demos/opt_go.py b/demos/opt_go.py index fec1a10..8a528bd 100644 --- a/demos/opt_go.py +++ b/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") diff --git a/demos/opt_hessian.py b/demos/opt_hessian.py index 49cf102..e5976f5 100644 --- a/demos/opt_hessian.py +++ b/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") diff --git a/demos/turbine/Makefile b/demos/turbine/Makefile index e5469b8..1818619 100644 --- a/demos/turbine/Makefile +++ b/demos/turbine/Makefile @@ -25,7 +25,7 @@ go: done plot: - python3 plot_qoi.py + python3 plot_results.py clean: rm *.log diff --git a/demos/turbine/plot_qoi.py b/demos/turbine/plot_qoi.py deleted file mode 100644 index d16a5aa..0000000 --- a/demos/turbine/plot_qoi.py +++ /dev/null @@ -1,48 +0,0 @@ -import glob - -import matplotlib.pyplot as plt -import numpy as np -from thetis import create_directory - - -def plot_progress(axes, method, ext, gradients=True, **kwargs): - """ - Plot the progress of an optimisation run, in terms - of control value vs. objective functional value. - - :arg axes: the figure axes to plot on - :arg method: the optimisation method - :kwargs: to be passed to matplotlib's plotting functions - """ - m = np.load(f"data/{method}_progress_m_{ext}.npy")[-1] - off = np.abs(m - 250.0) - J = -np.load(f"data/{method}_progress_J_{ext}.npy")[-1] / 1000 - axes.plot(off, J, "x", **kwargs) - - -fig, axes = plt.subplots(figsize=(5, 3)) -labels = { - "uniform": "Uniform meshing", - "hessian": "Hessian-based", - "go": "Goal-based", -} -for method, label in labels.items(): - cpu = [] - J = [] - for fname in glob.glob(f"data/{method}_*.log"): - ext = "_".join(fname.split("_")[1:])[:-4] - with open(fname, "r") as f: - words = f.readline().split() - if len(words) > 2 and "FAIL" in words[2]: - continue - cpu.append(float(words[1])) - J.append(-np.load(f"data/{method}_progress_J_{ext}.npy")[-1] / 1000) - axes.semilogx(cpu, J, "x", label=label) -axes.set_xlabel(r"CPU time ($\mathrm{s}$)") -axes.set_ylabel(r"Power output ($\mathrm{kW}$)") -axes.grid(True, which="both") -axes.legend() -plt.tight_layout() -create_directory("plots") -plt.savefig("plots/converged_qoi.pdf") -plt.savefig("plots/converged_qoi.jpg") diff --git a/demos/turbine/plot_results.py b/demos/turbine/plot_results.py new file mode 100644 index 0000000..303f28e --- /dev/null +++ b/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") diff --git a/demos/turbine/setup.py b/demos/turbine/setup.py index 486db88..fbf832d 100644 --- a/demos/turbine/setup.py +++ b/demos/turbine/setup.py @@ -3,14 +3,16 @@ 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 -from firedrake.functionspace import FunctionSpace, TensorFunctionSpace +from firedrake.functionspace import TensorFunctionSpace from firedrake.utility_meshes import RectangleMesh -from thetis.options import TidalTurbineFarmOptions +from thetis.options import DiscreteTidalTurbineFarmOptions from thetis.solver2d import FlowSolver2d -from thetis.utility import get_functionspace +from thetis.turbines import TurbineFunctionalCallback +from thetis.utility import domain_constant from opt_adapt.opt import get_state @@ -23,11 +25,10 @@ def initial_mesh(n=4): def initial_control(mesh): - R = FunctionSpace(mesh, "R", 0) - return Function(R).assign(250.0) + return domain_constant(260.0, mesh) -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. @@ -36,130 +37,97 @@ def forward_run(mesh, control=None, outfile=None, **model_options): """ x, y = ufl.SpatialCoordinate(mesh) - # Setup bathymetry - H = 40.0 - P1_2d = get_functionspace(mesh, "CG", 1) - bathymetry = Function(P1_2d) - bathymetry.assign(H) + # Specify bathymetry + channel_depth = 40.0 # Setup solver - solver_obj = FlowSolver2d(mesh, bathymetry) + solver_obj = FlowSolver2d(mesh, Constant(channel_depth)) options = solver_obj.options - options.element_family = "dg-cg" - options.timestep = 20.0 - options.simulation_export_time = 20.0 - options.simulation_end_time = 18.0 + options.element_family = "dg-dg" + options.timestep = 1.0 + options.simulation_export_time = 1.0 + options.simulation_end_time = 0.5 options.swe_timestepper_type = "SteadyState" options.swe_timestepper_options.solver_parameters = { - "mat_type": "aij", - "snes_type": "newtonls", - "snes_linesearch_type": "bt", - "snes_rtol": 1.0e-08, - "snes_max_it": 100, - "ksp_type": "preonly", - "pc_type": "lu", - "pc_factor_mat_solver_type": "mumps", + "snes_rtol": 1.0e-12, } - options.use_grad_div_viscosity_term = False + # options.use_grad_div_viscosity_term = False options.horizontal_viscosity = Constant(0.5) options.quadratic_drag_coefficient = Constant(0.0025) - options.use_lax_friedrichs_velocity = True - options.lax_friedrichs_velocity_scaling_factor = Constant(1.0) - options.use_grad_depth_viscosity_term = False + # options.use_grad_depth_viscosity_term = False options.update(model_options) - solver_obj.create_function_spaces() # 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 - bcs = { - 1: {"uv": u_in}, + solver_obj.bnd_functions["shallow_water"] = { + 1: {"uv": Constant((3.0, 0.0))}, 2: {"elev": Constant(0.0)}, 3: {"un": Constant(0.0)}, + 4: {"un": Constant(0.0)}, } - if 4 in mesh.exterior_facets.unique_markers: - bcs[4] = {"un": Constant(0.0)} - solver_obj.bnd_functions["shallow_water"] = bcs - - # Define turbine parameters - 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) - thrust_coefficient = 0.8 - turbine_diameter = 18.0 - turbine_footprint = turbine_diameter**2 - - 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) - integral = assemble(cond * ufl.dx) - assert integral > 0.0, f"Invalid area for {label}" - return cond / integral - - turbine_density = ( - bump(x1, y1, "turbine 1") - + bump(x2, y2, "turbine 2") - + bump(x3, y3, "turbine 3") - + bump(xc, yc, "control turbine") - ) - - # Apply thrust correction - vertical_slice = H * turbine_diameter - thrust_coefficient *= 4.0 / ( - 1.0 + ufl.sqrt(1.0 - thrust_coefficient * turbine_footprint / vertical_slice) - ) + solver_obj.create_function_spaces() # Setup tidal farm - farm_options = TidalTurbineFarmOptions() + farm_options = DiscreteTidalTurbineFarmOptions() + turbine_density = Function(solver_obj.function_spaces.P1_2d).assign(1.0) + farm_options.turbine_type = "constant" farm_options.turbine_density = turbine_density - farm_options.turbine_options.diameter = turbine_diameter - farm_options.turbine_options.thrust_coefficient = thrust_coefficient - options.tidal_turbine_farms["everywhere"] = [farm_options] - solver_obj.create_equations() - rho = Constant(1030.0) + farm_options.turbine_options.diameter = 18.0 + farm_options.turbine_options.thrust_coefficient = 0.8 + farm_options.quadrature_degree = 100 + farm_options.upwind_correction = False + farm_options.turbine_coordinates = [ + [domain_constant(x, mesh=mesh), domain_constant(y, mesh=mesh)] + for x, y in [[456, 250], [456, 310], [456, 190], [744, control or 250]] + ] + y2 = farm_options.turbine_coordinates[1][1] + y3 = farm_options.turbine_coordinates[2][1] + yc = farm_options.turbine_coordinates[3][1] + options.discrete_tidal_turbine_farms["everywhere"] = [farm_options] + + # Add a callback for computing the power output + cb = TurbineFunctionalCallback(solver_obj) + solver_obj.add_callback(cb, "timestep") # Apply initial conditions and solve - solver_obj.assign_initial_conditions(uv=u_in) + solver_obj.assign_initial_conditions(uv=(ufl.as_vector((1.0e-03, 0.0)))) solver_obj.iterate() if outfile is not None: u, eta = solver_obj.fields.solution_2d.subfunctions outfile.write(u, eta) - # Define objective function - 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 - # NOTE: negative because we want maximum + J_power = sum(cb.integrated_power) # Add a regularisation term for constraining the control - area = assemble(Constant(1.0) * ufl.dx(domain=mesh)) - alpha = 1.0 / area - J_reg = ( + area = assemble(domain_constant(1.0, mesh) * ufl.dx) + alpha = domain_constant(1.0 / area, mesh) + J_reg = assemble( alpha * ufl.conditional( - yc < y2, (yc - y2) ** 2, ufl.conditional(yc > y3, (yc - y3) ** 2, 0) + yc < y3, (yc - y3) ** 2, ufl.conditional(yc > y2, (yc - y2) ** 2, 0) ) * ufl.dx ) - J = assemble(J_power + J_reg, ad_block_tag="qoi") - return J, yc + # Sum the two components + # NOTE: We rescale the functional such that the gradients are ~ order magnitude 1 + # NOTE: We also multiply by -1 so that if we minimise the functional, we maximise + # power (maximize is also available from pyadjoint but currently broken) + scaling = 10000 + J = scaling * (-J_power + J_reg) + + print(f"DEBUG power: {J_power:.4e}, reg: {J_reg:.4e}") + + control_variable = yc + if debug: + # Perform a Taylor test + Jhat = pyadjoint.ReducedFunctional(J, pyadjoint.Control(control_variable)) + h = Function(control_variable) + np.random.seed(23) + h.dat.data[:] = np.random.random(h.dat.data.shape) + assert pyadjoint.taylor_test(Jhat, control, h) > 1.95 + print("Taylor test passed") + + return J, control_variable def hessian(mesh, **kwargs): @@ -199,7 +167,9 @@ def hessian_component(f): if __name__ == "__main__": from firedrake.output.vtk_output import VTKFile + pyadjoint.continue_annotation() 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) + pyadjoint.pause_annotation() diff --git a/opt_adapt/opt.py b/opt_adapt/opt.py index 1ee538b..c509bb8 100644 --- a/opt_adapt/opt.py +++ b/opt_adapt/opt.py @@ -1,4 +1,5 @@ from time import perf_counter +from warnings import warn import firedrake as fd import numpy as np @@ -217,6 +218,12 @@ def _gradient_descent(it, forward_run, m, params, u, u_, dJ_): dJ_diff = fd.assemble(ufl.inner(dJ_ - dJ, dJ_ - dJ) * ufl.dx) lr = abs(fd.assemble(ufl.inner(u_ - u, dJ_ - dJ) * ufl.dx) / dJ_diff) lr = max(lr, params.lr_min) + if np.isnan(lr) or lr <= 0.0: + warn( + "Barzilai-Borwein formula gave an invalid learning rate, using default", + stacklevel=1, + ) + lr = params.lr # Take a step downhill u += lr * P