Add BNS initial data executable and SpEC initial Data

src/Elliptic/Actions/InitializeFixedSources.hpp

@@ -82,15 +82,18 @@ struct InitializeFixedSources : tt::ConformsTo<::amr::protocols::Projector> {
       domain::Tags::Coordinates<Dim, Frame::Inertial>, BackgroundTag,
       elliptic::dg::Tags::Massive, domain::Tags::Mesh<Dim>,
       domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>,
+      domain::Tags::InverseJacobian<3, Frame::ElementLogical, Frame::Inertial>,
       Parallel::Tags::Metavariables>;
 
   template <typename Background, typename Metavariables, typename... AmrData>
   static void apply(
       const gsl::not_null<typename fixed_sources_tag::type*> fixed_sources,
       const tnsr::I<DataVector, Dim>& inertial_coords,
       const Background& background, const bool massive, const Mesh<Dim>& mesh,
-      const Scalar<DataVector>& det_inv_jacobian, const Metavariables& /*meta*/,
-      const AmrData&... /*amr_data*/) {
+      const Scalar<DataVector>& det_inv_jacobian,
+      const InverseJacobian<DataVector, 3, Frame::ElementLogical,
+                            Frame::Inertial>& inv_jacobian,
+      const Metavariables& /*meta*/, const AmrData&... /*amr_data*/) {
     // Retrieve the fixed-sources of the elliptic system from the background,
     // which (along with the boundary conditions) define the problem we want to
     // solve.
@@ -99,9 +102,21 @@ struct InitializeFixedSources : tt::ConformsTo<::amr::protocols::Projector> {
     *fixed_sources =
         call_with_dynamic_type<Variables<typename fixed_sources_tag::tags_list>,
                                tmpl::at<factory_classes, Background>>(
-            &background, [&inertial_coords](const auto* const derived) {
-              return variables_from_tagged_tuple(derived->variables(
-                  inertial_coords, typename fixed_sources_tag::tags_list{}));
+            &background, [&inertial_coords, &mesh,
+                          &inv_jacobian](const auto* const derived) {
+              // Classes with background fields take the mesh and inverse
+              // Jacobian to be able to compute numerical derivatives
+              if constexpr (std::is_same_v<typename System::background_fields,
+                                           tmpl::list<>>) {
+                // Classes without `background_fields` do  not compute numerical
+                // derivatives
+                return variables_from_tagged_tuple(derived->variables(
+                    inertial_coords, typename fixed_sources_tag::tags_list{}));
+              } else {
+                return variables_from_tagged_tuple(derived->variables(
+                    inertial_coords, mesh, inv_jacobian,
+                    typename fixed_sources_tag::tags_list{}));
+              }
             });
 
     // Apply DG mass matrix to the fixed sources if the DG operator is massive

src/Elliptic/Executables/BnsInitialData/CMakeLists.txt

@@ -0,0 +1,43 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+set(EXECUTABLE "SolveBnsInitialData")
+
+add_spectre_executable(
+  ${EXECUTABLE}
+  EXCLUDE_FROM_ALL
+  SolveBnsInitialData.cpp
+)
+
+set(LIBS_TO_LINK
+  Charmxx::main
+  Convergence
+  CoordinateMaps
+  DiscontinuousGalerkin
+  DomainCreators
+  Elliptic
+  EllipticDg
+  EllipticDgSubdomainOperator
+  EllipticSubdomainPreconditioners
+  Events
+  EventsAndTriggers
+  Informer
+  BnsInitialData
+  BnsInitialDataAnalyticData
+  BnsInitialDataBoundaryConditions
+  LinearOperators
+  MathFunctions
+  Observer
+  Options
+  Parallel
+  ParallelLinearSolver
+  ParallelMultigrid
+  ParallelSchwarz
+  PhaseControl
+  RelativisticEulerSolutions
+  Utilities
+)
+
+
+
+target_link_libraries(${EXECUTABLE} PRIVATE ${LIBS_TO_LINK})

src/Elliptic/Executables/BnsInitialData/SolveBnsInitialData.cpp

@@ -0,0 +1,29 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Elliptic/Executables/BnsInitialData/SolveBnsInitialData.hpp"
+
+#include <vector>
+
+#include "Domain/Creators/RegisterDerivedWithCharm.hpp"
+#include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp"
+#include "Elliptic/SubdomainPreconditioners/RegisterDerived.hpp"
+#include "Parallel/CharmMain.tpp"
+#include "PointwiseFunctions/Hydro/EquationsOfState/RegisterDerivedWithCharm.hpp"
+#include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp"
+
+// Parameters chosen in CMakeLists.txt
+using metavariables = Metavariables;
+
+extern "C" void CkRegisterMainModule() {
+  Parallel::charmxx::register_main_module<metavariables>();
+  Parallel::charmxx::register_init_node_and_proc(
+      {&domain::creators::register_derived_with_charm,
+       &domain::FunctionsOfTime::register_derived_with_charm,
+       &register_derived_classes_with_charm<
+           metavariables::solver::schwarz_smoother::subdomain_solver>,
+       &elliptic::subdomain_preconditioners::register_derived_with_charm,
+       &EquationsOfState::register_derived_with_charm,
+       &register_factory_classes_with_charm<metavariables>},
+      {});
+}

src/Elliptic/Executables/BnsInitialData/SolveBnsInitialData.hpp

@@ -0,0 +1,185 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "DataStructures/DataBox/PrefixHelpers.hpp"
+#include "DataStructures/Tensor/IndexType.hpp"
+#include "Domain/Creators/Factory1D.hpp"
+#include "Domain/Creators/Factory2D.hpp"
+#include "Domain/Creators/Factory3D.hpp"
+#include "Domain/RadiallyCompressedCoordinates.hpp"
+#include "Domain/Tags.hpp"
+#include "Elliptic/Actions/RunEventsAndTriggers.hpp"
+#include "Elliptic/BoundaryConditions/BoundaryCondition.hpp"
+#include "Elliptic/DiscontinuousGalerkin/DgElementArray.hpp"
+#include "Elliptic/Executables/Solver.hpp"
+#include "Elliptic/Systems/BnsInitialData/BoundaryConditions/Factory.hpp"
+#include "Elliptic/Systems/BnsInitialData/FirstOrderSystem.hpp"
+#include "Elliptic/Triggers/Factory.hpp"
+#include "IO/Observer/Actions/RegisterEvents.hpp"
+#include "IO/Observer/Helpers.hpp"
+#include "IO/Observer/ObserverComponent.hpp"
+#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
+#include "Options/Protocols/FactoryCreation.hpp"
+#include "Options/String.hpp"
+#include "Parallel/Phase.hpp"
+#include "Parallel/PhaseControl/ExecutePhaseChange.hpp"
+#include "Parallel/PhaseControl/VisitAndReturn.hpp"
+#include "Parallel/PhaseDependentActionList.hpp"
+#include "Parallel/Protocols/RegistrationMetavariables.hpp"
+#include "Parallel/Reduction.hpp"
+#include "ParallelAlgorithms/Actions/TerminatePhase.hpp"
+#include "ParallelAlgorithms/Amr/Actions/SendAmrDiagnostics.hpp"
+#include "ParallelAlgorithms/Amr/Criteria/Factory.hpp"
+#include "ParallelAlgorithms/Amr/Protocols/AmrMetavariables.hpp"
+#include "ParallelAlgorithms/Events/Factory.hpp"
+#include "ParallelAlgorithms/Events/Tags.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/ObserveVolumeData.hpp"
+#include "ParallelAlgorithms/LinearSolver/Multigrid/Tags.hpp"
+#include "PointwiseFunctions/AnalyticData/BnsInitialData/Factory.hpp"
+#include "PointwiseFunctions/Hydro/EquationsOfState/RegisterDerivedWithCharm.hpp"
+#include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/Background.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/InitialGuess.hpp"
+#include "PointwiseFunctions/MathFunctions/Factory.hpp"
+#include "Utilities/ProtocolHelpers.hpp"
+#include "Utilities/TMPL.hpp"
+
+/// \cond
+namespace PUP {
+class er;
+}  // namespace PUP
+struct Metavariables {
+  static constexpr Options::String help{
+      "Find the solution for the velocity potential of"
+      "an irrotational BNS system using fixed"
+      "problem on fixed space-time and with a fixed"
+      "rest-mass density (alternatively fixed "
+      "specific enthalpy) profile."};
+  static constexpr size_t volume_dim = 3;
+  using system = BnsInitialData::FirstOrderSystem;
+  using solver = elliptic::Solver<Metavariables>;
+
+  using solved_fields = typename system::primal_fields;
+  using deriv_fields = tmpl::list<::Tags::DerivTensorCompute<
+      tmpl::front<solved_fields>,
+      domain::Tags::InverseJacobian<3, Frame::ElementLogical, Frame::Inertial>,
+      domain::Tags::Mesh<3>>>;
+  using observe_fields = tmpl::append<
+      solved_fields, typename system::background_fields,
+      typename solver::observe_fields, deriv_fields,
+      tmpl::list<domain::Tags::Coordinates<volume_dim, Frame::Inertial>,
+                 ::Events::Tags::ObserverDetInvJacobianCompute<
+                     Frame::ElementLogical, Frame::Inertial>,
+                 domain::Tags::RadiallyCompressedCoordinatesCompute<
+                     volume_dim, Frame::Inertial>>>;
+  using observer_compute_tags =
+      tmpl::list<::Events::Tags::ObserverMeshCompute<volume_dim>>;
+
+  // Collect all items to store in the cache.
+  using const_global_cache_tags =
+      tmpl::list<domain::Tags::RadiallyCompressedCoordinatesOptions,
+                 BnsInitialData::Tags::EulerEnthalpyConstant>;
+  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,
+                   BnsInitialData::InitialData::all_initial_data>,
+        tmpl::pair<elliptic::analytic_data::InitialGuess,
+                   BnsInitialData::InitialData::all_initial_data>,
+        tmpl::pair<elliptic::analytic_data::AnalyticSolution, tmpl::list<>>,
+        tmpl::pair<elliptic::BoundaryConditions::BoundaryCondition<volume_dim>,
+                   BnsInitialData::BoundaryConditions::
+                       standard_boundary_conditions<system>>,
+        tmpl::pair<
+            ::amr::Criterion,
+            ::amr::Criteria::standard_criteria<
+                volume_dim, tmpl::list<BnsInitialData::Tags::VelocityPotential<
+                                DataVector>>>>,
+        tmpl::pair<Trigger, elliptic::Triggers::all_triggers<
+                                ::amr::OptionTags::AmrGroup>>,
+        tmpl::pair<Event,
+                   tmpl::flatten<tmpl::list<
+                       Events::Completion,
+                       dg::Events::field_observations<
+                           volume_dim, observe_fields, observer_compute_tags,
+                           LinearSolver::multigrid::Tags::IsFinestGrid>>>>,
+        tmpl::pair<
+            PhaseChange,
+            tmpl::list<
+                PhaseControl::VisitAndReturn<Parallel::Phase::BuildMatrix>,
+                // Phases for AMR
+                PhaseControl::VisitAndReturn<
+                    Parallel::Phase::EvaluateAmrCriteria>,
+                PhaseControl::VisitAndReturn<Parallel::Phase::AdjustDomain>,
+                PhaseControl::VisitAndReturn<Parallel::Phase::CheckDomain>>>>;
+  };
+
+  // Collect all reduction tags for observers
+  using observed_reduction_data_tags =
+      observers::collect_reduction_data_tags<tmpl::flatten<tmpl::list<
+          tmpl::at<typename 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>;
+
+  using solve_actions = typename solver::template solve_actions<tmpl::list<>>;
+
+  using dg_element_array = elliptic::DgElementArray<
+      Metavariables,
+      tmpl::list<
+          Parallel::PhaseActions<Parallel::Phase::Initialization,
+                                 initialization_actions>,
+          Parallel::PhaseActions<
+              Parallel::Phase::Register,
+              tmpl::push_back<register_actions,
+                              Parallel::Actions::TerminatePhase>>,
+          Parallel::PhaseActions<Parallel::Phase::Solve, solve_actions>,
+          Parallel::PhaseActions<Parallel::Phase::CheckDomain,
+                                 tmpl::list<::amr::Actions::SendAmrDiagnostics,
+                                            Parallel::Actions::TerminatePhase>>,
+          Parallel::PhaseActions<
+              Parallel::Phase::BuildMatrix,
+              tmpl::push_back<typename solver::build_matrix_actions,
+                              Parallel::Actions::TerminatePhase>>>,
+      LinearSolver::multigrid::ElementsAllocator<
+          volume_dim, typename solver::multigrid::options_group>>;
+
+  struct amr : tt::ConformsTo<::amr::protocols::AmrMetavariables> {
+    using element_array = dg_element_array;
+    using projectors = typename solver::amr_projectors;
+  };
+  struct registration
+      : tt::ConformsTo<Parallel::protocols::RegistrationMetavariables> {
+    using element_registrars =
+        tmpl::map<tmpl::pair<dg_element_array, register_actions>>;
+  };
+
+  // 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*/) {}
+};
+/// \endcond

src/Elliptic/Executables/CMakeLists.txt

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

src/Elliptic/Systems/BnsInitialData/BoundaryConditions/CMakeLists.txt

@@ -1,7 +1,7 @@
 # Distributed under the MIT License.
 # See LICENSE.txt for details.
 
-set(LIBRARY IrrotationalBnsBoundaryConditions)
+set(LIBRARY BnsInitialDataBoundaryConditions)
 
 add_spectre_library(${LIBRARY})
 

src/Elliptic/Systems/BnsInitialData/Equations.hpp

@@ -77,6 +77,8 @@ struct Fluxes {
                  Tags::RotationalShiftStress<DataVector>>;
   using volume_tags = tmpl::list<>;
   using const_global_cache_tags = tmpl::list<>;
+  static constexpr bool is_trivial = false;
+  static constexpr bool is_discontinuous = false;
   static void apply(gsl::not_null<tnsr::I<DataVector, 3>*> flux_for_potential,
                     const tnsr::II<DataVector, 3>& inverse_spatial_metric,
                     const tnsr::II<DataVector, 3>& rotational_shift_stress,

src/Elliptic/Systems/BnsInitialData/FirstOrderSystem.hpp

@@ -14,6 +14,7 @@
 #include "Elliptic/Systems/BnsInitialData/Tags.hpp"
 #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
 #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
+#include "PointwiseFunctions/Hydro/Tags.hpp"
 #include "Utilities/ProtocolHelpers.hpp"
 #include "Utilities/TMPL.hpp"
 
@@ -76,6 +77,7 @@ struct FirstOrderSystem
       ::Tags::Flux<velocity_potential, tmpl::size_t<3>, Frame::Inertial>>;
 
   using background_fields = tmpl::list<
+      hydro::Tags::RestMassDensity<DataVector>,
       gr::Tags::InverseSpatialMetric<DataVector, 3>,
       gr::Tags::SpatialChristoffelSecondKindContracted<DataVector, 3>,
       gr::Tags::Lapse<DataVector>,

src/Elliptic/Systems/BnsInitialData/Tags.hpp

@@ -7,6 +7,8 @@
 
 #include "DataStructures/DataBox/Tag.hpp"
 #include "DataStructures/Tensor/TypeAliases.hpp"
+#include "Options/Options.hpp"
+#include "Options/String.hpp"
 
 /// \cond
 class DataVector;
@@ -17,6 +19,14 @@ class DataVector;
  * \brief Items related to solving for irrotational bns initial data
  */
 namespace BnsInitialData::Tags {
+namespace OptionTags {
+struct EulerEnthalpyConstant {
+  using type = double;
+  static constexpr Options::String help =
+      "The Euler Enthalpy constant of the star";
+};
+
+}  // namespace OptionTags
 
 /*!
  * \brief The shift plus a spatial vector \f$ k^i\f$
@@ -65,6 +75,9 @@ struct DerivSpatialRotationalKillingVector : db::SimpleTag {
 
 struct EulerEnthalpyConstant : db::SimpleTag {
   using type = double;
+  using option_tags = tmpl::list<OptionTags::EulerEnthalpyConstant>;
+  static constexpr bool pass_metavariables = false;
+  static double create_from_options(const double value) { return value; };
 };
 
 }  // namespace BnsInitialData::Tags