diff --git a/src/Evolution/Systems/CurvedScalarWave/Actions/CMakeLists.txt b/src/Evolution/Systems/CurvedScalarWave/Actions/CMakeLists.txt new file mode 100644 index 0000000000000..134756a4a8ca0 --- /dev/null +++ b/src/Evolution/Systems/CurvedScalarWave/Actions/CMakeLists.txt @@ -0,0 +1,15 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +spectre_target_sources( + ${LIBRARY} + PRIVATE + NumericInitialData.cpp + ) + +spectre_target_headers( + ${LIBRARY} + INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src + HEADERS + NumericInitialData.hpp + ) diff --git a/src/Evolution/Systems/CurvedScalarWave/Actions/NumericInitialData.cpp b/src/Evolution/Systems/CurvedScalarWave/Actions/NumericInitialData.cpp new file mode 100644 index 0000000000000..3e24fef4e8e66 --- /dev/null +++ b/src/Evolution/Systems/CurvedScalarWave/Actions/NumericInitialData.cpp @@ -0,0 +1,48 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Systems/CurvedScalarWave/Actions/NumericInitialData.hpp" + +#include +#include +#include +#include + +#include "Utilities/PrettyType.hpp" + +namespace CurvedScalarWave { + +NumericInitialData::NumericInitialData( + std::string file_glob, std::string subfile_name, + std::variant observation_value, + const bool enable_interpolation, ScalarVars selected_variables) + : importer_options_(std::move(file_glob), std::move(subfile_name), + observation_value, enable_interpolation), + selected_variables_(std::move(selected_variables)) {} + +NumericInitialData::NumericInitialData(CkMigrateMessage* msg) + : InitialData(msg) {} + +PUP::able::PUP_ID NumericInitialData::my_PUP_ID = 0; + +size_t NumericInitialData::volume_data_id() const { + size_t hash = 0; + boost::hash_combine(hash, pretty_type::get_name()); + boost::hash_combine(hash, + get(importer_options_)); + boost::hash_combine(hash, + get(importer_options_)); + return hash; +} + +void NumericInitialData::pup(PUP::er& p) { + p | importer_options_; + p | selected_variables_; +} + +bool operator==(const NumericInitialData& lhs, const NumericInitialData& rhs) { + return lhs.importer_options_ == rhs.importer_options_ and + lhs.selected_variables_ == rhs.selected_variables_; +} + +} // namespace CurvedScalarWave diff --git a/src/Evolution/Systems/CurvedScalarWave/Actions/NumericInitialData.hpp b/src/Evolution/Systems/CurvedScalarWave/Actions/NumericInitialData.hpp new file mode 100644 index 0000000000000..42a6a736c2dbf --- /dev/null +++ b/src/Evolution/Systems/CurvedScalarWave/Actions/NumericInitialData.hpp @@ -0,0 +1,247 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include + +#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/Tensor/EagerMath/DotProduct.hpp" +#include "DataStructures/Tensor/EagerMath/RaiseOrLowerIndex.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/Tags.hpp" +#include "Evolution/Initialization/InitialData.hpp" +#include "Evolution/Systems/CurvedScalarWave/Tags.hpp" +#include "IO/Importers/Actions/ReadVolumeData.hpp" +#include "IO/Importers/ElementDataReader.hpp" +#include "IO/Importers/Tags.hpp" +#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "Options/String.hpp" +#include "Parallel/AlgorithmExecution.hpp" +#include "Parallel/GlobalCache.hpp" +#include "Parallel/Invoke.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp" +#include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp" +#include "Utilities/CallWithDynamicType.hpp" +#include "Utilities/ErrorHandling/Error.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/Serialization/CharmPupable.hpp" +#include "Utilities/SetNumberOfGridPoints.hpp" +#include "Utilities/TMPL.hpp" +#include "Utilities/TaggedTuple.hpp" + +namespace CurvedScalarWave { + +/*! + * \brief Numeric initial data loaded from volume data files + */ +class NumericInitialData : public evolution::initial_data::InitialData { + public: + template + struct VarName { + using tag = Tag; + static std::string name() { return db::tag_name(); } + using type = std::string; + static constexpr Options::String help = + "Name of the variable in the volume data file"; + }; + + // These are the scalar variables that we support loading from volume + // data files + using all_vars = + tmpl::list; + using optional_primitive_vars = tmpl::list<>; + + struct ScalarVars : tuples::tagged_tuple_from_typelist< + db::wrap_tags_in> { + static constexpr Options::String help = + "Scalar variables: 'Psi', 'Pi' and 'Phi'."; + using options = tags_list; + using TaggedTuple::TaggedTuple; + }; + + // Input-file options + struct Variables { + using type = ScalarVars; + static constexpr Options::String help = + "Set of initial data variables from which the Valencia evolution " + "variables are computed."; + }; + + using options = + tmpl::list; + + static constexpr Options::String help = + "Numeric initial data loaded from volume data files"; + + NumericInitialData() = default; + NumericInitialData(const NumericInitialData& rhs) = default; + NumericInitialData& operator=(const NumericInitialData& rhs) = default; + NumericInitialData(NumericInitialData&& /*rhs*/) = default; + NumericInitialData& operator=(NumericInitialData&& /*rhs*/) = default; + ~NumericInitialData() = default; + + /// \cond + explicit NumericInitialData(CkMigrateMessage* msg); + using PUP::able::register_constructor; + WRAPPED_PUPable_decl_template(NumericInitialData); + /// \endcond + + std::unique_ptr get_clone() + const override { + return std::make_unique(*this); + } + + NumericInitialData( + std::string file_glob, std::string subfile_name, + std::variant observation_value, + bool enable_interpolation, ScalarVars selected_variables); + + const importers::ImporterOptions& importer_options() const { + return importer_options_; + } + + const ScalarVars& selected_variables() const { return selected_variables_; } + + size_t volume_data_id() const; + + template + void select_for_import( + const gsl::not_null*> fields) const { + // Select the subset of the available variables that we want to read from + // the volume data file + using selected_vars = std::decay_t; + tmpl::for_each( + [&fields, this](const auto tag_v) { + using tag = typename std::decay_t::type::tag; + get>(*fields) = + get>(selected_variables_); + }); + } + + template + void set_initial_data( + const gsl::not_null*> psi_scalar, + const gsl::not_null*> pi_scalar, + const gsl::not_null*> phi_scalar, + const gsl::not_null*> numeric_data, + const Mesh<3>& mesh, + const InverseJacobian& inv_jacobian) const { + *psi_scalar = std::move(get(*numeric_data)); + *pi_scalar = std::move(get(*numeric_data)); + // Set Phi to the numerical spatial derivative of the scalar + partial_derivative(phi_scalar, *psi_scalar, mesh, inv_jacobian); + } + + void pup(PUP::er& p) override; + + friend bool operator==(const NumericInitialData& lhs, + const NumericInitialData& rhs); + + private: + importers::ImporterOptions importer_options_{}; + ScalarVars selected_variables_{}; +}; + +namespace Actions { + +/*! + * \brief Dispatch loading numeric initial data from files. + * + * Place this action before + * CurvedScalarWave::Actions::SetNumericInitialData in the action list. + * See importers::Actions::ReadAllVolumeDataAndDistribute for details, which is + * invoked by this action. + */ +struct ReadNumericInitialData { + using const_global_cache_tags = + tmpl::list; + + template + static Parallel::iterable_action_return_t apply( + db::DataBox& box, + const tuples::TaggedTuple& /*inboxes*/, + Parallel::GlobalCache& cache, + const ArrayIndex& /*array_index*/, const ActionList /*meta*/, + const ParallelComponent* const /*meta*/) { + // Select the subset of the available variables that we want to read from + // the volume data file + const auto& initial_data = dynamic_cast( + db::get(box)); + tuples::tagged_tuple_from_typelist> + selected_fields{}; + initial_data.select_for_import(make_not_null(&selected_fields)); + auto& reader_component = Parallel::get_parallel_component< + importers::ElementDataReader>(cache); + Parallel::simple_action>( + reader_component, initial_data.importer_options(), + initial_data.volume_data_id(), std::move(selected_fields)); + return {Parallel::AlgorithmExecution::Continue, std::nullopt}; + } +}; + +/*! + * \brief Receive numeric initial data loaded by + * CurvedScalarWave::Actions::ReadNumericInitialData. + */ +struct SetNumericInitialData { + static constexpr size_t Dim = 3; + using inbox_tags = + tmpl::list>; + + template + static Parallel::iterable_action_return_t apply( + db::DataBox& box, tuples::TaggedTuple& inboxes, + const Parallel::GlobalCache& /*cache*/, + const ElementId& /*element_id*/, const ActionList /*meta*/, + const ParallelComponent* const /*meta*/) { + auto& inbox = + tuples::get>( + inboxes); + const auto& initial_data = dynamic_cast( + db::get(box)); + const size_t volume_data_id = initial_data.volume_data_id(); + if (inbox.find(volume_data_id) == inbox.end()) { + return {Parallel::AlgorithmExecution::Retry, std::nullopt}; + } + auto numeric_data = std::move(inbox.extract(volume_data_id).mapped()); + + const auto& mesh = db::get>(box); + const auto& inv_jacobian = + db::get>(box); + + db::mutate>( + [&initial_data, &numeric_data, &mesh, &inv_jacobian]( + const gsl::not_null*> psi_scalar, + const gsl::not_null*> pi_scalar, + const gsl::not_null*> phi_scalar) { + initial_data.set_initial_data(psi_scalar, pi_scalar, phi_scalar, + make_not_null(&numeric_data), mesh, + inv_jacobian); + }, + make_not_null(&box)); + + return {Parallel::AlgorithmExecution::Continue, std::nullopt}; + } +}; + +} // namespace Actions + +} // namespace CurvedScalarWave diff --git a/src/Evolution/Systems/CurvedScalarWave/CMakeLists.txt b/src/Evolution/Systems/CurvedScalarWave/CMakeLists.txt index f29fe6e63e378..820077d67b5b7 100644 --- a/src/Evolution/Systems/CurvedScalarWave/CMakeLists.txt +++ b/src/Evolution/Systems/CurvedScalarWave/CMakeLists.txt @@ -42,6 +42,7 @@ target_link_libraries( GeneralRelativity ) +add_subdirectory(Actions) add_subdirectory(BoundaryConditions) add_subdirectory(BoundaryCorrections) add_subdirectory(Worldtube) diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Actions/Test_NumericInitialData.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Actions/Test_NumericInitialData.cpp new file mode 100644 index 0000000000000..0e95720584917 --- /dev/null +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Actions/Test_NumericInitialData.cpp @@ -0,0 +1,220 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "DataStructures/VariablesTag.hpp" +#include "Domain/CoordinateMaps/CoordinateMap.tpp" +#include "Domain/CoordinateMaps/Wedge.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/Tags.hpp" +#include "Evolution/Systems/CurvedScalarWave/Actions/NumericInitialData.hpp" +#include "Framework/ActionTesting.hpp" +#include "Framework/TestCreation.hpp" +#include "IO/Importers/Actions/ReadVolumeData.hpp" +#include "IO/Importers/ElementDataReader.hpp" +#include "IO/Importers/Tags.hpp" +#include "NumericalAlgorithms/Spectral/Basis.hpp" +#include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "NumericalAlgorithms/Spectral/Quadrature.hpp" +#include "Options/Protocols/FactoryCreation.hpp" +#include "Parallel/Phase.hpp" +#include "PointwiseFunctions/AnalyticData/CurvedWaveEquation/PureSphericalHarmonic.hpp" +#include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp" +#include "Utilities/TaggedTuple.hpp" + +namespace CurvedScalarWave::Actions { + +namespace { + +using all_scalar_vars = + tmpl::list>; + +template +struct MockElementArray { + using component_being_mocked = void; // Not needed + using metavariables = Metavariables; + using chare_type = ActionTesting::MockArrayChare; + using array_index = ElementId<3>; + using phase_dependent_action_list = tmpl::list< + Parallel::PhaseActions< + Parallel::Phase::Initialization, + tmpl::list, domain::Tags::Mesh<3>, + domain::Tags::Coordinates<3, Frame::Inertial>, + domain::Tags::InverseJacobian<3, Frame::ElementLogical, + Frame::Inertial>>>>>, + Parallel::PhaseActions< + Parallel::Phase::Testing, + tmpl::list>>; +}; + +struct MockReadVolumeData { + template + static void apply( + DataBox& /*box*/, Parallel::GlobalCache& cache, + const ArrayIndex& /*array_index*/, + const importers::ImporterOptions& options, const size_t volume_data_id, + tuples::tagged_tuple_from_typelist> + selected_fields) { + const auto& initial_data = dynamic_cast( + get(cache)); + CHECK(options == initial_data.importer_options()); + CHECK(volume_data_id == initial_data.volume_data_id()); + CHECK(get>( + selected_fields) == "CustomPsi"); + CHECK(get>( + selected_fields) == "CustomPi"); + } +}; + +template +struct MockVolumeDataReader { + using component_being_mocked = importers::ElementDataReader; + using metavariables = Metavariables; + using chare_type = ActionTesting::MockNodeGroupChare; + using array_index = size_t; + using phase_dependent_action_list = tmpl::list< + Parallel::PhaseActions>>; + using replace_these_simple_actions = + tmpl::list>>; + using with_these_simple_actions = tmpl::list; +}; + +struct Metavariables { + static constexpr size_t volume_dim = 3; + using component_list = tmpl::list, + MockVolumeDataReader>; + + struct factory_creation + : tt::ConformsTo { + using factory_classes = + tmpl::map>>; + }; +}; + +void test_numeric_initial_data(const NumericInitialData& initial_data, + const std::string& option_string) { + { + INFO("Factory creation"); + const auto created = TestHelpers::test_creation< + std::unique_ptr, Metavariables>( + option_string); + CHECK(dynamic_cast(*created) == initial_data); + } + + using reader_component = MockVolumeDataReader; + using element_array = MockElementArray; + + ActionTesting::MockRuntimeSystem runner{ + initial_data.get_clone()}; + + const double wave_radius = 2.0; + + CurvedScalarWave::AnalyticData::PureSphericalHarmonic csw_sh_wave{ + wave_radius, 1.0, std::pair{0, 0}}; + + // Setup mock data file reader + ActionTesting::emplace_nodegroup_component( + make_not_null(&runner)); + + // Setup element + const ElementId<3> element_id{0, {{{1, 0}, {1, 0}, {1, 0}}}}; + const Mesh<3> mesh{6, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + const size_t num_points = mesh.number_of_grid_points(); + const auto map = + domain::make_coordinate_map( + domain::CoordinateMaps::Wedge<3>{ + wave_radius / 2., wave_radius * 2., 1., 1., {}, true}); + const auto logical_coords = logical_coordinates(mesh); + const auto coords = map(logical_coords); + const auto inv_jacobian = map.inv_jacobian(logical_coords); + ActionTesting::emplace_component_and_initialize( + make_not_null(&runner), element_id, + {Variables{num_points}, mesh, coords, inv_jacobian}); + + const auto get_element_tag = [&runner, + &element_id](auto tag_v) -> decltype(auto) { + using tag = std::decay_t; + return ActionTesting::get_databox_tag(runner, + element_id); + }; + + ActionTesting::set_phase(make_not_null(&runner), Parallel::Phase::Testing); + + // ReadNumericInitialData + ActionTesting::next_action(make_not_null(&runner), element_id); + REQUIRE_FALSE(ActionTesting::next_action_if_ready( + make_not_null(&runner), element_id)); + + // MockReadVolumeData + ActionTesting::invoke_queued_simple_action( + make_not_null(&runner), 0); + + // Insert Csw data into the inbox + auto& inbox = ActionTesting::get_inbox_tag< + element_array, importers::Tags::VolumeData, + Metavariables>(make_not_null(&runner), + element_id)[initial_data.volume_data_id()]; + const auto& selected_vars = initial_data.selected_variables(); + auto csw_sh_wave_vars = csw_sh_wave.variables(coords, all_scalar_vars{}); + get(inbox) = + get(csw_sh_wave_vars); + get(inbox) = + get(csw_sh_wave_vars); + + // SetNumericInitialData + ActionTesting::next_action(make_not_null(&runner), element_id); + + // Check result + tmpl::for_each([&get_element_tag, + &csw_sh_wave_vars](const auto tag_v) { + using tag = tmpl::type_from>; + CAPTURE(pretty_type::name()); + CHECK_ITERABLE_APPROX(get_element_tag(tag{}), (get(csw_sh_wave_vars))); + }); +} + +} // namespace + +SPECTRE_TEST_CASE("Unit.Evolution.Systems.CurvedScalarWave.SetInitialData", + "[Unit][Evolution]") { + register_factory_classes_with_charm(); + test_numeric_initial_data(NumericInitialData{"TestInitialData.h5", + "VolumeData", + 0., + false, + {"CustomPsi", "CustomPi"}}, + "NumericInitialData:\n" + " FileGlob: TestInitialData.h5\n" + " Subgroup: VolumeData\n" + " ObservationValue: 0.\n" + " Interpolate: False\n" + " Variables:\n" + " Psi: CustomPsi\n" + " Pi: CustomPi"); +} + +} // namespace CurvedScalarWave::Actions diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt b/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt index f7e16ff0583b0..76f5522e2f326 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/CMakeLists.txt @@ -6,6 +6,7 @@ add_subdirectory(Worldtube) set(LIBRARY "Test_CurvedScalarWave") set(LIBRARY_SOURCES + Actions/Test_NumericInitialData.cpp BoundaryCorrections/Test_UpwindPenalty.cpp BoundaryConditions/Test_AnalyticConstant.cpp BoundaryConditions/Test_ConstraintPreservingSphericalRadiation.cpp @@ -34,6 +35,7 @@ target_link_libraries( DomainBoundaryConditions DomainBoundaryConditionsHelpers GeneralRelativityHelpers + Importers MathFunctions Time Utilities