From c87fede41c102f77c00b7d7b63bb7d9cd2212dc6 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Mon, 6 Jan 2025 08:28:22 +0000 Subject: [PATCH 1/3] Re-enable turbines in test_anisotropic --- test/swe2d/test_anisotropic.py | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/test/swe2d/test_anisotropic.py b/test/swe2d/test_anisotropic.py index ce3ccaaa1..95b84286a 100644 --- a/test/swe2d/test_anisotropic.py +++ b/test/swe2d/test_anisotropic.py @@ -73,21 +73,7 @@ def run(solve_adjoint=False, mesh=None, **model_options): options.use_grad_depth_viscosity_term = False options.no_exports = True options.update(model_options) - solver_obj.create_equations() - - # recover Hessian - if not solve_adjoint: - if 'hessian_2d' in field_metadata: - field_metadata.pop('hessian_2d') - P1t_2d = get_functionspace(mesh2d, 'CG', 1, tensor=True) - hessian_2d = Function(P1t_2d) - u_2d = solver_obj.fields.uv_2d[0] - hessian_calculator = thetis.diagnostics.HessianRecoverer2D(u_2d, hessian_2d) - solver_obj.add_new_field(hessian_2d, - 'hessian_2d', - 'Hessian of x-velocity', - 'Hessian2d', - preproc_func=hessian_calculator) + solver_obj.create_function_spaces() # apply boundary conditions solver_obj.bnd_functions['shallow_water'] = { @@ -135,7 +121,22 @@ def bump(mesh, locs, scale=1.0): farm_options.turbine_options.diameter = D C_T = thrust_coefficient * correction farm_options.turbine_options.thrust_coefficient = C_T - solver_obj.options.tidal_turbine_farms['everywhere'] = farm_options + solver_obj.options.tidal_turbine_farms['everywhere'] = [farm_options] + solver_obj.create_equations() + + # recover Hessian + if not solve_adjoint: + if 'hessian_2d' in field_metadata: + field_metadata.pop('hessian_2d') + P1t_2d = get_functionspace(mesh2d, 'CG', 1, tensor=True) + hessian_2d = Function(P1t_2d) + u_2d = solver_obj.fields.uv_2d[0] + hessian_calculator = thetis.diagnostics.HessianRecoverer2D(u_2d, hessian_2d) + solver_obj.add_new_field(hessian_2d, + 'hessian_2d', + 'Hessian of x-velocity', + 'Hessian2d', + preproc_func=hessian_calculator) # apply initial guess of inflow velocity, solve and return number of nonlinear solver iterations solver_obj.assign_initial_conditions(uv=inflow_velocity) From ba2bf56f6465f69fbbe9430f3771654ae233cafe Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Mon, 6 Jan 2025 09:08:37 +0000 Subject: [PATCH 2/3] Add a check that turbines were picked up --- test/swe2d/test_anisotropic.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/test/swe2d/test_anisotropic.py b/test/swe2d/test_anisotropic.py index 95b84286a..319d71b1d 100644 --- a/test/swe2d/test_anisotropic.py +++ b/test/swe2d/test_anisotropic.py @@ -138,10 +138,15 @@ def bump(mesh, locs, scale=1.0): 'Hessian2d', preproc_func=hessian_calculator) - # apply initial guess of inflow velocity, solve and return number of nonlinear solver iterations + # Apply initial guess of inflow velocity and solve solver_obj.assign_initial_conditions(uv=inflow_velocity) solver_obj.iterate() + if not solve_adjoint: + # Check that turbines have been picked up + assert not np.allclose(solver_obj.fields.uv_2d.dat.data[0], 0.5, atol=1e-3) + + # Return number of nonlinear solver iterations return solver_obj.timestepper.solver.snes.getIterationNumber() # quantity of interest: power output From 50467b7d9d7cfac508d8a8ddc0c87474565d2735 Mon Sep 17 00:00:00 2001 From: Joe Wallwork Date: Mon, 6 Jan 2025 17:27:24 +0000 Subject: [PATCH 3/3] Add extra check for total_turbine_density --- test/swe2d/test_anisotropic.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/swe2d/test_anisotropic.py b/test/swe2d/test_anisotropic.py index 319d71b1d..41a2cd6cc 100644 --- a/test/swe2d/test_anisotropic.py +++ b/test/swe2d/test_anisotropic.py @@ -146,6 +146,10 @@ def bump(mesh, locs, scale=1.0): # Check that turbines have been picked up assert not np.allclose(solver_obj.fields.uv_2d.dat.data[0], 0.5, atol=1e-3) + # Check the turbine density has been set up correctly + total_turbine_density = assemble(solver_obj.tidal_farms[0].turbine_density * dx) + assert np.isclose(total_turbine_density, 2, atol=0.01), f"Expected 2, but got {total_turbine_density}" + # Return number of nonlinear solver iterations return solver_obj.timestepper.solver.snes.getIterationNumber()