diff --git a/src/Elliptic/Actions/InitializeFixedSources.hpp b/src/Elliptic/Actions/InitializeFixedSources.hpp index dd2227d669f2..dcfc0b1e3ffe 100644 --- a/src/Elliptic/Actions/InitializeFixedSources.hpp +++ b/src/Elliptic/Actions/InitializeFixedSources.hpp @@ -82,6 +82,7 @@ struct InitializeFixedSources : tt::ConformsTo<::amr::protocols::Projector> { domain::Tags::Coordinates, BackgroundTag, elliptic::dg::Tags::Massive, domain::Tags::Mesh, domain::Tags::DetInvJacobian, + domain::Tags::InverseJacobian<3, Frame::ElementLogical, Frame::Inertial>, Parallel::Tags::Metavariables>; template @@ -89,8 +90,10 @@ struct InitializeFixedSources : tt::ConformsTo<::amr::protocols::Projector> { const gsl::not_null fixed_sources, const tnsr::I& inertial_coords, const Background& background, const bool massive, const Mesh& mesh, - const Scalar& det_inv_jacobian, const Metavariables& /*meta*/, - const AmrData&... /*amr_data*/) { + const Scalar& det_inv_jacobian, + const InverseJacobian& 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, tmpl::at>( - &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>) { + // 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 diff --git a/src/Elliptic/Executables/BnsInitialData/CMakeLists.txt b/src/Elliptic/Executables/BnsInitialData/CMakeLists.txt new file mode 100644 index 000000000000..c3d5b9b3b9c8 --- /dev/null +++ b/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}) diff --git a/src/Elliptic/Executables/BnsInitialData/SolveBnsInitialData.cpp b/src/Elliptic/Executables/BnsInitialData/SolveBnsInitialData.cpp new file mode 100644 index 000000000000..cbb1925296a7 --- /dev/null +++ b/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 + +#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(); + Parallel::charmxx::register_init_node_and_proc( + {&domain::creators::register_derived_with_charm, + &domain::FunctionsOfTime::register_derived_with_charm, + ®ister_derived_classes_with_charm< + metavariables::solver::schwarz_smoother::subdomain_solver>, + &elliptic::subdomain_preconditioners::register_derived_with_charm, + &EquationsOfState::register_derived_with_charm, + ®ister_factory_classes_with_charm}, + {}); +} diff --git a/src/Elliptic/Executables/BnsInitialData/SolveBnsInitialData.hpp b/src/Elliptic/Executables/BnsInitialData/SolveBnsInitialData.hpp new file mode 100644 index 000000000000..4115f58b28a2 --- /dev/null +++ b/src/Elliptic/Executables/BnsInitialData/SolveBnsInitialData.hpp @@ -0,0 +1,185 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#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; + + using solved_fields = typename system::primal_fields; + using deriv_fields = tmpl::list<::Tags::DerivTensorCompute< + tmpl::front, + 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, + ::Events::Tags::ObserverDetInvJacobianCompute< + Frame::ElementLogical, Frame::Inertial>, + domain::Tags::RadiallyCompressedCoordinatesCompute< + volume_dim, Frame::Inertial>>>; + using observer_compute_tags = + tmpl::list<::Events::Tags::ObserverMeshCompute>; + + // Collect all items to store in the cache. + using const_global_cache_tags = + tmpl::list; + struct factory_creation + : tt::ConformsTo { + using factory_classes = tmpl::map< + tmpl::pair, domain_creators>, + tmpl::pair, + tmpl::pair, + tmpl::pair>, + tmpl::pair, + BnsInitialData::BoundaryConditions:: + standard_boundary_conditions>, + tmpl::pair< + ::amr::Criterion, + ::amr::Criteria::standard_criteria< + volume_dim, tmpl::list>>>, + tmpl::pair>, + tmpl::pair>>>, + tmpl::pair< + PhaseChange, + tmpl::list< + PhaseControl::VisitAndReturn, + // Phases for AMR + PhaseControl::VisitAndReturn< + Parallel::Phase::EvaluateAmrCriteria>, + PhaseControl::VisitAndReturn, + PhaseControl::VisitAndReturn>>>; + }; + + // Collect all reduction tags for observers + using observed_reduction_data_tags = + observers::collect_reduction_data_tags, + solver>>>; + + using initialization_actions = + tmpl::push_back; + + using register_actions = + tmpl::push_back; + + using solve_actions = typename solver::template solve_actions>; + + using dg_element_array = elliptic::DgElementArray< + Metavariables, + tmpl::list< + Parallel::PhaseActions, + Parallel::PhaseActions< + Parallel::Phase::Register, + tmpl::push_back>, + Parallel::PhaseActions, + Parallel::PhaseActions>, + Parallel::PhaseActions< + Parallel::Phase::BuildMatrix, + tmpl::push_back>>, + 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 { + using element_registrars = + tmpl::map>; + }; + + // Specify all parallel components that will execute actions at some point. + using component_list = tmpl::flatten< + tmpl::list, + observers::ObserverWriter>>; + + static constexpr std::array 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 diff --git a/src/Elliptic/Executables/CMakeLists.txt b/src/Elliptic/Executables/CMakeLists.txt index 1d414f05ba03..144d93794459 100644 --- a/src/Elliptic/Executables/CMakeLists.txt +++ b/src/Elliptic/Executables/CMakeLists.txt @@ -5,3 +5,4 @@ add_subdirectory(Elasticity) add_subdirectory(Poisson) add_subdirectory(Punctures) add_subdirectory(Xcts) +add_subdirectory(BnsInitialData) diff --git a/src/Elliptic/Systems/BnsInitialData/BoundaryConditions/CMakeLists.txt b/src/Elliptic/Systems/BnsInitialData/BoundaryConditions/CMakeLists.txt index 5c2d5fe67dd4..fc419b9b65f6 100644 --- a/src/Elliptic/Systems/BnsInitialData/BoundaryConditions/CMakeLists.txt +++ b/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}) diff --git a/src/Elliptic/Systems/BnsInitialData/Equations.hpp b/src/Elliptic/Systems/BnsInitialData/Equations.hpp index 80b0ff950f97..4d797f625d28 100644 --- a/src/Elliptic/Systems/BnsInitialData/Equations.hpp +++ b/src/Elliptic/Systems/BnsInitialData/Equations.hpp @@ -77,6 +77,8 @@ struct Fluxes { Tags::RotationalShiftStress>; 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*> flux_for_potential, const tnsr::II& inverse_spatial_metric, const tnsr::II& rotational_shift_stress, diff --git a/src/Elliptic/Systems/BnsInitialData/FirstOrderSystem.hpp b/src/Elliptic/Systems/BnsInitialData/FirstOrderSystem.hpp index 59b11e3449b5..d16a810a9f2a 100644 --- a/src/Elliptic/Systems/BnsInitialData/FirstOrderSystem.hpp +++ b/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, Frame::Inertial>>; using background_fields = tmpl::list< + hydro::Tags::RestMassDensity, gr::Tags::InverseSpatialMetric, gr::Tags::SpatialChristoffelSecondKindContracted, gr::Tags::Lapse, diff --git a/src/Elliptic/Systems/BnsInitialData/Tags.hpp b/src/Elliptic/Systems/BnsInitialData/Tags.hpp index c381d1e7c432..48f49174eaa0 100644 --- a/src/Elliptic/Systems/BnsInitialData/Tags.hpp +++ b/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; + static constexpr bool pass_metavariables = false; + static double create_from_options(const double value) { return value; }; }; } // namespace BnsInitialData::Tags diff --git a/src/PointwiseFunctions/AnalyticData/BnsInitialData/CMakeLists.txt b/src/PointwiseFunctions/AnalyticData/BnsInitialData/CMakeLists.txt new file mode 100644 index 000000000000..72691f86e04d --- /dev/null +++ b/src/PointwiseFunctions/AnalyticData/BnsInitialData/CMakeLists.txt @@ -0,0 +1,36 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. +set(LIBRARY BnsInitialDataAnalyticData) + +add_spectre_library(${LIBRARY}) + + +target_link_libraries( + ${LIBRARY} + PUBLIC + DataStructures + ErrorHandling + GeneralRelativity + Options + RelativisticEulerSolutions + Utilities + ) + if (TARGET SpEC::Exporter) + spectre_target_sources( + ${LIBRARY} + PRIVATE + SpecInitialData.cpp + ) + spectre_target_headers( + ${LIBRARY} + INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src + HEADERS + SpecInitialData.hpp + ) + target_link_libraries( + ${LIBRARY} + PUBLIC + GeneralRelativityAnalyticData + SpEC::Exporter + ) + endif() diff --git a/src/PointwiseFunctions/AnalyticData/BnsInitialData/Factory.hpp b/src/PointwiseFunctions/AnalyticData/BnsInitialData/Factory.hpp new file mode 100644 index 000000000000..48b01cd9102a --- /dev/null +++ b/src/PointwiseFunctions/AnalyticData/BnsInitialData/Factory.hpp @@ -0,0 +1,10 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include "PointwiseFunctions/AnalyticData/BnsInitialData/SpecInitialData.hpp" +namespace BnsInitialData::InitialData { +using all_initial_data = tmpl::list>; + +} // namespace BnsInitialData::InitialData diff --git a/src/PointwiseFunctions/AnalyticData/BnsInitialData/SpecInitialData.cpp b/src/PointwiseFunctions/AnalyticData/BnsInitialData/SpecInitialData.cpp new file mode 100644 index 000000000000..7835ab52da92 --- /dev/null +++ b/src/PointwiseFunctions/AnalyticData/BnsInitialData/SpecInitialData.cpp @@ -0,0 +1,126 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "PointwiseFunctions/AnalyticData/BnsInitialData/SpecInitialData.hpp" + +#include // The SpEC Exporter +#include +#include +#include +#include +#include + +#include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp" +#include "DataStructures/Tensor/EagerMath/DotProduct.hpp" +#include "DataStructures/Tensor/EagerMath/RaiseOrLowerIndex.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "IO/External/InterpolateFromSpec.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "PointwiseFunctions/Hydro/SpecificEnthalpy.hpp" +#include "PointwiseFunctions/Hydro/Tags.hpp" +#include "Utilities/ContainerHelpers.hpp" +#include "Utilities/ErrorHandling/Error.hpp" +#include "Utilities/GenerateInstantiations.hpp" +#include "Utilities/System/ParallelInfo.hpp" +#include "Utilities/TaggedTuple.hpp" + +namespace BnsInitialData::InitialData { + +template +SpecInitialData::SpecInitialData( + std::string data_directory, + std::unique_ptr equation_of_state, + const double density_cutoff, const double orbital_angular_velocity, + const double euler_enthalpy_constant) + : data_directory_(std::move(data_directory)), + equation_of_state_(std::move(equation_of_state)), + density_cutoff_(density_cutoff), + orbital_angular_velocity_(orbital_angular_velocity), + euler_enthalpy_constant_(euler_enthalpy_constant), + spec_exporter_(std::make_unique( + sys::procs_on_node(sys::my_node()), data_directory_, + vars_to_interpolate_)) {} + +template +SpecInitialData& SpecInitialData::operator=( + const SpecInitialData& rhs) { + data_directory_ = rhs.data_directory_; + equation_of_state_ = rhs.equation_of_state_->get_clone(); + density_cutoff_ = rhs.density_cutoff_; + orbital_angular_velocity_ = rhs.orbital_angular_velocity_; + euler_enthalpy_constant_ = rhs.euler_enthalpy_constant_; + spec_exporter_ = + std::make_unique(sys::procs_on_node(sys::my_node()), + data_directory_, vars_to_interpolate_); + return *this; +} + +template +SpecInitialData::SpecInitialData(const SpecInitialData& rhs) { + data_directory_ = rhs.data_directory_; + equation_of_state_ = rhs.equation_of_state_->get_clone(); + density_cutoff_ = rhs.density_cutoff_; + orbital_angular_velocity_ = rhs.orbital_angular_velocity_; + euler_enthalpy_constant_ = rhs.euler_enthalpy_constant_; + spec_exporter_ = + std::make_unique(sys::procs_on_node(sys::my_node()), + data_directory_, vars_to_interpolate_); +} + +template +std::unique_ptr +SpecInitialData::get_clone() const { + return std::make_unique(*this); +} + +template +SpecInitialData::SpecInitialData(CkMigrateMessage* msg) + : elliptic::analytic_data::Background(msg) {} + +template +void SpecInitialData::pup(PUP::er& p) { + elliptic::analytic_data::Background::pup(p); + p | data_directory_; + p | equation_of_state_; + p | density_cutoff_; + p | orbital_angular_velocity_; + p | euler_enthalpy_constant_; + if (p.isUnpacking()) { + spec_exporter_ = + std::make_unique(sys::procs_on_node(sys::my_node()), + data_directory_, vars_to_interpolate_); + } +} + +template +PUP::able::PUP_ID SpecInitialData::my_PUP_ID = 0; + +template +template +tuples::tagged_tuple_from_typelist::template interpolated_tags> +SpecInitialData::interpolate_from_spec( + const tnsr::I& x) const { + return io::interpolate_from_spec>( + make_not_null(spec_exporter_.get()), x, + static_cast(sys::my_local_rank())); +} + +#define THERMODIM(data) BOOST_PP_TUPLE_ELEM(0, data) + +#define INSTANTIATION(r, data) \ + template class SpecInitialData; \ + template tuples::tagged_tuple_from_typelist::template interpolated_tags> \ + SpecInitialData::interpolate_from_spec( \ + const tnsr::I& x) const; \ + template tuples::tagged_tuple_from_typelist::template interpolated_tags> \ + SpecInitialData::interpolate_from_spec( \ + const tnsr::I& x) const; + +GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) + +#undef INSTANTIATION +#undef THERMODIM +} // namespace BnsInitialData::InitialData diff --git a/src/PointwiseFunctions/AnalyticData/BnsInitialData/SpecInitialData.hpp b/src/PointwiseFunctions/AnalyticData/BnsInitialData/SpecInitialData.hpp new file mode 100644 index 000000000000..1bfb4b125cda --- /dev/null +++ b/src/PointwiseFunctions/AnalyticData/BnsInitialData/SpecInitialData.hpp @@ -0,0 +1,335 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include // The SpEC Exporter +#include +#include + +#include "DataStructures/CachedTempBuffer.hpp" +#include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp" +#include "Elliptic/Systems/BnsInitialData/FirstOrderSystem.hpp" +#include "Elliptic/Systems/BnsInitialData/Tags.hpp" +#include "Evolution/NumericInitialData.hpp" +#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp" +#include "Options/String.hpp" +#include "PointwiseFunctions/AnalyticData/BnsInitialData/TovStar.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp" +#include "PointwiseFunctions/Hydro/InitialData/IrrotationalBns.hpp" +#include "PointwiseFunctions/Hydro/Tags.hpp" +#include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp" +#include "Utilities/Serialization/CharmPupable.hpp" +#include "Utilities/TMPL.hpp" +#include "Utilities/TaggedTuple.hpp" + +/// \cond +namespace PUP { +class er; +} // namespace PUP +/// \endcond + +namespace BnsInitialData::InitialData { + +/*! + * \brief Hydro initial data generated by SpEC. + * + * This class loads numerical data written out by the SpEC initial data solver. + * It uses the `spec::Exporter` linked in from SpEC to interpolate to arbitrary + * grid points. The coordinates are assumed to be in SpEC's "grid" frame. + * We interpolate the following quantities: + * + * - "g": spatial metric + * - "K": (lower) extrinsic curvature + * - "Lapse": lapse + * - "Shift": (upper) shift + * - "BaryonDensity": rest mass density + * - "u_i": lower spatial four-velocity + * + * The remaining hydro quantities are computed from the interpolated data and + * the equation of state. + * The magnetic field is set to zero and the electron fraction is set to a + * constant read from the input file. + */ +template +class SpecInitialData : public elliptic::analytic_data::Background, + public elliptic::analytic_data::InitialGuess, + public evolution::NumericInitialData { + public: + using equation_of_state_type = + EquationsOfState::EquationOfState; + + template + using variable_tags = + tmpl::list>; + template + using background_tags = tmpl::list< + hydro::Tags::RestMassDensity, + gr::Tags::InverseSpatialMetric, + gr::Tags::SpatialChristoffelSecondKindContracted, + gr::Tags::Lapse, + ::Tags::deriv, + tmpl::integral_constant, Frame::Inertial>, + gr::Tags::Shift, + ::Tags::deriv, + tmpl::integral_constant, Frame::Inertial>, + Tags::RotationalShift, + Tags::DerivLogLapseTimesDensityOverSpecificEnthalpy, + Tags::RotationalShiftStress>; + template + using tags = tmpl::append, variable_tags>; + + struct DataDirectory { + using type = std::string; + static constexpr Options::String help = { + "Path to a directory of data produced by SpEC. The directory is " + "expected to contain 'GrDomain.input' and 'Vars*.h5' files for all the " + "subdomains in GrDomain.input."}; + }; + + struct DensityCutoff { + using type = double; + static constexpr Options::String help = + "Where the density is below this cutoff the fluid variables are set to " + "vacuum (zero density, pressure, energy and velocity, unit Lorentz " + "factor and enthalpy). " + "During the evolution, atmosphere treatment will typically kick in and " + "fix the value of the fluid variables in these regions. Therefore, " + "it makes sense to set this density cutoff to the same value as the " + "atmosphere density cutoff."; + }; + + using options = tmpl::list< + DataDirectory, + hydro::OptionTags::InitialDataEquationOfState, + DensityCutoff, + BnsInitialData::InitialData::TovStar::OrbitalAngularVelocity, + BnsInitialData::Tags::OptionTags::EulerEnthalpyConstant>; + + static constexpr Options::String help = {"Initial data generated by SpEC"}; + + SpecInitialData() = default; + SpecInitialData(const SpecInitialData& rhs); + SpecInitialData& operator=(const SpecInitialData& rhs); + SpecInitialData(SpecInitialData&& /*rhs*/) = default; + SpecInitialData& operator=(SpecInitialData&& /*rhs*/) = default; + ~SpecInitialData() override = default; + + SpecInitialData(std::string data_directory, + std::unique_ptr equation_of_state, + double density_cutoff, double orbital_angular_velocity, + double euler_enthalpy_constant); + + auto get_clone() const + -> std::unique_ptr; + + /// \cond + explicit SpecInitialData(CkMigrateMessage* msg); + using PUP::able::register_constructor; + WRAPPED_PUPable_decl_template(SpecInitialData); + /// \endcond + const equation_of_state_type& equation_of_state() const { + return *equation_of_state_; + } + // The velocity potential is used in the initial guess + template + tuples::TaggedTuple> variables( + const tnsr::I& x, + tmpl::list> /*meta*/) const { + // return velocity potential (only a guess) + Scalar velocity_potential = + make_with_value>(x, orbital_angular_velocity_); + // This is not a good guess, but it's better than nothing. + // It's not clear to me this should actually be used + get(velocity_potential) *= (get<1>(x) * get<0>(x)); + + tuples::TaggedTuple> result{ + velocity_potential}; + return result; + } + // The fixed sources are used in initialization + template + tuples::TaggedTuple<::Tags::FixedSource>> + variables(const tnsr::I& x, const Mesh<3>& mesh, + const InverseJacobian& inv_jacobian, + tmpl::list<::Tags::FixedSource< + Tags::VelocityPotential>> /*meta*/) const { + auto background_values = + variables(x, mesh, inv_jacobian, background_tags{}); + tuples::TaggedTuple<::Tags::FixedSource>> + result{}; + const auto& lapse = get>(background_values); + const auto& rotational_shift = + get>(background_values); + const auto& deriv_log_lapse_times_density_over_specific_enthalpy = + get>( + background_values); + const auto& deriv_of_lapse = + get<::Tags::deriv, + tmpl::integral_constant, Frame::Inertial>>( + background_values); + + const auto& deriv_of_shift = + get<::Tags::deriv, + tmpl::integral_constant, Frame::Inertial>>( + background_values); + ::tenex::evaluate<>( + make_not_null( + &get<::Tags::FixedSource>>( + result)), + euler_enthalpy_constant_ * + (1.0 / square(lapse()) * rotational_shift(ti::I) * + deriv_log_lapse_times_density_over_specific_enthalpy(ti::i) - + 2.0 / square(lapse()) * rotational_shift(ti::I) * + deriv_of_lapse(ti::i) + + 1.0 / square(lapse()) * deriv_of_shift(ti::i, ti::I))); + + return result; + } + template + tuples::TaggedTuple> variables( + const tnsr::I& x, + tmpl::list> /*meta*/) const { + // interpolate from spec, then set gamma + auto interpolated_vars = interpolate_from_spec(x); + + auto spatial_metric = + get>(interpolated_vars); + tuples::TaggedTuple> result{}; + get>(result) = + determinant_and_inverse(spatial_metric).second; + return result; + } + template + tuples::tagged_tuple_from_typelist> variables( + const tnsr::I& x, const Mesh<3>& mesh, + const InverseJacobian& inv_jacobian, + background_tags /*meta*/) const { + // interpolate from spec, take num derivatives, return + // Shift, lapse spatial metric imported + auto result = + tuples::tagged_tuple_from_typelist>{}; + auto interpolated_vars = interpolate_from_spec(x); + auto spatial_metric = + get>(interpolated_vars); + auto spatial_metric_determinant_and_inverse = + determinant_and_inverse(spatial_metric); + auto inv_spatial_metric = spatial_metric_determinant_and_inverse.second; + get>(result) = + inv_spatial_metric; + auto sqrt_det_spatial_metric = Scalar{ + sqrt(get(spatial_metric_determinant_and_inverse.first))}; + // Get the one contracted Christoffel needed for fluxes (is there a more + // uniform way to do this?) + const auto spatial_christoffel_second_kind_contracted = + tenex::evaluate( + 0.5 * + partial_derivative( + get>(interpolated_vars), + mesh, inv_jacobian)(ti::i, ti::j, ti::k) * + inv_spatial_metric(ti::J, ti::K)); + get>(result) = + spatial_christoffel_second_kind_contracted; + get>(result) = + get>(interpolated_vars); + + // Get Lapse and shift derivatives + get<::Tags::deriv, + tmpl::integral_constant, Frame::Inertial>>( + result) = + partial_derivative(get>(interpolated_vars), + mesh, inv_jacobian); + get>(result) = + get>(interpolated_vars); + get<::Tags::deriv, + tmpl::integral_constant, Frame::Inertial>>( + result) = + partial_derivative(get>(interpolated_vars), + mesh, inv_jacobian); + // Get the rotational shift + deriv of log lapse over enthalpy + stress + const auto spatial_rotational_killing_vector = hydro::initial_data:: + irrotational_bns::spatial_rotational_killing_vector( + x, orbital_angular_velocity_, sqrt_det_spatial_metric); + const auto rotational_shift = + hydro::initial_data::irrotational_bns::rotational_shift( + get>(interpolated_vars), + spatial_rotational_killing_vector); + get>(result) = rotational_shift; + auto& rest_mass_density = + get>(interpolated_vars); + get>(result) = rest_mass_density; + // The SpEC solution should have e.g. B&S Eq. 15.76 satisfied + // so this is equivalent to some expression in the other hydro variables + get(rest_mass_density) += 1.0e-15; + const DataType enthalpy_density = + get(equation_of_state_->pressure_from_density(rest_mass_density)) + + get(rest_mass_density) * + (1.0 + + get(equation_of_state_->specific_internal_energy_from_density( + rest_mass_density))); + const auto deriv_log_lapse_times_density_over_specific_enthalpy = + partial_derivative( + Scalar{ + get(get>(interpolated_vars)) * + square(get(rest_mass_density)) / enthalpy_density}, + mesh, inv_jacobian); + get>(result) = + deriv_log_lapse_times_density_over_specific_enthalpy; + const auto rotational_shift_stress = + hydro::initial_data::irrotational_bns::rotational_shift_stress( + rotational_shift, + get>(interpolated_vars)); + get>(result) = + rotational_shift_stress; + // Compute fixed sources + return result; + } + + // NOLINTNEXTLINE(google-runtime-references) + void pup(PUP::er& /*p*/) override; + + private: + /// These quantities are supported for interpolation from SpEC + template + using interpolated_tags = tmpl::list< + // GR quantities + gr::Tags::SpatialMetric, + gr::Tags::ExtrinsicCurvature, gr::Tags::Lapse, + gr::Tags::Shift, + // Hydro quantities + hydro::Tags::RestMassDensity, + hydro::Tags::LowerSpatialFourVelocity>; + + /// These are the names in SpEC datasets corresponding to the quantities above + static const inline std::vector vars_to_interpolate_{ + // GR quantities + "g", "K", "Lapse", "Shift", + // Hydro quantities + "BaryonDensity", "u_i"}; + + template + tuples::tagged_tuple_from_typelist> + interpolate_from_spec(const tnsr::I& x) const; + + /// This cache computes all derived quantities from the interpolated + /// quantities on demand + template + using VariablesCache = cached_temp_buffer_from_typelist, interpolated_tags>, + hydro::Tags::LorentzFactorTimesSpatialVelocity, + gr::Tags::InverseSpatialMetric>>; + + std::string data_directory_{}; + std::unique_ptr equation_of_state_{nullptr}; + double density_cutoff_ = std::numeric_limits::signaling_NaN(); + double orbital_angular_velocity_ = + std::numeric_limits::signaling_NaN(); + double euler_enthalpy_constant_ = + std::numeric_limits::signaling_NaN(); + std::unique_ptr spec_exporter_{nullptr}; +}; + +} // namespace BnsInitialData::InitialData diff --git a/src/PointwiseFunctions/AnalyticData/CMakeLists.txt b/src/PointwiseFunctions/AnalyticData/CMakeLists.txt index d2fb24d21604..fe600ebb4253 100644 --- a/src/PointwiseFunctions/AnalyticData/CMakeLists.txt +++ b/src/PointwiseFunctions/AnalyticData/CMakeLists.txt @@ -14,6 +14,7 @@ spectre_target_headers( ) add_subdirectory(Burgers) +add_subdirectory(BnsInitialData) add_subdirectory(CurvedWaveEquation) add_subdirectory(ForceFree) add_subdirectory(GeneralRelativity) diff --git a/src/PointwiseFunctions/Hydro/InitialData/IrrotationalBns.cpp b/src/PointwiseFunctions/Hydro/InitialData/IrrotationalBns.cpp index beee635c4a62..72e5f1df5a16 100644 --- a/src/PointwiseFunctions/Hydro/InitialData/IrrotationalBns.cpp +++ b/src/PointwiseFunctions/Hydro/InitialData/IrrotationalBns.cpp @@ -116,10 +116,9 @@ void spatial_rotational_killing_vector( const Scalar& sqrt_det_spatial_metric) { // Cross product involves volume element in arbitrary coordinates set_number_of_grid_points(result, sqrt_det_spatial_metric); - get<0>(*result) = - -get(sqrt_det_spatial_metric) * get<1>(x) * orbital_angular_velocity; - get<1>(*result) = - get(sqrt_det_spatial_metric) * get<0>(x) * orbital_angular_velocity; + // Should this get a factor of sqrt metric determinant? + get<0>(*result) = -get<1>(x) * orbital_angular_velocity; + get<1>(*result) = get<0>(x) * orbital_angular_velocity; get<2>(*result) = 0.0; }