Merge pull request #5575 from nilsvu/elliptic_time_dependent_maps

Support domain deformations in elliptic solver

cmake/SetupCxxFlags.cmake

@@ -103,6 +103,16 @@ set_property(TARGET SpectreFlags
   INTERFACE_COMPILE_OPTIONS
   $<$<COMPILE_LANGUAGE:CXX>:-ftemplate-backtrace-limit=0>)
 
+# Increase bracket depth for fold expressions
+# (see also https://github.com/llvm/llvm-project/issues/48973)
+if (CMAKE_CXX_COMPILER_ID MATCHES "Clang"
+    AND CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 12.0)
+  set_property(TARGET SpectreFlags
+    APPEND PROPERTY
+    INTERFACE_COMPILE_OPTIONS
+    $<$<COMPILE_LANGUAGE:CXX>:-fbracket-depth=1024>)
+endif()
+
 # Disable cmath setting the error flag. This allows the compiler to more
 # aggressively vectorize code since it doesn't need to respect some global
 # state.

src/Domain/CMakeLists.txt

@@ -70,15 +70,14 @@ target_link_libraries(
   DomainBoundaryConditions
   DomainStructure
   ErrorHandling
+  FunctionsOfTime
   Options
   Spectral
   Utilities
   PRIVATE
   LinearOperators
   RootFinding
   SphericalHarmonics
-  INTERFACE
-  FunctionsOfTime
 )
 
 add_subdirectory(Amr)

src/Domain/CoordinateMaps/Composition.cpp

@@ -365,6 +365,15 @@ bool Composition<Frames, Dim, std::index_sequence<Is...>>::is_equal_to(
 #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
 
 #define INSTANTIATE(_, data)                                                   \
+  template class Composition<                                                  \
+      tmpl::list<Frame::BlockLogical, Frame::Grid, Frame::Inertial>,           \
+      DIM(data)>;                                                              \
+  template class Composition<                                                  \
+      tmpl::list<Frame::BlockLogical, Frame::Grid, Frame::Distorted>,          \
+      DIM(data)>;                                                              \
+  template class Composition<tmpl::list<Frame::BlockLogical, Frame::Grid,      \
+                                        Frame::Distorted, Frame::Inertial>,    \
+                             DIM(data)>;                                       \
   template class Composition<                                                  \
       tmpl::list<Frame::ElementLogical, Frame::BlockLogical, Frame::Inertial>, \
       DIM(data)>;                                                              \
@@ -390,6 +399,48 @@ GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))
 
 #if defined(__clang__) && __clang_major__ >= 15
 #define INSTANTIATE2(_, data)                                                  \
+  template domain::CoordinateMaps::Composition<                                \
+      brigand::list<Frame::BlockLogical, Frame::Grid, Frame::Distorted>,       \
+      DIM(data), std::integer_sequence<unsigned long, 0ul, 1ul>>::             \
+      Composition(                                                             \
+          std::unique_ptr<domain::CoordinateMapBase<Frame::BlockLogical,       \
+                                                    Frame::Grid, DIM(data)>,   \
+                          std::default_delete<domain::CoordinateMapBase<       \
+                              Frame::BlockLogical, Frame::Grid, DIM(data)>>>,  \
+          std::unique_ptr<domain::CoordinateMapBase<                           \
+                              Frame::Grid, Frame::Distorted, DIM(data)>,       \
+                          std::default_delete<domain::CoordinateMapBase<       \
+                              Frame::Grid, Frame::Distorted, DIM(data)>>>);    \
+  template domain::CoordinateMaps::Composition<                                \
+      brigand::list<Frame::BlockLogical, Frame::Grid, Frame::Inertial>,        \
+      DIM(data), std::integer_sequence<unsigned long, 0ul, 1ul>>::             \
+      Composition(                                                             \
+          std::unique_ptr<domain::CoordinateMapBase<Frame::BlockLogical,       \
+                                                    Frame::Grid, DIM(data)>,   \
+                          std::default_delete<domain::CoordinateMapBase<       \
+                              Frame::BlockLogical, Frame::Grid, DIM(data)>>>,  \
+          std::unique_ptr<domain::CoordinateMapBase<                           \
+                              Frame::Grid, Frame::Inertial, DIM(data)>,        \
+                          std::default_delete<domain::CoordinateMapBase<       \
+                              Frame::Grid, Frame::Inertial, DIM(data)>>>);     \
+  template domain::CoordinateMaps::Composition<                                \
+      brigand::list<Frame::BlockLogical, Frame::Grid, Frame::Distorted,        \
+                    Frame::Inertial>,                                          \
+      DIM(data), std::integer_sequence<unsigned long, 0ul, 1ul, 2ul>>::        \
+      Composition(                                                             \
+          std::unique_ptr<domain::CoordinateMapBase<Frame::BlockLogical,       \
+                                                    Frame::Grid, DIM(data)>,   \
+                          std::default_delete<domain::CoordinateMapBase<       \
+                              Frame::BlockLogical, Frame::Grid, DIM(data)>>>,  \
+          std::unique_ptr<domain::CoordinateMapBase<                           \
+                              Frame::Grid, Frame::Distorted, DIM(data)>,       \
+                          std::default_delete<domain::CoordinateMapBase<       \
+                              Frame::Grid, Frame::Distorted, DIM(data)>>>,     \
+          std::unique_ptr<                                                     \
+              domain::CoordinateMapBase<Frame::Distorted, Frame::Inertial,     \
+                                        DIM(data)>,                            \
+              std::default_delete<domain::CoordinateMapBase<                   \
+                  Frame::Distorted, Frame::Inertial, DIM(data)>>>);            \
   template domain::CoordinateMaps::Composition<                                \
       brigand::list<Frame::ElementLogical, Frame::BlockLogical,                \
                     Frame::Inertial>,                                          \

