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..0e72f3e2ab8 --- /dev/null +++ b/source/mesh_deformation/external_tool_interface.cc @@ -0,0 +1,167 @@ +/* + 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 +#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: + // 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 (const types::global_dof_index index : constrained_dofs) + { + if (mesh_velocity_constraints.can_store_line(index)) + if (mesh_velocity_constraints.is_constrained(index)==false) + { +#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 + } + } + } + + + 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 + (void)velocities; + + + return vector_with_surface_velocities; + } + } + + + + namespace MeshDeformation + { +#define INSTANTIATE(dim) \ + template class ExternalToolInterface; + + ASPECT_INSTANTIATE(INSTANTIATE) + +#undef INSTANTIATE + } +}