Re-initialize Schwarz subdomain before each solve

To account for changes in domain geometry (by AMR).

src/Elliptic/DiscontinuousGalerkin/Initialization.hpp

@@ -103,6 +103,7 @@ struct InitializeGeometry {
                  domain::Tags::InitialRefinementLevels<Dim>,
                  domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime,
                  elliptic::dg::Tags::Quadrature>;
+  using volume_tags = argument_tags;
   static void apply(
       gsl::not_null<Mesh<Dim>*> mesh, gsl::not_null<Element<Dim>*> element,
       gsl::not_null<DirectionalIdMap<Dim, Mesh<Dim>>*> neighbor_meshes,
@@ -140,6 +141,8 @@ struct ProjectGeometry : tt::ConformsTo<::amr::protocols::Projector> {
   using argument_tags =
       tmpl::list<domain::Tags::Mesh<Dim>, domain::Tags::Element<Dim>,
                  domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime>;
+  using volume_tags =
+      tmpl::list<domain::Tags::Domain<Dim>, domain::Tags::FunctionsOfTime>;
 
   // p-refinement
   static void apply(

src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/InitializeSubdomain.hpp

@@ -107,19 +107,36 @@ struct InitializeOverlapGeometry {
  * Note that the geometry depends on the system and on the choice of background
  * through the normalization of face normals, which involves a metric.
  *
+ * The `FromInitialDomain` template parameter controls how the
+ * next-to-nearest-neighbor information is acquired:
+ * - If `FromInitialDomain = true`, uses the initial parameters from the input
+ *   file. This can be used at the beginning of a simulation or if the domain
+ *   never changes (no AMR).
+ * - If `FromInitialDomain = false`, uses the mesh, element, and neighbor meshes
+ *   of each neighbor. This typically requires communication to exchange this
+ *   data between neighbors each time the domain changes, e.g. with AMR.
+ *
  * DataBox:
  * - Uses:
  *   - `BackgroundTag`
  *   - `domain::Tags::Element<Dim>`
- *   - `domain::Tags::InitialExtents<Dim>`
- *   - `domain::Tags::InitialRefinementLevels<Dim>`
+ *   - If `FromInitialDomain = true`:
+ *     - `domain::Tags::InitialExtents<Dim>`
+ *     - `domain::Tags::InitialRefinementLevels<Dim>`
+ *     - `elliptic::dg::Tags::Quadrature`
+ *   - If `FromInitialDomain = false`:
+ *     - `Overlaps<domain::Tags::Mesh<Dim>>`
+ *     - `Overlaps<domain::Tags::Element<Dim>>`
+ *     - `Overlaps<domain::Tags::NeighborMesh<Dim>>`
  *   - `domain::Tags::Domain<Dim>`
+ *   - `domain::Tags::FunctionsOfTime`
  *   - `LinearSolver::Schwarz::Tags::MaxOverlap<OptionsGroup>`
  * - Adds: Many tags prefixed with `LinearSolver::Schwarz::Tags::Overlaps`. See
  *   `elliptic::dg::Actions::InitializeDomain` and
  *   `elliptic::dg::Actions::initialize_operator` for a complete list.
  */
-template <typename System, typename BackgroundTag, typename OptionsGroup>
+template <typename System, typename BackgroundTag, typename OptionsGroup,
+          bool FromInitialDomain = true>
 struct InitializeSubdomain {
  private:
   static constexpr size_t Dim = System::volume_dim;
@@ -128,7 +145,10 @@ struct InitializeSubdomain {
   static constexpr bool has_background_fields =
       not std::is_same_v<typename System::background_fields, tmpl::list<>>;
 
-  using InitializeGeometry = elliptic::dg::InitializeGeometry<Dim>;
+  using InitializeGeometry =
+      tmpl::conditional_t<FromInitialDomain,
+                          elliptic::dg::InitializeGeometry<Dim>,
+                          elliptic::dg::ProjectGeometry<Dim>>;
   using InitializeOverlapGeometry = detail::InitializeOverlapGeometry<Dim>;
   using InitializeFacesAndMortars = elliptic::dg::InitializeFacesAndMortars<
       Dim, typename System::inv_metric_tag, BackgroundTag>;
@@ -186,8 +206,13 @@ struct InitializeSubdomain {
         elliptic::util::mutate_apply_at<
             db::wrap_tags_in<overlaps_tag,
                              typename InitializeGeometry::return_tags>,
-            typename InitializeGeometry::argument_tags,
-            typename InitializeGeometry::argument_tags>(
+            tmpl::append<
+                db::wrap_tags_in<overlaps_tag,
+                                 tmpl::list_difference<
+                                     typename InitializeGeometry::argument_tags,
+                                     typename InitializeGeometry::volume_tags>>,
+                typename InitializeGeometry::volume_tags>,
+            typename InitializeGeometry::volume_tags>(
             InitializeGeometry{}, make_not_null(&box), overlap_id, neighbor_id);
         // Initialize subdomain-specific tags on overlaps
         elliptic::util::mutate_apply_at<

src/Elliptic/Executables/Solver.hpp

@@ -213,8 +213,6 @@ struct Solver {
                        typename system::primal_fluxes>,
           elliptic::analytic_data::AnalyticSolution>,
       elliptic::dg::Actions::initialize_operator<system, background_tag>,
-      elliptic::dg::subdomain_operator::Actions::InitializeSubdomain<
-          system, background_tag, typename schwarz_smoother::options_group>,
       tmpl::conditional_t<
           is_linear, tmpl::list<>,
           ::Initialization::Actions::AddComputeTags<tmpl::list<
@@ -251,6 +249,14 @@ struct Solver {
   using smooth_actions = typename schwarz_smoother::template solve<
       typename dg_operator<true>::apply_actions, Label>;
 
+  // These tags are communicated on subdomain overlaps to initialize the
+  // subdomain geometry. AMR updates these tags, so we have to communicate them
+  // after each AMR step.
+  using subdomain_init_tags =
+      tmpl::list<domain::Tags::Mesh<volume_dim>,
+                 domain::Tags::Element<volume_dim>,
+                 domain::Tags::NeighborMesh<volume_dim>>;
+
   /// This data needs to be communicated on subdomain overlap regions
   using communicated_overlap_tags = tmpl::flatten<tmpl::list<
       // For linearized sources
@@ -261,6 +267,11 @@ struct Solver {
           volume_dim, db::wrap_tags_in<::Tags::NormalDotFlux,
                                        typename system::primal_fields>>>>;
 
+  using init_subdomain_action =
+      elliptic::dg::subdomain_operator::Actions::InitializeSubdomain<
+          system, background_tag, typename schwarz_smoother::options_group,
+          false>;
+
   template <typename StepActions>
   using linear_solve_actions = typename linear_solver::template solve<
       tmpl::list<
@@ -289,10 +300,12 @@ struct Solver {
           // Communicate data on subdomain overlap regions
           LinearSolver::Schwarz::Actions::SendOverlapFields<
               communicated_overlap_tags,
-              typename schwarz_smoother::options_group, false>,
+              typename schwarz_smoother::options_group, false,
+              nonlinear_solver_iteration_id>,
           LinearSolver::Schwarz::Actions::ReceiveOverlapFields<
               volume_dim, communicated_overlap_tags,
-              typename schwarz_smoother::options_group>,
+              typename schwarz_smoother::options_group, false,
+              nonlinear_solver_iteration_id>,
           // Reset Schwarz subdomain solver
           LinearSolver::Schwarz::Actions::ResetSubdomainSolver<
               typename schwarz_smoother::options_group>,
@@ -301,22 +314,33 @@ struct Solver {
       StepActions>;
 
   template <typename StepActions>
-  using solve_actions = tmpl::conditional_t<
-      is_linear,
-      // Linear solve
-      tmpl::list<
-          // Apply the DG operator to the initial guess
-          typename elliptic::dg::Actions::DgOperator<
-              system, true, linear_solver_iteration_id, fields_tag,
-              fluxes_vars_tag, operator_applied_to_fields_tag, vars_tag,
-              fluxes_vars_tag>::apply_actions,
-          // Modify fixed sources with boundary conditions
-          elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource<
-              system, fixed_sources_tag>,
-          // Krylov solve
-          linear_solve_actions<StepActions>>,
-      // Nonlinear solve
-      nonlinear_solve_actions<StepActions>>;
+  using solve_actions = tmpl::flatten<tmpl::list<
+      // Communicate subdomain geometry and initialize subdomain to account for
+      // domain changes
+      LinearSolver::Schwarz::Actions::SendOverlapFields<
+          subdomain_init_tags, typename schwarz_smoother::options_group, false>,
+      LinearSolver::Schwarz::Actions::ReceiveOverlapFields<
+          volume_dim, subdomain_init_tags,
+          typename schwarz_smoother::options_group, false>,
+      init_subdomain_action,
+      // Linear or nonlinear solve
+      tmpl::conditional_t<
+          is_linear,
+          // Linear solve
+          tmpl::list<
+              // Apply the DG operator to the initial guess
+              typename elliptic::dg::Actions::DgOperator<
+                  system, true, linear_solver_iteration_id, fields_tag,
+                  fluxes_vars_tag, operator_applied_to_fields_tag, vars_tag,
+                  fluxes_vars_tag>::apply_actions,
+              // Modify fixed sources with boundary conditions
+              elliptic::dg::Actions::
+                  ImposeInhomogeneousBoundaryConditionsOnSource<
+                      system, fixed_sources_tag>,
+              // Krylov solve
+              linear_solve_actions<StepActions>>,
+          // Nonlinear solve
+          nonlinear_solve_actions<StepActions>>>>;
 
   using component_list = tmpl::flatten<
       tmpl::list<tmpl::conditional_t<is_linear, tmpl::list<>,

src/ParallelAlgorithms/LinearSolver/Schwarz/Actions/CommunicateOverlapFields.hpp

@@ -70,14 +70,15 @@ struct OverlapFieldsTag
  * This actions should be followed by
  * `LinearSolver::Schwarz::Actions::ReceiveOverlapFields` in the action list.
  */
-template <typename OverlapFields, typename OptionsGroup, bool RestrictToOverlap>
+template <typename OverlapFields, typename OptionsGroup, bool RestrictToOverlap,
+          typename TemporalIdTag = Convergence::Tags::IterationId<OptionsGroup>>
 struct SendOverlapFields;
 
 /// \cond
 template <typename... OverlapFields, typename OptionsGroup,
-          bool RestrictToOverlap>
+          bool RestrictToOverlap, typename TemporalIdTag>
 struct SendOverlapFields<tmpl::list<OverlapFields...>, OptionsGroup,
-                         RestrictToOverlap> {
+                         RestrictToOverlap, TemporalIdTag> {
   using const_global_cache_tags =
       tmpl::list<Tags::MaxOverlap<OptionsGroup>,
                  logging::Tags::Verbosity<OptionsGroup>>;
@@ -92,15 +93,16 @@ struct SendOverlapFields<tmpl::list<OverlapFields...>, OptionsGroup,
       const ParallelComponent* const /*meta*/) {
     const auto& element = get<domain::Tags::Element<Dim>>(box);
 
-    // Skip communicating if the overlap is empty
-    if (UNLIKELY(db::get<Tags::MaxOverlap<OptionsGroup>>(box) == 0 or
+    // Skip communicating if the overlap is empty and we RestrictToOverlap,
+    // because we would end up with empty tensor data.
+    if (UNLIKELY((RestrictToOverlap and
+                  db::get<Tags::MaxOverlap<OptionsGroup>>(box) == 0) or
                  element.number_of_neighbors() == 0)) {
       return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     }
 
     // Do some logging
-    const auto& iteration_id =
-        get<Convergence::Tags::IterationId<OptionsGroup>>(box);
+    const auto& iteration_id = get<TemporalIdTag>(box);
     if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
                  ::Verbosity::Debug)) {
       Parallel::printf("%s %s(%zu): Send overlap fields\n", element_id,
@@ -172,12 +174,16 @@ struct SendOverlapFields<tmpl::list<OverlapFields...>, OptionsGroup,
  * This actions should be preceded by
  * `LinearSolver::Schwarz::Actions::SendOverlapFields` in the action list.
  */
-template <size_t Dim, typename OverlapFields, typename OptionsGroup>
+template <size_t Dim, typename OverlapFields, typename OptionsGroup,
+          bool RestrictToOverlap,
+          typename TemporalIdTag = Convergence::Tags::IterationId<OptionsGroup>>
 struct ReceiveOverlapFields;
 
 /// \cond
-template <size_t Dim, typename... OverlapFields, typename OptionsGroup>
-struct ReceiveOverlapFields<Dim, tmpl::list<OverlapFields...>, OptionsGroup> {
+template <size_t Dim, typename... OverlapFields, typename OptionsGroup,
+          bool RestrictToOverlap, typename TemporalIdTag>
+struct ReceiveOverlapFields<Dim, tmpl::list<OverlapFields...>, OptionsGroup,
+                            RestrictToOverlap, TemporalIdTag> {
  private:
   using overlap_fields_tag =
       detail::OverlapFieldsTag<Dim, tmpl::list<OverlapFields...>, OptionsGroup>;
@@ -197,12 +203,12 @@ struct ReceiveOverlapFields<Dim, tmpl::list<OverlapFields...>, OptionsGroup> {
       const Parallel::GlobalCache<Metavariables>& /*cache*/,
       const ElementId<Dim>& element_id, const ActionList /*meta*/,
       const ParallelComponent* const /*meta*/) {
-    const auto& iteration_id =
-        get<Convergence::Tags::IterationId<OptionsGroup>>(box);
+    const auto& iteration_id = get<TemporalIdTag>(box);
     const auto& element = get<domain::Tags::Element<Dim>>(box);
 
-    // Nothing to receive if overlap is empty
-    if (UNLIKELY(db::get<Tags::MaxOverlap<OptionsGroup>>(box) == 0 or
+    // Nothing to receive if nothing was sent in the action above
+    if (UNLIKELY((RestrictToOverlap and
+                  db::get<Tags::MaxOverlap<OptionsGroup>>(box) == 0) or
                  element.number_of_neighbors() == 0)) {
       return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     }
@@ -227,6 +233,7 @@ struct ReceiveOverlapFields<Dim, tmpl::list<OverlapFields...>, OptionsGroup> {
     db::mutate<Tags::Overlaps<OverlapFields, Dim, OptionsGroup>...>(
         [&received_overlap_fields](const auto... local_overlap_fields) {
           if constexpr (sizeof...(OverlapFields) > 1) {
+            (local_overlap_fields->clear(), ...);
             for (auto& [overlap_id, overlap_fields] : received_overlap_fields) {
               expand_pack((*local_overlap_fields)[overlap_id] =
                               std::move(get<OverlapFields>(overlap_fields))...);