From 5a086368f99d61b9f984bc731df6b560b7729a65 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 29 Jun 2025 22:00:49 -0600 Subject: [PATCH 1/6] Implement an interface that external surface deformation tools can use. --- .../external_tool_interface.h | 191 ++++++++++++++++++ .../external_tool_interface.cc | 159 +++++++++++++++ 2 files changed, 350 insertions(+) create mode 100644 include/aspect/mesh_deformation/external_tool_interface.h create mode 100644 source/mesh_deformation/external_tool_interface.cc diff --git a/include/aspect/mesh_deformation/external_tool_interface.h b/include/aspect/mesh_deformation/external_tool_interface.h new file mode 100644 index 00000000000..8133ec5c085 --- /dev/null +++ b/include/aspect/mesh_deformation/external_tool_interface.h @@ -0,0 +1,191 @@ +/* + Copyright (C) 2011 - 2024 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + + +#ifndef _aspect_mesh_deformation_external_tool_interface_h +#define _aspect_mesh_deformation_external_tool_interface_h + +#include + + +namespace aspect +{ + namespace MeshDeformation + { + /** + * This class provides some support for writing classes derived + * from MeshDeformation::Interface when implementing surface + * deformation models that are based on external tools such as + * Fastscape, Landlab, OpenLEM, etc. The primary complication of + * interacting with these external tools is that they generally + * use their own discretization of surface processes, using a mesh + * that, in most cases, will be different from the one used by + * ASPECT. (Of course, ASPECT's use is a *volume* mesh, whereas + * surface processes use *surface* meshes. ASPECT's own surface + * diffusion, for example, works on the surface faces of ASPECT's + * mesh, but in general, external tools use a mesh that is + * *entirely unrelated* to the surface cells of ASPECT's volume + * mesh.) In these cases, one needs to implement transfer + * operations to take ASPECT's solution and interpolate it to the + * external tool's mesh, run some surface evolution steps, and + * then interpolate back to the surface nodes of ASPECT's + * mesh. These interpolation operations are difficult to implement + * in their own right, but are made particularly cumbersome in + * parallel because then, the points at which the external tool + * wants to know the solution on any given MPI process will, in + * general, not be located on that part of ASPECT's mesh that is + * owned by that MPI process. As a consequence, implementing the + * interpolation requires finding the MPI process that owns the + * cell on which that point is located, and then communication + * back and forth. + * + * @sect3{Helper functions} + * + * This class implements the interpolation functionality for + * derived classes to use. It works through a two-stage approach: + * First, derived classes declare at which points they require + * ASPECT's solution to be evaluated, via the + * set_evaluation_points() function. This function then finds + * which process owns the cell around this point, and sets up + * communication structures that will make later evaluation + * efficient. Second, this class provides the + * evaluate_aspect_variables_at_points() function that uses these + * communication structures to evaluate ASPECT's current solution + * at the points previously set. The class also provides the + * interpolate_vertical_velocities_to_surface_points() function + * that takes a set of (vertical) velocity values at these points + * and uses them to interpolate the information back onto ASPECT's + * mesh. + * + * All three of these functions are `protected` member functions + * of this class, ready to be called by derived classes. + * + * + * @sect3{A high-level function} + * + * All classes derived from Interface need to implement the + * Interface::compute_velocity_constraints_on_boundary() + * function. Because the basic outline of this function looks + * essentially the same for all external tools, this class also + * provides a high-level implementation of this function is also + * provided as part of this class. In essence, it looks as + * follows (pseudo-code), relying on the + * `compute_updated_velocities_at_points()` function that gets + * ASPECT's solution at the evaluation points and returns velocity + * vectors at all of these evaluation points: + * @code + * // Interpolate ASPECT's solution at evaluation points: + * aspect_surface_velocities = evaluate_aspect_variables_at_points(); + * + * // Call derive class's method to compute updated velocities: + * external_surface_velocities = compute_updated_velocities_at_points(aspect_surface_velocities); + * + * // Turn the result into the constraints that + * // compute_velocity_constraints_on_boundary is supposed + * // to return. + * @endcode + */ + template + class ExternalToolInterface : public Interface, public SimulatorAccess + { + public: + virtual + void + compute_velocity_constraints_on_boundary(const DoFHandler &mesh_deformation_dof_handler, + AffineConstraints &mesh_velocity_constraints, + const std::set &boundary_ids) const override; + + virtual + std::vector> + compute_updated_velocities_at_points (const std::vector> ¤t_solution_at_points) const = 0; + + + protected: + // Helper functions to be used in derived classes: + + /** + * Declare at which points the external tool driven by a class + * derived from the current one needs to evaluate the ASPECT + * solution. Each process will call this function with the points + * it needs. Different processes will, in general, provide + * distinct sets of points. + * + * The points provided here are given in the undeformed + * configuration that corresponds to the initial mesh. For + * example, if the mesh describes a box geometry, then the + * points should like in the $x$-$y$-plane that forms the top + * surface, rather than on the currently deformed top surface + * that describes the current elevation map. + * + * @note This function sets up communication structures that + * encode, among other things, which process owns which + * points. This information becomes outdated when the + * ASPECT mesh changes for mesh refinement. As a + * consequence, the current function sets up a process that + * invalidates the communication structures upon mesh + * refinement. Derived classes therefore have to set up + * their own code to call set_evaluation_points() again when + * the ASPECT mesh changes, or perhaps simply check whether + * the information needs to be updated in an overload of the + * update() function called at the beginning of each time + * step. + */ + void + set_evaluation_points (const std::vector> &evaluation_points); + + /** + * Return the value of the ASPECT solution at the set of points + * previously set via set_evaluation_points(). + * + * For each point, the return object contains a vector with as + * many components as there are ASPECT solution components. + */ + std::vector> + evaluate_aspect_variables_at_points () const; + + /** + * Given velocities (typically vertical, but not always) at + * all of the points this process has previously set via + * set_evaluation_points(), compute a (global) finite element + * field vector that in the velocity components of surface + * nodes corresponds to an interpolation of the velocities + * provided. + * + * The output of this function can then be used as input for a + * function that implements the + * compute_velocity_constraints_on_boundary() function of the + * base class. + */ + LinearAlgebra::Vector + interpolate_velocities_to_surface_points (const std::vector> &vertical_velocities) const; + + + private: + + std::vector> evaluation_points; + + // Timo: Replace by whatever type you need here + class X {}; + std::unique_ptr remote_point_evaluator; + }; + } +} + +#endif diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc new file mode 100644 index 00000000000..1f3931183f8 --- /dev/null +++ b/source/mesh_deformation/external_tool_interface.cc @@ -0,0 +1,159 @@ +/* + Copyright (C) 2014 - 2024 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . + */ + + +#include + +namespace aspect +{ + namespace MeshDeformation + { + template + void + ExternalToolInterface:: + compute_velocity_constraints_on_boundary(const DoFHandler &mesh_deformation_dof_handler, + AffineConstraints &mesh_velocity_constraints, + const std::set &boundary_ids) const + { + // First compute a (global) vector that has the correct velocities + // set at all boundary nodes: + const std::vector> aspect_surface_velocities = evaluate_aspect_variables_at_points(); + + const std::vector> external_surface_velocities + = compute_updated_velocities_at_points(aspect_surface_velocities); + + const LinearAlgebra::Vector v_interpolated + = interpolate_velocities_to_surface_points(external_surface_velocities); + + // Turn v_interpolated into constraints. For this, loop over all + // boundary DoFs and if a boundary DoF is locally owned, create a + // constraint. We later make that consistent across processors to + // ensure we also know about the locally relevant DoFs' + // constraints: + std::vector face_dof_indices (mesh_deformation_dof_handler.get_fe().dofs_per_face); + for (const auto &cell : mesh_deformation_dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + for (const auto &face : cell->face_iterators()) + if (face->at_boundary() && + (boundary_ids.find(face->boundary_id()) != boundary_ids.end())) + { + face->get_dof_indices(face_dof_indices); + + for (types::global_dof_index i : face_dof_indices) + if (v_interpolated.locally_owned_elements().is_element(i)) + mesh_velocity_constraints.add_constraint (i, {}, v_interpolated(i)); + } + mesh_velocity_constraints.make_consistent_in_parallel(mesh_deformation_dof_handler.locally_owned_dofs(), + DoFTools::extract_locally_relevant_dofs(mesh_deformation_dof_handler), + this->get_mpi_communicator()); + + } + + + template + void + ExternalToolInterface:: + set_evaluation_points (const std::vector> &evaluation_points) + { + // First, save a copy of the points at which we need the solution, + // among other reasons so that we can track that input arguments + // for later function calls describe the same number of points. + this->evaluation_points = evaluation_points; + + // Then also invalidate the previous evaluator: + remote_point_evaluator.reset(); + + // Timo: + // TODO: Implement set-up phase via MPIRemotePointEvaluation + + + // Finally, also ensure that upon mesh refinement, all of the + // information set herein is invalidated: + this->get_signals().pre_refinement_store_user_data + .connect([this](typename parallel::distributed::Triangulation &) + { + this->evaluation_points.clear(); + this->remote_point_evaluator.reset(); + } + ); + } + + + + template + std::vector> + ExternalToolInterface:: + evaluate_aspect_variables_at_points () const + { + Assert (remote_point_evaluator != nullptr, + ExcMessage("You can only call this function if you have previously " + "set the evaluation points by calling set_evaluation_points(), " + "and if the evaluator has not been invalidated by a mesh " + "refinement step.")); + + std::vector> solution_at_points (evaluation_points.size(), + std::vector(this->introspection().n_components)); + // Timo: Implement via MPIRemoteEvaluation + + return solution_at_points; + } + + + + template + LinearAlgebra::Vector + ExternalToolInterface:: + interpolate_velocities_to_surface_points (const std::vector> &velocities) const + { + Assert (remote_point_evaluator != nullptr, + ExcMessage("You can only call this function if you have previously " + "set the evaluation points by calling set_evaluation_points(), " + "and if the evaluator has not been invalidated by a mesh " + "refinement step.")); + AssertDimension(velocities.size(), evaluation_points.size()); + + // Create the output vector. TODO: We need to get access to the + // locally owned DoFs index set. Look up how this is done in the + // other implementations of the Interface base class, if they do + // it this way at all. + + // LinearAlgebra::Vector vector_with_surface_velocities(this->get_mesh_deformation_handler().mesh_deformation_dof_handler.locally_owned_dofs(), + // mpi_communicator) + LinearAlgebra::Vector vector_with_surface_velocities; + + + // Timo: Implement via MPIRemoteEvaluation + + return vector_with_surface_velocities; + } + } + + + + namespace MeshDeformation + { +#define INSTANTIATE(dim) \ + template class ExternalToolInterface; + + ASPECT_INSTANTIATE(INSTANTIATE) + +#undef INSTANTIATE + } +} From b8ac48f452c01e309c9fd4e42d0684a496a0e341 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 30 Jun 2025 08:17:23 -0600 Subject: [PATCH 2/6] Use the same style to compute constraints as in all other mesh deformation plugins. --- .../external_tool_interface.cc | 35 +++++++++++-------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc index 1f3931183f8..27292688781 100644 --- a/source/mesh_deformation/external_tool_interface.cc +++ b/source/mesh_deformation/external_tool_interface.cc @@ -47,23 +47,28 @@ namespace aspect // constraint. We later make that consistent across processors to // ensure we also know about the locally relevant DoFs' // constraints: - std::vector face_dof_indices (mesh_deformation_dof_handler.get_fe().dofs_per_face); - for (const auto &cell : mesh_deformation_dof_handler.active_cell_iterators()) - if (cell->is_locally_owned()) - for (const auto &face : cell->face_iterators()) - if (face->at_boundary() && - (boundary_ids.find(face->boundary_id()) != boundary_ids.end())) + // now insert the relevant part of the solution into the mesh constraints + const IndexSet constrained_dofs = + DoFTools::extract_boundary_dofs(mesh_deformation_dof_handler, + ComponentMask(dim, true), + boundary_ids); + + for (unsigned int i = 0; i < constrained_dofs.n_elements(); ++i) + { + types::global_dof_index index = constrained_dofs.nth_index_in_set(i); + if (mesh_velocity_constraints.can_store_line(index)) + if (mesh_velocity_constraints.is_constrained(index)==false) { - face->get_dof_indices(face_dof_indices); - - for (types::global_dof_index i : face_dof_indices) - if (v_interpolated.locally_owned_elements().is_element(i)) - mesh_velocity_constraints.add_constraint (i, {}, v_interpolated(i)); +#if DEAL_II_VERSION_GTE(9,6,0) + mesh_velocity_constraints.add_constraint(index, + {}, + v_interpolated(index)); +#else + mesh_velocity_constraints.add_line(index); + mesh_velocity_constraints.set_inhomogeneity(index, v_interpolated(index)); +#endif } - mesh_velocity_constraints.make_consistent_in_parallel(mesh_deformation_dof_handler.locally_owned_dofs(), - DoFTools::extract_locally_relevant_dofs(mesh_deformation_dof_handler), - this->get_mpi_communicator()); - + } } From c6a72d1f5b1ff29afca257ee1737b1a533eea08f Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 30 Jun 2025 08:23:15 -0600 Subject: [PATCH 3/6] Make index handling more efficient. --- source/mesh_deformation/external_tool_interface.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc index 27292688781..74225a52fc3 100644 --- a/source/mesh_deformation/external_tool_interface.cc +++ b/source/mesh_deformation/external_tool_interface.cc @@ -53,9 +53,8 @@ namespace aspect ComponentMask(dim, true), boundary_ids); - for (unsigned int i = 0; i < constrained_dofs.n_elements(); ++i) + for (const types::global_dof_index index : constrained_dofs) { - types::global_dof_index index = constrained_dofs.nth_index_in_set(i); if (mesh_velocity_constraints.can_store_line(index)) if (mesh_velocity_constraints.is_constrained(index)==false) { From 1925254d1efa7fe403cdcd9faa03d6b7d4705525 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 30 Jun 2025 08:40:33 -0600 Subject: [PATCH 4/6] Silence a compiler warning. --- source/mesh_deformation/external_tool_interface.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc index 74225a52fc3..e31ab5c2dfb 100644 --- a/source/mesh_deformation/external_tool_interface.cc +++ b/source/mesh_deformation/external_tool_interface.cc @@ -144,6 +144,8 @@ namespace aspect // Timo: Implement via MPIRemoteEvaluation + (void)velocities; + return vector_with_surface_velocities; } From b2ae4a38534742608e58cbf8bb02f5f11307cd8f Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 30 Jun 2025 11:06:19 -0600 Subject: [PATCH 5/6] Fix whitespace issue. --- source/mesh_deformation/external_tool_interface.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc index e31ab5c2dfb..3502366a723 100644 --- a/source/mesh_deformation/external_tool_interface.cc +++ b/source/mesh_deformation/external_tool_interface.cc @@ -145,7 +145,7 @@ namespace aspect // Timo: Implement via MPIRemoteEvaluation (void)velocities; - + return vector_with_surface_velocities; } From 776a36df54af91671894fd2936b1040d12a2932e Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 30 Jun 2025 11:07:08 -0600 Subject: [PATCH 6/6] Add necessary header include. --- source/mesh_deformation/external_tool_interface.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc index 3502366a723..0e72f3e2ab8 100644 --- a/source/mesh_deformation/external_tool_interface.cc +++ b/source/mesh_deformation/external_tool_interface.cc @@ -20,6 +20,8 @@ #include +#include + namespace aspect {