Add SolvePunctures executable

src/Elliptic/Executables/CMakeLists.txt

@@ -3,4 +3,5 @@
 
 add_subdirectory(Elasticity)
 add_subdirectory(Poisson)
+add_subdirectory(Punctures)
 add_subdirectory(Xcts)

src/Elliptic/Executables/Punctures/CMakeLists.txt

@@ -0,0 +1,42 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+set(LIBS_TO_LINK
+  Convergence
+  CoordinateMaps
+  DataStructures
+  DiscontinuousGalerkin
+  Domain
+  DomainCreators
+  Elliptic
+  EllipticDg
+  EllipticDgSubdomainOperator
+  EllipticSubdomainPreconditioners
+  ErrorHandling
+  Events
+  EventsAndTriggers
+  Informer
+  Initialization
+  LinearOperators
+  MathFunctions
+  Observer
+  Options
+  Parallel
+  ParallelLinearSolver
+  ParallelMultigrid
+  ParallelNonlinearSolver
+  ParallelSchwarz
+  Utilities
+  Punctures
+  PuncturesAnalyticData
+  PuncturesBoundaryConditions
+  PuncturesSolutions
+  )
+
+add_spectre_parallel_executable(
+  SolvePunctures
+  SolvePunctures
+  Elliptic/Executables/Punctures
+  "Metavariables"
+  "${LIBS_TO_LINK}"
+  )

src/Elliptic/Executables/Punctures/SolvePunctures.hpp

