From d88ab1a22360b934cd3b3e3b4b45514c9f3d1464 Mon Sep 17 00:00:00 2001 From: Stephan Kramer Date: Mon, 13 Jan 2025 19:15:29 +0000 Subject: [PATCH 1/3] DO NOT MERGE Revert to 8c6f9d1ed954a05a8d1013f79ab5b7a702c9a4a2 --- test/examples/test_examples.py | 1 - test/swe2d/test_anisotropic.py | 44 ++++++++++---------------- test_adjoint/examples/test_examples.py | 4 +-- thetis/callback.py | 25 +++------------ thetis/exporter.py | 6 ++-- thetis/utility.py | 41 ++---------------------- thetis/utility3d.py | 20 ++++++------ 7 files changed, 38 insertions(+), 103 deletions(-) diff --git a/test/examples/test_examples.py b/test/examples/test_examples.py index c9ebfe1b8..1a53238e0 100644 --- a/test/examples/test_examples.py +++ b/test/examples/test_examples.py @@ -39,7 +39,6 @@ 'tidalfarm/tidalfarm.py', 'tidal_barrage/plotting.py', 'channel_inversion/plot_elevation_progress.py', - 'channel_inversion/inverse_problem.py', 'tohoku_inversion/okada.py', 'tohoku_inversion/plot_convergence.py', 'tohoku_inversion/plot_elevation_initial_guess.py', diff --git a/test/swe2d/test_anisotropic.py b/test/swe2d/test_anisotropic.py index 41a2cd6cc..ce3ccaaa1 100644 --- a/test/swe2d/test_anisotropic.py +++ b/test/swe2d/test_anisotropic.py @@ -73,7 +73,21 @@ 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_function_spaces() + 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 boundary conditions solver_obj.bnd_functions['shallow_water'] = { @@ -121,36 +135,12 @@ 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.create_equations() + solver_obj.options.tidal_turbine_farms['everywhere'] = farm_options - # 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 and solve + # apply initial guess of inflow velocity, solve and return number of nonlinear solver iterations 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) - - # 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() # quantity of interest: power output diff --git a/test_adjoint/examples/test_examples.py b/test_adjoint/examples/test_examples.py index 18d021eea..84658fbf6 100644 --- a/test_adjoint/examples/test_examples.py +++ b/test_adjoint/examples/test_examples.py @@ -15,7 +15,7 @@ # list of all adjoint examples to run adjoint_files = [ 'tidalfarm/tidalfarm.py', - 'channel_inversion/inverse_problem.py', + # 'channel_inversion/inverse_problem.py', # FIXME requires obs time series ] cwd = os.path.abspath(os.path.dirname(__file__)) @@ -34,8 +34,6 @@ def example_file(request): def test_examples(example_file, tmpdir, monkeypatch): assert os.path.isfile(example_file), 'File not found {:}'.format(example_file) - if 'examples/channel_inversion/inverse_problem.py' in example_file: - pytest.xfail("Known issue with Firedrake and mixed function spaces. See Firedrake issue #3368.") # copy mesh files source = os.path.dirname(example_file) for f in glob.glob(os.path.join(source, '*.msh')): diff --git a/thetis/callback.py b/thetis/callback.py index 878dc65a7..52195e926 100644 --- a/thetis/callback.py +++ b/thetis/callback.py @@ -529,9 +529,6 @@ def __init__(self, solver_obj, self.field_names = field_names self._name = name - # initialise interpolation functions using vom_interpolator_functions - self.interp_functions = vom_interpolator_functions(solver_obj, field_names, detector_locations) - @property def name(self): return self._name @@ -542,8 +539,7 @@ def variable_names(self): def _values_per_field(self, values): """ - Given all values evaulated in a detector location, return the values per field - """ + Given all values evaulated in a detector location, return the values per field""" i = 0 result = [] for dim in self.field_dims: @@ -558,11 +554,7 @@ def message_str(self, *args): for name, values in zip(self.detector_names, args)) def _evaluate_field(self, field_name): - field = self.solver_obj.fields[field_name] - f_at_points, f_at_input_points = self.interp_functions[field_name] - f_at_points.interpolate(field) - f_at_input_points.interpolate(f_at_points) - return f_at_input_points.dat.data_ro[:] + return self.solver_obj.fields[field_name](self.detector_locations) def __call__(self): """ @@ -692,16 +684,10 @@ def _initialize(self): xyz = (self.x, self.y, self.z) if self.on_sphere else (self.x, self.y) self.xyz = numpy.array([xyz]) - # initialise interpolation functions using vom_interpolator_functions - self.interp_functions = vom_interpolator_functions(self.solver_obj, self.fieldnames, self.xyz) - # test evaluation try: if self.eval_func is None: - field = self.solver_obj.fields[self.fieldnames[0]] - f_at_points, f_at_input_points = self.interp_functions[self.fieldnames[0]] - f_at_points.interpolate(field) - f_at_input_points.interpolate(f_at_points) + self.solver_obj.fields.bathymetry_2d.at(self.xyz, tolerance=self.tolerance) else: self.eval_func(self.solver_obj.fields.bathymetry_2d, self.xyz, tolerance=self.tolerance) except PointNotInDomainError as e: @@ -721,10 +707,7 @@ def __call__(self): try: field = self.solver_obj.fields[fieldname] if self.eval_func is None: - f_at_points, f_at_input_points = self.interp_functions[fieldname] - f_at_points.interpolate(field) - f_at_input_points.interpolate(f_at_points) - val = f_at_input_points.dat.data_ro[:] + val = field.at(self.xyz, tolerance=self.tolerance) else: val = self.eval_func(field, self.xyz, tolerance=self.tolerance) arr = numpy.array(val) diff --git a/thetis/exporter.py b/thetis/exporter.py index fd4d1337c..5c83668c5 100644 --- a/thetis/exporter.py +++ b/thetis/exporter.py @@ -22,11 +22,11 @@ def get_visu_space(fs): """ mesh = fs.mesh() family = 'Lagrange' if is_cg(fs) else 'Discontinuous Lagrange' - if len(fs.value_shape) == 1: - dim = fs.value_shape[0] + if len(fs.ufl_element().value_shape) == 1: + dim = fs.ufl_element().value_shape[0] visu_fs = get_functionspace(mesh, family, 1, family, 1, vector=True, dim=dim) - elif len(fs.value_shape) == 2: + elif len(fs.ufl_element().value_shape) == 2: visu_fs = get_functionspace(mesh, family, 1, family, 1, tensor=True) else: diff --git a/thetis/utility.py b/thetis/utility.py index ade3a2694..d822f432e 100755 --- a/thetis/utility.py +++ b/thetis/utility.py @@ -328,7 +328,7 @@ def get_facet_mask(function_space, facet='bottom'): Here we assume that the mesh has been extruded upwards (along positive z axis). """ - from finat.element_factory import create_element as create_finat_element + from tsfc.finatinterface import create_element as create_finat_element # get base element elem = get_extruded_base_element(function_space.ufl_element()) @@ -605,8 +605,8 @@ def compute_elem_height(zcoord, output): } } }""" % {'nodes': zcoord.cell_node_map().arity, - 'func_dim': zcoord.function_space().block_size, - 'output_dim': output.function_space().block_size}, + 'func_dim': zcoord.function_space().value_size, + 'output_dim': output.function_space().value_size}, 'my_kernel') op2.par_loop( kernel, fs_out.mesh().cell_set, @@ -1154,38 +1154,3 @@ def form2indicator(F): }, ) return indicator - - -@PETSc.Log.EventDecorator("thetis.vom_interpolator_functions") -def vom_interpolator_functions(solver_obj, field_names, locations): - r""" - Creates function spaces and associated Functions for interpolation - on a VertexOnlyMesh (VOM) and returns them for reuse. - - :arg solver_obj: Thetis solver object - :arg field_names: List of field names to create functions for. - :arg locations: List of locations for interpolation. - :return: A dictionary mapping field names to a tuple of (f_at_points, f_at_input_points) - which are Functions for interpolation. - """ - vom = VertexOnlyMesh(solver_obj.mesh2d, locations, redundant=True) - - functions_dict = {} - - for field_name in field_names: - field = solver_obj.fields[field_name] - - if isinstance(field.function_space().ufl_element(), VectorElement): - P0DG = VectorFunctionSpace(vom, "DG", 0) - P0DG_input_ordering = VectorFunctionSpace(vom.input_ordering, "DG", 0) - else: - P0DG = FunctionSpace(vom, "DG", 0) - P0DG_input_ordering = FunctionSpace(vom.input_ordering, "DG", 0) - - f_at_points = Function(P0DG) - f_at_input_points = Function(P0DG_input_ordering) - - # Store the Functions in the dictionary keyed by field name - functions_dict[field_name] = (f_at_points, f_at_input_points) - - return functions_dict diff --git a/thetis/utility3d.py b/thetis/utility3d.py index d6f092415..32c8f5f07 100644 --- a/thetis/utility3d.py +++ b/thetis/utility3d.py @@ -527,8 +527,8 @@ def __init__(self, input_2d, output_3d, elem_height=None): } } }""" % {'nodes': self.fs_2d.finat_element.space_dimension(), - 'func2d_dim': self.input_2d.function_space().block_size, - 'func3d_dim': self.fs_3d.block_size, + 'func2d_dim': self.input_2d.function_space().value_size, + 'func3d_dim': self.fs_3d.value_size, 'v_nodes': n_vert_nodes}, 'my_kernel') @@ -664,8 +664,8 @@ def __init__(self, input_3d, output_2d, } } }""" % {'nodes': self.output_2d.cell_node_map().arity, - 'func2d_dim': self.output_2d.function_space().block_size, - 'func3d_dim': self.fs_3d.block_size}, + 'func2d_dim': self.output_2d.function_space().value_size, + 'func3d_dim': self.fs_3d.value_size}, 'my_kernel') else: self.kernel = op2.Kernel(""" @@ -676,8 +676,8 @@ def __init__(self, input_3d, output_2d, } } }""" % {'nodes': self.output_2d.cell_node_map().arity, - 'func2d_dim': self.output_2d.function_space().block_size, - 'func3d_dim': self.fs_3d.block_size}, + 'func2d_dim': self.output_2d.function_space().value_size, + 'func3d_dim': self.fs_3d.value_size}, 'my_kernel') if self.do_hdiv_scaling: @@ -771,8 +771,8 @@ def __init__(self, solver): } } }""" % {'nodes': self.fs_2d.finat_element.space_dimension(), - 'func2d_dim': self.fs_2d.block_size, - 'func3d_dim': self.fs_3d.block_size, + 'func2d_dim': self.fs_2d.value_size, + 'func3d_dim': self.fs_3d.value_size, 'v_nodes': n_vert_nodes}, 'my_kernel') @@ -790,8 +790,8 @@ def __init__(self, solver): } } }""" % {'nodes': self.fs_2d.finat_element.space_dimension(), - 'func2d_dim': self.fs_2d.block_size, - 'func3d_dim': self.fs_3d.block_size, + 'func2d_dim': self.fs_2d.value_size, + 'func3d_dim': self.fs_3d.value_size, 'v_nodes': n_vert_nodes}, 'my_kernel') From 042be0d25ef0e4615132b9ab9562747e395b4735 Mon Sep 17 00:00:00 2001 From: Stephan Kramer Date: Mon, 13 Jan 2025 14:00:28 +0000 Subject: [PATCH 2/3] DO NOT MERGE Rerun CI in firedrake-vanilla/2024-10 --- .github/workflows/build.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index dcfb136f4..178f909e0 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -18,7 +18,7 @@ jobs: runs-on: self-hosted # The docker container to use. container: - image: firedrakeproject/firedrake-vanilla:latest + image: firedrakeproject/firedrake-vanilla:2024-10 steps: - uses: actions/checkout@v2 - name: Cleanup From bc3212bf3288461b77fdd01feaa86f8f3c106a8f Mon Sep 17 00:00:00 2001 From: Stephan Kramer Date: Thu, 16 Jan 2025 11:24:51 +0000 Subject: [PATCH 3/3] Print timings of CI tests --- .github/workflows/build.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 178f909e0..49b0e71a4 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -35,12 +35,12 @@ jobs: run: | . /home/firedrake/firedrake/bin/activate python $(which firedrake-clean) - python -m pytest -n 12 -v test + python -m pytest --durations=0 -n 12 -v test - name: Test Thetis adjoint run: | . /home/firedrake/firedrake/bin/activate python $(which firedrake-clean) - python -m pytest -n 12 -v test_adjoint + python -m pytest --durations=0 -n 12 -v test_adjoint - name: Lint if: ${{ always() }} run: |