src/Domain/ElementDistribution.cpp

@@ -63,10 +63,7 @@ double get_num_points_and_grid_spacing_cost(
       initial_extents, element_id, quadrature);
   Element<Dim> element = ::domain::Initialization::create_initial_element(
       element_id, block, initial_refinement_levels);
-  ElementMap<Dim, Frame::Grid> element_map{
-      element_id, block.is_time_dependent()
-                      ? block.moving_mesh_logical_to_grid_map().get_clone()
-                      : block.stationary_map().get_to_grid_frame()};
+  ElementMap<Dim, Frame::Grid> element_map{element_id, block};
   const tnsr::I<DataVector, Dim, Frame::ElementLogical> logical_coords =
       logical_coordinates(mesh);
   const tnsr::I<DataVector, Dim, Frame::Grid> grid_coords =

src/Domain/ElementMap.cpp

@@ -3,8 +3,10 @@
 
 #include "Domain/ElementMap.hpp"
 
+#include "Domain/CoordinateMaps/Composition.hpp"
 #include "Domain/CoordinateMaps/CoordinateMap.hpp"  // IWYU pragma: keep
 #include "Domain/Structure/Side.hpp"
+#include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/OptimizerHacks.hpp"
 
@@ -53,6 +55,40 @@ ElementMap<Dim, TargetFrame>::ElementMap(
       jacobian_{map_slope_},
       inverse_jacobian_{map_inverse_slope_} {}
 
+// We could refactor the ElementMap class to use a `Composition` internally also
+// for the element-to-block-logical map, so the whole map is just one
+// composition.
+template <size_t Dim, typename TargetFrame>
+ElementMap<Dim, TargetFrame>::ElementMap(ElementId<Dim> element_id,
+                                         const Block<Dim>& block)
+    : ElementMap(
+          std::move(element_id),
+          [&block]() -> std::unique_ptr<domain::CoordinateMapBase<
+                         Frame::BlockLogical, TargetFrame, Dim>> {
+            if constexpr (std::is_same_v<TargetFrame, Frame::Inertial>) {
+              if (block.is_time_dependent()) {
+                using CompositionType = domain::CoordinateMaps::Composition<
+                    tmpl::list<Frame::BlockLogical, Frame::Grid,
+                               Frame::Inertial>,
+                    Dim>;
+                return std::make_unique<CompositionType>(
+                    block.moving_mesh_logical_to_grid_map().get_clone(),
+                    block.moving_mesh_grid_to_inertial_map().get_clone());
+              } else {
+                return block.stationary_map().get_clone();
+              }
+            } else if constexpr (std::is_same_v<TargetFrame, Frame::Grid>) {
+              if (block.is_time_dependent()) {
+                return block.moving_mesh_logical_to_grid_map().get_clone();
+              } else {
+                return block.stationary_map().get_to_grid_frame();
+              }
+            }
+          }()) {
+  ASSERT(element_id_.block_id() == block.id(),
+         "Element " << element_id_ << " is not in block " << block.id() << ".");
+}
+
 template <size_t Dim, typename TargetFrame>
 void ElementMap<Dim, TargetFrame>::pup(PUP::er& p) {
   p | block_map_;

src/Domain/ElementMap.hpp

@@ -9,7 +9,9 @@
 #include <memory>
 
 #include "DataStructures/Tensor/Tensor.hpp"
+#include "Domain/Block.hpp"
 #include "Domain/CoordinateMaps/CoordinateMap.hpp"
+#include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
 #include "Domain/Structure/ElementId.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/MakeArray.hpp"
@@ -45,6 +47,14 @@ class ElementMap {
                                                        TargetFrame, Dim>>
                  block_map);
 
+  /// Construct from an `element_id` within the `block`. The (affine)
+  /// ElementLogical to BlockLogical map is determined by the `element_id`. The
+  /// BlockLogical to TargetFrame map is determined by the `block`:
+  /// - If the block is time-independent: the `block.stationary_map()` is used.
+  /// - If the block is time-dependent: The `block.moving_mesh_*_map()` maps
+  ///   are used. Which maps are used depends on the TargetFrame.
+  ElementMap(ElementId<Dim> element_id, const Block<Dim>& block);
+
   const domain::CoordinateMapBase<Frame::BlockLogical, TargetFrame, Dim>&
   block_map() const {
     return *block_map_;
@@ -54,17 +64,23 @@ class ElementMap {
 
   template <typename T>
   tnsr::I<T, Dim, TargetFrame> operator()(
-      const tnsr::I<T, Dim, Frame::ElementLogical>& source_point) const {
+      const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
+      const double time = std::numeric_limits<double>::signaling_NaN(),
+      const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
     auto block_source_point =
         apply_affine_transformation_to_point(source_point);
-    return block_map_->operator()(std::move(block_source_point));
+    return block_map_->operator()(std::move(block_source_point), time,
+                                  functions_of_time);
   }
 
   template <typename T>
   tnsr::I<T, Dim, Frame::ElementLogical> inverse(
-      tnsr::I<T, Dim, TargetFrame> target_point) const {
+      tnsr::I<T, Dim, TargetFrame> target_point,
+      const double time = std::numeric_limits<double>::signaling_NaN(),
+      const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
     auto block_source_point{
-        block_map_->inverse(std::move(target_point)).value()};
+        block_map_->inverse(std::move(target_point), time, functions_of_time)
+            .value()};
     // Apply the affine map to the points
     tnsr::I<T, Dim, Frame::ElementLogical> source_point;
     for (size_t d = 0; d < Dim; ++d) {
@@ -77,11 +93,13 @@ class ElementMap {
 
   template <typename T>
   InverseJacobian<T, Dim, Frame::ElementLogical, TargetFrame> inv_jacobian(
-      const tnsr::I<T, Dim, Frame::ElementLogical>& source_point) const {
+      const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
+      const double time = std::numeric_limits<double>::signaling_NaN(),
+      const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
     auto block_source_point =
         apply_affine_transformation_to_point(source_point);
-    auto block_inv_jac =
-        block_map_->inv_jacobian(std::move(block_source_point));
+    auto block_inv_jac = block_map_->inv_jacobian(std::move(block_source_point),
+                                                  time, functions_of_time);
     InverseJacobian<T, Dim, Frame::ElementLogical, TargetFrame> inv_jac;
     for (size_t d = 0; d < Dim; ++d) {
       for (size_t i = 0; i < Dim; ++i) {
@@ -94,10 +112,13 @@ class ElementMap {
 
   template <typename T>
   Jacobian<T, Dim, Frame::ElementLogical, TargetFrame> jacobian(
-      const tnsr::I<T, Dim, Frame::ElementLogical>& source_point) const {
+      const tnsr::I<T, Dim, Frame::ElementLogical>& source_point,
+      const double time = std::numeric_limits<double>::signaling_NaN(),
+      const domain::FunctionsOfTimeMap& functions_of_time = {}) const {
     auto block_source_point =
         apply_affine_transformation_to_point(source_point);
-    auto block_jac = block_map_->jacobian(std::move(block_source_point));
+    auto block_jac = block_map_->jacobian(std::move(block_source_point), time,
+                                          functions_of_time);
     Jacobian<T, Dim, Frame::ElementLogical, TargetFrame> jac;
     for (size_t d = 0; d < Dim; ++d) {
       for (size_t i = 0; i < Dim; ++i) {

src/Domain/FaceNormal.cpp

@@ -198,32 +198,40 @@ GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3),
 
 #undef INSTANTIATION
 
-#define INSTANTIATION(_, data)                                              \
-  template void unnormalized_face_normal(                                   \
-      const gsl::not_null<                                                  \
-          tnsr::i<DataVector, GET_DIM(data), Frame::Inertial>*>             \
-          result,                                                           \
-      const Mesh<GET_DIM(data) - 1>& interface_mesh,                        \
-      const ElementMap<GET_DIM(data), Frame::Grid>& logical_to_grid_map,    \
-      const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial,         \
-                                      GET_DIM(data)>& grid_to_inertial_map, \
-      const double time,                                                    \
-      const std::unordered_map<                                             \
-          std::string,                                                      \
-          std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&        \
-          functions_of_time,                                                \
-      const Direction<GET_DIM(data)>& direction);                           \
-  template tnsr::i<DataVector, GET_DIM(data), Frame::Inertial>              \
-  unnormalized_face_normal(                                                 \
-      const Mesh<GET_DIM(data) - 1>& interface_mesh,                        \
-      const ElementMap<GET_DIM(data), Frame::Grid>& logical_to_grid_map,    \
-      const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial,         \
-                                      GET_DIM(data)>& grid_to_inertial_map, \
-      const double time,                                                    \
-      const std::unordered_map<                                             \
-          std::string,                                                      \
-          std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&        \
-          functions_of_time,                                                \
+#define INSTANTIATION(_, data)                                                \
+  template void unnormalized_face_normal(                                     \
+      const gsl::not_null<                                                    \
+          tnsr::i<DataVector, GET_DIM(data), Frame::Inertial>*>               \
+          result,                                                             \
+      const Mesh<GET_DIM(data) - 1>& interface_mesh,                          \
+      const ElementMap<GET_DIM(data), Frame::Grid>& logical_to_grid_map,      \
+      const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial,           \
+                                      GET_DIM(data)>& grid_to_inertial_map,   \
+      const double time,                                                      \
+      const std::unordered_map<                                               \
+          std::string,                                                        \
+          std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&          \
+          functions_of_time,                                                  \
+      const Direction<GET_DIM(data)>& direction);                             \
+  template void unnormalized_face_normal(                                     \
+      const gsl::not_null<                                                    \
+          tnsr::i<DataVector, GET_DIM(data), Frame::Inertial>*>               \
+          result,                                                             \
+      const Mesh<GET_DIM(data) - 1>& interface_mesh,                          \
+      const InverseJacobian<DataVector, GET_DIM(data), Frame::ElementLogical, \
+                            Frame::Inertial>& inv_jacobian_on_interface,      \
+      const Direction<GET_DIM(data)>& direction);                             \
+  template tnsr::i<DataVector, GET_DIM(data), Frame::Inertial>                \
+  unnormalized_face_normal(                                                   \
+      const Mesh<GET_DIM(data) - 1>& interface_mesh,                          \
+      const ElementMap<GET_DIM(data), Frame::Grid>& logical_to_grid_map,      \
+      const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial,           \
+                                      GET_DIM(data)>& grid_to_inertial_map,   \
+      const double time,                                                      \
+      const std::unordered_map<                                               \
+          std::string,                                                        \
+          std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&          \
+          functions_of_time,                                                  \
       const Direction<GET_DIM(data)>& direction);
 
 GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3))

src/Domain/FunctionsOfTime/FunctionOfTime.hpp

@@ -84,4 +84,8 @@ class FunctionOfTime : public PUP::able {
   WRAPPED_PUPable_abstract(FunctionOfTime);  // NOLINT
 };
 }  // namespace FunctionsOfTime
+
+using FunctionsOfTimeMap = std::unordered_map<
+    std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>;
+
 }  // namespace domain

src/Domain/FunctionsOfTime/FunctionsOfTimeAreReady.hpp

@@ -46,7 +46,7 @@ bool functions_of_time_are_ready_impl(
     const Component* /*meta*/, const double time,
     const std::optional<std::unordered_set<std::string>>& functions_to_check,
     Args&&... args) {
-  if constexpr (Parallel::is_in_global_cache<Metavariables, CacheTag>) {
+  if constexpr (Parallel::is_in_mutable_global_cache<Metavariables, CacheTag>) {
     const auto& proxy =
         ::Parallel::get_parallel_component<Component>(cache)[array_index];
     const Parallel::ArrayComponentId array_component_id =