@@ -0,0 +1,145 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <array>
+#include <cstddef>
+
+#include "Domain/Creators/Factory3D.hpp"
+#include "Domain/Creators/RegisterDerivedWithCharm.hpp"
+#include "Domain/Tags.hpp"
+#include "Elliptic/DiscontinuousGalerkin/DgElementArray.hpp"
+#include "Elliptic/Executables/NonlinearEllipticSolver.hpp"
+#include "Elliptic/SubdomainPreconditioners/RegisterDerived.hpp"
+#include "Elliptic/Systems/Punctures/BoundaryConditions/Flatness.hpp"
+#include "Elliptic/Systems/Punctures/FirstOrderSystem.hpp"
+#include "Elliptic/Triggers/Factory.hpp"
+#include "IO/Observer/Actions/RegisterEvents.hpp"
+#include "IO/Observer/Helpers.hpp"
+#include "IO/Observer/ObserverComponent.hpp"
+#include "Options/Options.hpp"
+#include "Options/Protocols/FactoryCreation.hpp"
+#include "Parallel/GlobalCache.hpp"
+#include "Parallel/InitializationFunctions.hpp"
+#include "Parallel/Phase.hpp"
+#include "Parallel/PhaseDependentActionList.hpp"
+#include "Parallel/RegisterDerivedClassesWithCharm.hpp"
+#include "ParallelAlgorithms/Actions/TerminatePhase.hpp"
+#include "ParallelAlgorithms/Events/Factory.hpp"
+#include "ParallelAlgorithms/Events/Tags.hpp"
+#include "ParallelAlgorithms/EventsAndTriggers/Actions/RunEventsAndTriggers.hpp"
+#include "ParallelAlgorithms/EventsAndTriggers/Completion.hpp"
+#include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
+#include "ParallelAlgorithms/EventsAndTriggers/Trigger.hpp"
+#include "ParallelAlgorithms/LinearSolver/Multigrid/ElementsAllocator.hpp"
+#include "ParallelAlgorithms/LinearSolver/Multigrid/Tags.hpp"
+#include "PointwiseFunctions/AnalyticData/Punctures/MultiplePunctures.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/Punctures/Flatness.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/Background.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/InitialGuess.hpp"
+#include "Utilities/Blas.hpp"
+#include "Utilities/ErrorHandling/FloatingPointExceptions.hpp"
+#include "Utilities/MemoryHelpers.hpp"
+#include "Utilities/ProtocolHelpers.hpp"
+#include "Utilities/TMPL.hpp"
+
+/// \cond
+struct Metavariables {
+  static constexpr Options::String help{"Solve for puncture initial data"};
+
+  static constexpr size_t volume_dim = 3;
+  using system = Punctures::FirstOrderSystem;
+  using solver = elliptic::nonlinear_solver::Solver<Metavariables>;
+
+  using observe_fields = tmpl::append<
+      typename system::primal_fields, typename system::background_fields,
+      // Add ADM quantities here
+      tmpl::list<domain::Tags::Coordinates<volume_dim, Frame::Inertial>>>;
+  using observer_compute_tags =
+      tmpl::list<::Events::Tags::ObserverMeshCompute<volume_dim>>;
+
+  struct factory_creation
+      : tt::ConformsTo<Options::protocols::FactoryCreation> {
+    using factory_classes = tmpl::map<
+        tmpl::pair<DomainCreator<volume_dim>, domain_creators<volume_dim>>,
+        tmpl::pair<elliptic::analytic_data::Background,
+                   tmpl::list<Punctures::AnalyticData::MultiplePunctures>>,
+        tmpl::pair<elliptic::analytic_data::InitialGuess,
+                   tmpl::list<Punctures::Solutions::Flatness>>,
+        tmpl::pair<elliptic::analytic_data::AnalyticSolution, tmpl::list<>>,
+        tmpl::pair<elliptic::BoundaryConditions::BoundaryCondition<volume_dim>,
+                   tmpl::list<Punctures::BoundaryConditions::Flatness>>,
+        tmpl::pair<
+            Event,
+            tmpl::flatten<tmpl::list<
+                Events::Completion,
+                dg::Events::field_observations<
+                    volume_dim, typename solver::nonlinear_solver_iteration_id,
+                    observe_fields, observer_compute_tags,
+                    LinearSolver::multigrid::Tags::IsFinestGrid>>>>,
+        tmpl::pair<Trigger,
+                   elliptic::Triggers::all_triggers<
+                       typename solver::nonlinear_solver::options_group>>>;
+  };
+
+  // Additional items to store in the global cache
+  using const_global_cache_tags = tmpl::list<>;
+
+  // Collect all reduction tags for observers
+  using observed_reduction_data_tags =
+      observers::collect_reduction_data_tags<tmpl::push_back<
+          tmpl::at<factory_creation::factory_classes, Event>, solver>>;
+
+  using initialization_actions =
+      tmpl::push_back<typename solver::initialization_actions,
+                      Parallel::Actions::TerminatePhase>;
+
+  using register_actions =
+      tmpl::push_back<typename solver::register_actions,
+                      observers::Actions::RegisterEventsWithObservers,
+                      Parallel::Actions::TerminatePhase>;
+
+  using step_actions = tmpl::list<Actions::RunEventsAndTriggers>;
+
+  using solve_actions =
+      tmpl::push_back<typename solver::template solve_actions<step_actions>,
+                      Parallel::Actions::TerminatePhase>;
+
+  using dg_element_array = elliptic::DgElementArray<
+      Metavariables,
+      tmpl::list<
+          Parallel::PhaseActions<Parallel::Phase::Initialization,
+                                 initialization_actions>,
+          Parallel::PhaseActions<Parallel::Phase::Register, register_actions>,
+          Parallel::PhaseActions<Parallel::Phase::Solve, solve_actions>>,
+      LinearSolver::multigrid::ElementsAllocator<
+          volume_dim, typename solver::multigrid::options_group>>;
+
+  // Specify all parallel components that will execute actions at some point.
+  using component_list = tmpl::flatten<
+      tmpl::list<dg_element_array, typename solver::component_list,
+                 observers::Observer<Metavariables>,
+                 observers::ObserverWriter<Metavariables>>>;
+
+  static constexpr std::array<Parallel::Phase, 4> default_phase_order{
+      {Parallel::Phase::Initialization, Parallel::Phase::Register,
+       Parallel::Phase::Solve, Parallel::Phase::Exit}};
+
+  // NOLINTNEXTLINE(google-runtime-references)
+  void pup(PUP::er& /*p*/) {}
+};
+
+static const std::vector<void (*)()> charm_init_node_funcs{
+    &setup_error_handling,
+    &setup_memory_allocation_failure_reporting,
+    &disable_openblas_multithreading,
+    &domain::creators::register_derived_with_charm,
+    &Parallel::register_derived_classes_with_charm<
+        metavariables::solver::schwarz_smoother::subdomain_solver>,
+    &elliptic::subdomain_preconditioners::register_derived_with_charm,
+    &Parallel::register_factory_classes_with_charm<metavariables>};
+static const std::vector<void (*)()> charm_init_proc_funcs{
+    &enable_floating_point_exceptions};
+/// \endcond

src/Elliptic/Executables/Punctures/SolvePuncturesFwd.hpp

@@ -0,0 +1,10 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+/// \cond
+struct Metavariables;
+/// \endcond

tests/InputFiles/Punctures/MultiplePunctures.yaml

@@ -0,0 +1,102 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+# Executable: SolvePunctures
+# Check: parse;execute
+# Timeout: 10
+# ExpectedOutput:
+#   PuncturesReductions.h5
+#   PuncturesVolume0.h5
+
+ResourceInfo:
+  AvoidGlobalProc0: false
+  Singletons: Auto
+
+Background:
+  MultiplePunctures:
+    Punctures:
+      - Position: [0., 0., 0.]
+        Mass: 1.
+        Momentum: [0., 0., 0.]
+        Spin: [0.1, 0.2, 0.3]
+
+InitialGuess: Flatness
+
+DomainCreator:
+  Sphere:
+    InnerRadius: 2.
+    OuterRadius: 10.
+    InnerCubeSphericity: 0.
+    InitialRefinement: 0
+    InitialGridPoints: [4, 4]
+    UseEquiangularMap: True
+    TimeDependence: None
+    BoundaryCondition: Flatness
+
+Discretization:
+  DiscontinuousGalerkin:
+    PenaltyParameter: 1.
+    Massive: True
+
+Observers:
+  VolumeFileName: "PuncturesVolume"
+  ReductionFileName: "PuncturesReductions"
+
+NonlinearSolver:
+  NewtonRaphson:
+    ConvergenceCriteria:
+      MaxIterations: 20
+      RelativeResidual: 0.
+      AbsoluteResidual: 1.e-10
+    SufficientDecrease: 1.e-4
+    MaxGlobalizationSteps: 40
+    DampingFactor: 1.
+    Verbosity: Quiet
+
+LinearSolver:
+  Gmres:
+    ConvergenceCriteria:
+      MaxIterations: 30
+      RelativeResidual: 1.e-4
+      AbsoluteResidual: 1.e-12
+    Verbosity: Quiet
+
+  Multigrid:
+    Iterations: 1
+    MaxLevels: Auto
+    PreSmoothing: True
+    PostSmoothingAtBottom: False
+    Verbosity: Silent
+    OutputVolumeData: False
+
+  SchwarzSmoother:
+    Iterations: 3
+    MaxOverlap: 2
+    Verbosity: Silent
+    SubdomainSolver:
+      Gmres:
+        ConvergenceCriteria:
+          MaxIterations: 3
+          RelativeResidual: 1.e-4
+          AbsoluteResidual: 1.e-12
+        Verbosity: Silent
+        Restart: None
+        Preconditioner:
+          MinusLaplacian:
+            Solver: ExplicitInverse
+            BoundaryConditions: Auto
+    SkipResets: True
+    ObservePerCoreReductions: False
+
+EventsAndTriggers:
+  - - HasConverged
+    - - ObserveFields:
+          SubfileName: VolumeData
+          VariablesToObserve:
+            - Field
+            - Alpha
+            - Beta
+          InterpolateToMesh: None
+          CoordinatesFloatingPointType: Double
+          FloatingPointTypes: [Double]
+          OverrideObservationValue: None