From cba96a9845061be751263a78491fe4abe8c6949a Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Fri, 16 Feb 2024 17:12:00 -0800 Subject: [PATCH 1/9] Add AMR logging tag to component Before it wasn't reliably added. --- src/ParallelAlgorithms/Amr/Actions/Component.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ParallelAlgorithms/Amr/Actions/Component.hpp b/src/ParallelAlgorithms/Amr/Actions/Component.hpp index a267a46df030..cc770f9317d0 100644 --- a/src/ParallelAlgorithms/Amr/Actions/Component.hpp +++ b/src/ParallelAlgorithms/Amr/Actions/Component.hpp @@ -31,7 +31,8 @@ struct Component { using chare_type = Parallel::Algorithms::Singleton; using const_global_cache_tags = - tmpl::list; + tmpl::list>; using phase_dependent_action_list = tmpl::list< Parallel::PhaseActions>>; From b093a03b6cc767351e055722b12348d930aadba9 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Fri, 26 Jan 2024 13:05:16 -0800 Subject: [PATCH 2/9] Ignore a clang-tidy check --- .clang-tidy | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.clang-tidy b/.clang-tidy index abd22e691ec5..e467891b616d 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -78,6 +78,8 @@ Checks: '*, -readability-magic-numbers, # Same as llvm-qualified-auto above., -readability-qualified-auto, +# Access specifiers can be useful to structure code, +-readability-redundant-access-specifiers, # We can have two of: methods are static when possible, static, # methods are not called through instances, and methods of, # calling, e.g., x.size(), are consistent across classes. We, From 6414f7325b438213b90421a8189061acedbfde97 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Thu, 25 Jan 2024 17:09:00 -0800 Subject: [PATCH 3/9] Add AMR projector for elliptic sources --- .../Actions/InitializeFixedSources.hpp | 60 ++++++++++++------- src/Elliptic/CMakeLists.txt | 1 + .../Actions/Test_InitializeFixedSources.cpp | 26 +++++--- 3 files changed, 57 insertions(+), 30 deletions(-) diff --git a/src/Elliptic/Actions/InitializeFixedSources.hpp b/src/Elliptic/Actions/InitializeFixedSources.hpp index 0dda06aa5cd8..dd2227d669f2 100644 --- a/src/Elliptic/Actions/InitializeFixedSources.hpp +++ b/src/Elliptic/Actions/InitializeFixedSources.hpp @@ -14,10 +14,12 @@ #include "Domain/Structure/ElementId.hpp" #include "Domain/Tags.hpp" #include "Elliptic/DiscontinuousGalerkin/Tags.hpp" -#include "Elliptic/Utilities/GetAnalyticData.hpp" #include "NumericalAlgorithms/DiscontinuousGalerkin/ApplyMassMatrix.hpp" #include "Parallel/AlgorithmExecution.hpp" +#include "Parallel/Tags/Metavariables.hpp" +#include "ParallelAlgorithms/Amr/Protocols/Projector.hpp" #include "ParallelAlgorithms/Initialization/MutateAssign.hpp" +#include "Utilities/CallWithDynamicType.hpp" #include "Utilities/MakeWithValue.hpp" #include "Utilities/TMPL.hpp" #include "Utilities/TaggedTuple.hpp" @@ -49,50 +51,64 @@ namespace elliptic::Actions { * - Adds: * - `db::wrap_tags_in<::Tags::FixedSource, primal_fields>` */ -template -struct InitializeFixedSources { +template +struct InitializeFixedSources : tt::ConformsTo<::amr::protocols::Projector> { private: using fixed_sources_tag = ::Tags::Variables< db::wrap_tags_in<::Tags::FixedSource, typename System::primal_fields>>; - public: + public: // Iterable action using const_global_cache_tags = tmpl::list; using simple_tags = tmpl::list; using compute_tags = tmpl::list<>; template + typename ActionList, typename ParallelComponent> static Parallel::iterable_action_return_t apply( db::DataBox& box, const tuples::TaggedTuple& /*inboxes*/, const Parallel::GlobalCache& /*cache*/, const ElementId& /*array_index*/, const ActionList /*meta*/, const ParallelComponent* const /*meta*/) { - const auto& inertial_coords = - get>(box); - const auto& background = db::get(box); + db::mutate_apply(make_not_null(&box)); + return {Parallel::AlgorithmExecution::Continue, std::nullopt}; + } + + public: // DataBox mutator, amr::protocols::Projector + using return_tags = tmpl::list; + using argument_tags = tmpl::list< + domain::Tags::Coordinates, BackgroundTag, + elliptic::dg::Tags::Massive, domain::Tags::Mesh, + domain::Tags::DetInvJacobian, + Parallel::Tags::Metavariables>; + template + static void apply( + 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*/) { // Retrieve the fixed-sources of the elliptic system from the background, // which (along with the boundary conditions) define the problem we want to // solve. - auto fixed_sources = elliptic::util::get_analytic_data< - typename fixed_sources_tag::tags_list>(background, box, - inertial_coords); + using factory_classes = + typename std::decay_t::factory_creation::factory_classes; + *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{})); + }); // Apply DG mass matrix to the fixed sources if the DG operator is massive - if (db::get(box)) { - const auto& mesh = db::get>(box); - const auto& det_inv_jacobian = db::get< - domain::Tags::DetInvJacobian>( - box); - fixed_sources /= get(det_inv_jacobian); - ::dg::apply_mass_matrix(make_not_null(&fixed_sources), mesh); + if (massive) { + *fixed_sources /= get(det_inv_jacobian); + ::dg::apply_mass_matrix(fixed_sources, mesh); } - - ::Initialization::mutate_assign(make_not_null(&box), - std::move(fixed_sources)); - return {Parallel::AlgorithmExecution::Continue, std::nullopt}; } }; diff --git a/src/Elliptic/CMakeLists.txt b/src/Elliptic/CMakeLists.txt index 6e27b89cc6c3..a7295fbaba14 100644 --- a/src/Elliptic/CMakeLists.txt +++ b/src/Elliptic/CMakeLists.txt @@ -23,6 +23,7 @@ target_link_libraries( Domain EventsAndTriggers LinearOperators + ParallelAmr Serialization Utilities ) diff --git a/tests/Unit/Elliptic/Actions/Test_InitializeFixedSources.cpp b/tests/Unit/Elliptic/Actions/Test_InitializeFixedSources.cpp index 2786d7e31ff5..417c07724874 100644 --- a/tests/Unit/Elliptic/Actions/Test_InitializeFixedSources.cpp +++ b/tests/Unit/Elliptic/Actions/Test_InitializeFixedSources.cpp @@ -28,6 +28,7 @@ #include "Framework/ActionTesting.hpp" #include "Options/Protocols/FactoryCreation.hpp" #include "Parallel/Phase.hpp" +#include "ParallelAlgorithms/Amr/Protocols/Projector.hpp" #include "PointwiseFunctions/InitialDataUtilities/Background.hpp" #include "Utilities/ProtocolHelpers.hpp" #include "Utilities/Serialization/CharmPupable.hpp" @@ -42,6 +43,7 @@ struct ScalarFieldTag : db::SimpleTag { }; struct System { + static constexpr size_t volume_dim = 1; using primal_fields = tmpl::list; }; @@ -105,6 +107,13 @@ struct Metavariables { SPECTRE_TEST_CASE("Unit.Elliptic.Actions.InitializeFixedSources", "[Unit][Elliptic][Actions]") { + using background_tag = + elliptic::Tags::Background; + static_assert( + tt::assert_conforms_to_v< + elliptic::Actions::InitializeFixedSources, + amr::protocols::Projector>); + domain::creators::register_derived_with_charm(); register_factory_classes_with_charm(); // Which element we work with does not matter for this test @@ -113,13 +122,14 @@ SPECTRE_TEST_CASE("Unit.Elliptic.Actions.InitializeFixedSources", {{-0.5}}, {{1.5}}, {{2}}, {{4}}}; using element_array = ElementArray; - ActionTesting::MockRuntimeSystem runner{tuples::TaggedTuple< - elliptic::Tags::Background, - domain::Tags::Domain<1>, domain::Tags::FunctionsOfTimeInitialize, - elliptic::dg::Tags::Massive, elliptic::dg::Tags::Quadrature>{ - std::make_unique(), domain_creator.create_domain(), - domain_creator.functions_of_time(), false, - Spectral::Quadrature::GaussLobatto}}; + ActionTesting::MockRuntimeSystem runner{ + tuples::TaggedTuple, + domain::Tags::FunctionsOfTimeInitialize, + elliptic::dg::Tags::Massive, + elliptic::dg::Tags::Quadrature>{ + std::make_unique(), domain_creator.create_domain(), + domain_creator.functions_of_time(), false, + Spectral::Quadrature::GaussLobatto}}; ActionTesting::emplace_component_and_initialize( &runner, element_id, {domain_creator.initial_refinement_levels(), @@ -127,7 +137,7 @@ SPECTRE_TEST_CASE("Unit.Elliptic.Actions.InitializeFixedSources", ActionTesting::next_action(make_not_null(&runner), element_id); ActionTesting::set_phase(make_not_null(&runner), Parallel::Phase::Testing); ActionTesting::next_action(make_not_null(&runner), element_id); - const auto get_tag = [&runner, &element_id ](auto tag_v) -> const auto& { + const auto get_tag = [&runner, &element_id](auto tag_v) -> const auto& { using tag = std::decay_t; return ActionTesting::get_databox_tag(runner, element_id); From 794eb414d40298e8a78aa315824c006fdd24fd84 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Fri, 26 Jan 2024 10:29:12 -0800 Subject: [PATCH 4/9] Add AMR projector for elliptic analytic solutions --- .../Actions/InitializeAnalyticSolution.hpp | 64 ++++++++++++++----- src/Elliptic/Executables/Solver.hpp | 2 +- .../Test_InitializeAnalyticSolution.cpp | 15 ++++- .../DiscontinuousGalerkin/Test_DgOperator.cpp | 2 +- 4 files changed, 64 insertions(+), 19 deletions(-) diff --git a/src/Elliptic/Actions/InitializeAnalyticSolution.hpp b/src/Elliptic/Actions/InitializeAnalyticSolution.hpp index eae32d1d5f55..9c5abcbf38a4 100644 --- a/src/Elliptic/Actions/InitializeAnalyticSolution.hpp +++ b/src/Elliptic/Actions/InitializeAnalyticSolution.hpp @@ -12,14 +12,16 @@ #include "DataStructures/DataBox/PrefixHelpers.hpp" #include "DataStructures/DataBox/Prefixes.hpp" #include "DataStructures/Tensor/Tensor.hpp" -#include "DataStructures/VariablesTag.hpp" +#include "DataStructures/Variables.hpp" #include "Domain/Structure/ElementId.hpp" #include "Domain/Tags.hpp" -#include "Elliptic/Utilities/GetAnalyticData.hpp" #include "Parallel/AlgorithmExecution.hpp" #include "Parallel/GlobalCache.hpp" +#include "Parallel/Tags/Metavariables.hpp" +#include "ParallelAlgorithms/Amr/Protocols/Projector.hpp" #include "ParallelAlgorithms/Initialization/MutateAssign.hpp" #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp" +#include "Utilities/CallWithDynamicType.hpp" #include "Utilities/TMPL.hpp" #include "Utilities/TaggedTuple.hpp" @@ -49,37 +51,69 @@ namespace elliptic::Actions { * - Adds: * - `::Tags::AnalyticSolutionsBase` */ -template -struct InitializeOptionalAnalyticSolution { +struct InitializeOptionalAnalyticSolution + : tt::ConformsTo<::amr::protocols::Projector> { private: using analytic_fields_tag = ::Tags::AnalyticSolutions; - public: + public: // Iterable action using simple_tags = tmpl::list; using compute_tags = tmpl::list<>; using const_global_cache_tags = tmpl::list; template + typename ActionList, typename ParallelComponent> static Parallel::iterable_action_return_t apply( db::DataBox& box, const tuples::TaggedTuple& /*inboxes*/, const Parallel::GlobalCache& /*cache*/, const ElementId& /*array_index*/, const ActionList /*meta*/, const ParallelComponent* const /*meta*/) { + db::mutate_apply(make_not_null(&box)); + return {Parallel::AlgorithmExecution::Continue, std::nullopt}; + } + + public: // DataBox mutator, amr::protocols::Projector + using return_tags = tmpl::list; + using argument_tags = + tmpl::list, + domain::Tags::Coordinates, BackgroundTag, + Parallel::Tags::Metavariables>; + + template + static void apply(const gsl::not_null + analytic_solution_fields, + const Mesh& mesh, + const tnsr::I inertial_coords, + const Background& background, const Metavariables& /*meta*/, + const AmrData&... amr_data) { + if constexpr (sizeof...(AmrData) == 1) { + if constexpr (std::is_same_v, Element>>) { + if (((mesh == amr_data.first) and ...)) { + // This element hasn't changed during AMR. Nothing to do. + return; + } + } + } + const auto analytic_solution = - dynamic_cast(&db::get(box)); + dynamic_cast(&background); if (analytic_solution != nullptr) { - const auto& inertial_coords = - get>(box); - auto analytic_fields = - elliptic::util::get_analytic_data( - *analytic_solution, box, inertial_coords); - Initialization::mutate_assign(make_not_null(&box), - std::move(analytic_fields)); + using factory_classes = typename std::decay_t< + Metavariables>::factory_creation::factory_classes; + *analytic_solution_fields = call_with_dynamic_type< + Variables, + tmpl::at>( + analytic_solution, [&inertial_coords](const auto* const derived) { + return variables_from_tagged_tuple( + derived->variables(inertial_coords, AnalyticSolutionFields{})); + }); + } else { + *analytic_solution_fields = std::nullopt; } - return {Parallel::AlgorithmExecution::Continue, std::nullopt}; } }; /// @} diff --git a/src/Elliptic/Executables/Solver.hpp b/src/Elliptic/Executables/Solver.hpp index 47866e547697..985385de867a 100644 --- a/src/Elliptic/Executables/Solver.hpp +++ b/src/Elliptic/Executables/Solver.hpp @@ -208,7 +208,7 @@ struct Solver { ::Actions::RandomizeVariables, elliptic::Actions::InitializeFixedSources, elliptic::Actions::InitializeOptionalAnalyticSolution< - background_tag, + volume_dim, background_tag, tmpl::append, elliptic::analytic_data::AnalyticSolution>, diff --git a/tests/Unit/Elliptic/Actions/Test_InitializeAnalyticSolution.cpp b/tests/Unit/Elliptic/Actions/Test_InitializeAnalyticSolution.cpp index d15624600e8b..1f8c36f3c17c 100644 --- a/tests/Unit/Elliptic/Actions/Test_InitializeAnalyticSolution.cpp +++ b/tests/Unit/Elliptic/Actions/Test_InitializeAnalyticSolution.cpp @@ -19,6 +19,7 @@ #include "Options/Protocols/FactoryCreation.hpp" #include "Parallel/ParallelComponentHelpers.hpp" #include "Parallel/Phase.hpp" +#include "ParallelAlgorithms/Amr/Protocols/Projector.hpp" #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp" #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp" #include "PointwiseFunctions/InitialDataUtilities/Background.hpp" @@ -76,11 +77,13 @@ struct ElementArray { Parallel::PhaseActions< Parallel::Phase::Initialization, tmpl::list>>>>, + tmpl::list, + domain::Tags::Coordinates<1, Frame::Inertial>>>>>, Parallel::PhaseActions< Parallel::Phase::Testing, tmpl::list, tmpl::list, elliptic::analytic_data::AnalyticSolution>>>>; @@ -115,7 +118,7 @@ void test_initialize_analytic_solution( {std::move(analytic_solution_or_data)}}; const ElementId<1> element_id{0}; ActionTesting::emplace_component_and_initialize( - &runner, element_id, {inertial_coords}); + &runner, element_id, {Mesh<1>{}, inertial_coords}); ActionTesting::set_phase(make_not_null(&runner), Parallel::Phase::Testing); for (size_t i = 0; i < 2; ++i) { @@ -148,6 +151,14 @@ void test_initialize_analytic_solution( SPECTRE_TEST_CASE("Unit.Elliptic.Actions.InitializeAnalyticSolution", "[Unit][Elliptic][Actions]") { + static_assert( + tt::assert_conforms_to_v< + elliptic::Actions::InitializeOptionalAnalyticSolution< + 1, + elliptic::Tags::Background, + tmpl::list, + elliptic::analytic_data::AnalyticSolution>, + amr::protocols::Projector>); test_initialize_analytic_solution( tnsr::I{{{{1., 2., 3., 4.}}}}, Scalar{{{{2., 4., 6., 8.}}}}); diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp index 7684170faf5b..5e6cd94243de 100644 --- a/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp +++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp @@ -131,7 +131,7 @@ struct ElementArray { domain::Tags::InitialExtents>>, ::elliptic::dg::Actions::InitializeDomain, ::elliptic::Actions::InitializeOptionalAnalyticSolution< - analytic_solution_tag, + Dim, analytic_solution_tag, tmpl::append, typename metavariables::analytic_solution>, From b9dfb5ff3c2a2a1b38e9bf378b5c4ed102a735eb Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Fri, 26 Jan 2024 13:41:47 -0800 Subject: [PATCH 5/9] Add AMR projector for elliptic domain --- .../DiscontinuousGalerkin/Initialization.cpp | 61 ++++++++---- .../DiscontinuousGalerkin/Initialization.hpp | 96 +++++++++++++++++++ 2 files changed, 139 insertions(+), 18 deletions(-) diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp index 4e19e8892c36..6e5289fed5e2 100644 --- a/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp +++ b/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp @@ -36,6 +36,44 @@ namespace elliptic::dg { +namespace detail{ +template +void initialize_coords_and_jacobians( + gsl::not_null*> + logical_coords, + gsl::not_null*> inertial_coords, + gsl::not_null*> + inv_jacobian, + gsl::not_null*> det_inv_jacobian, + gsl::not_null*> det_jacobian, + gsl::not_null*> + det_times_inv_jacobian, + const Mesh& mesh, + const ElementMap& element_map, + const domain::FunctionsOfTimeMap& functions_of_time) { + // Coordinates + *logical_coords = logical_coordinates(mesh); + *inertial_coords = + element_map(*logical_coords, 0., functions_of_time); + // Jacobian + // Note: we can try to use `::dg::metric_identity_jacobian_quantities` here. + // When I tried (NV, Dec 2023) the DG residuals diverged on a sphere domain + // with a large outer boundary (1e9). This was possibly because no + // metric-identity Jacobians were used on faces, though I also tried slicing + // the metric-identity Jacobian to the faces and that didn't help. + *inv_jacobian = + element_map.inv_jacobian(*logical_coords, 0., functions_of_time); + *det_inv_jacobian = determinant(*inv_jacobian); + get(*det_jacobian) = 1. / get(*det_inv_jacobian); + *det_times_inv_jacobian = *inv_jacobian; + for (auto& component : *det_times_inv_jacobian) { + component *= get(*det_jacobian); + } + } +} // namespace detail + template void InitializeGeometry::apply( const gsl::not_null*> mesh, @@ -81,24 +119,11 @@ void InitializeGeometry::apply( } // Element map *element_map = ElementMap{element_id, block}; - // Coordinates - *logical_coords = logical_coordinates(*mesh); - *inertial_coords = - element_map->operator()(*logical_coords, 0., functions_of_time); - // Jacobian - // Note: we can try to use `::dg::metric_identity_jacobian_quantities` here. - // When I tried (NV, Dec 2023) the DG residuals diverged on a sphere domain - // with a large outer boundary (1e9). This was possibly because no - // metric-identity Jacobians were used on faces, though I also tried slicing - // the metric-identity Jacobian to the faces and that didn't help. - *inv_jacobian = - element_map->inv_jacobian(*logical_coords, 0., functions_of_time); - *det_inv_jacobian = determinant(*inv_jacobian); - get(*det_jacobian) = 1. / get(*det_inv_jacobian); - *det_times_inv_jacobian = *inv_jacobian; - for (auto& component : *det_times_inv_jacobian) { - component *= get(*det_jacobian); - } + // Coordinates and Jacobians + detail::initialize_coords_and_jacobians( + logical_coords, inertial_coords, inv_jacobian, det_inv_jacobian, + det_jacobian, det_times_inv_jacobian, *mesh, *element_map, + functions_of_time); } namespace detail { diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp index 16ac1d70961f..df74fa700295 100644 --- a/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp +++ b/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp @@ -45,14 +45,34 @@ #include "NumericalAlgorithms/Spectral/Projection.hpp" #include "NumericalAlgorithms/Spectral/Quadrature.hpp" #include "Parallel/Tags/Metavariables.hpp" +#include "ParallelAlgorithms/Amr/Protocols/Projector.hpp" #include "Utilities/CallWithDynamicType.hpp" #include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/Gsl.hpp" +#include "Utilities/ProtocolHelpers.hpp" #include "Utilities/TMPL.hpp" #include "Utilities/TaggedTuple.hpp" namespace elliptic::dg { +namespace detail { +template +void initialize_coords_and_jacobians( + gsl::not_null*> + logical_coords, + gsl::not_null*> inertial_coords, + gsl::not_null*> + inv_jacobian, + gsl::not_null*> det_inv_jacobian, + gsl::not_null*> det_jacobian, + gsl::not_null*> + det_times_inv_jacobian, + const Mesh& mesh, const ElementMap& element_map, + const domain::FunctionsOfTimeMap& functions_of_time); +} // namespace detail + /*! * \brief Initialize the background-independent geometry for the elliptic DG * operator @@ -105,6 +125,82 @@ struct InitializeGeometry { Spectral::Quadrature quadrature, const ElementId& element_id); }; +template +struct ProjectGeometry : tt::ConformsTo<::amr::protocols::Projector> { + using return_tags = tmpl::list< + domain::Tags::ElementMap, + domain::Tags::Coordinates, + domain::Tags::Coordinates, + domain::Tags::InverseJacobian, + domain::Tags::DetInvJacobian, + domain::Tags::DetJacobian, + domain::Tags::DetTimesInvJacobian>; + using argument_tags = + tmpl::list, domain::Tags::Element, + domain::Tags::Domain, domain::Tags::FunctionsOfTime>; + + // p-refinement + static void apply( + const gsl::not_null*> element_map, + const gsl::not_null*> + logical_coords, + const gsl::not_null*> + inertial_coords, + const gsl::not_null*> + inv_jacobian, + const gsl::not_null*> det_inv_jacobian, + const gsl::not_null*> det_jacobian, + const gsl::not_null*> + det_times_inv_jacobian, + const Mesh& mesh, const Element& /*element*/, + const Domain& /*domain*/, + const domain::FunctionsOfTimeMap& functions_of_time, + const std::pair, Element>& old_mesh_and_element) { + if (mesh == old_mesh_and_element.first) { + // Neighbors may have changed, but this element hasn't. Nothing to do. + return; + } + // The element map doesn't change with p-refinement + detail::initialize_coords_and_jacobians( + logical_coords, inertial_coords, inv_jacobian, det_inv_jacobian, + det_jacobian, det_times_inv_jacobian, mesh, *element_map, + functions_of_time); + } + + // h-refinement + template + static void apply( + const gsl::not_null*> element_map, + const gsl::not_null*> + logical_coords, + const gsl::not_null*> + inertial_coords, + const gsl::not_null*> + inv_jacobian, + const gsl::not_null*> det_inv_jacobian, + const gsl::not_null*> det_jacobian, + const gsl::not_null*> + det_times_inv_jacobian, + const Mesh& mesh, const Element& element, + const Domain& domain, + const domain::FunctionsOfTimeMap& functions_of_time, + const ParentOrChildrenItemsType&... /*parent_or_children_items*/) { + const auto& element_id = element.id(); + const auto& block = domain.blocks()[element_id.block_id()]; + *element_map = ElementMap{element_id, block}; + detail::initialize_coords_and_jacobians( + logical_coords, inertial_coords, inv_jacobian, det_inv_jacobian, + det_jacobian, det_times_inv_jacobian, mesh, *element_map, + functions_of_time); + } +}; + namespace detail { // Compute derivative of the Jacobian numerically template From 2406070d30ad34bc1712df9fad957675d3121bd6 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Sun, 28 Jan 2024 16:20:02 -0800 Subject: [PATCH 6/9] Add AMR projector for elliptic faces & mortars --- .../DiscontinuousGalerkin/Initialization.hpp | 50 +++++++++++++++++-- 1 file changed, 45 insertions(+), 5 deletions(-) diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp index df74fa700295..d6c8ff9be1ba 100644 --- a/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp +++ b/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp @@ -234,7 +234,7 @@ tnsr::I mortar_logical_coordinates( /// The `::Tags::deriv>` is only added /// on external boundaries, for use by boundary conditions. template -struct InitializeFacesAndMortars { +struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> { using return_tags = tmpl::append< domain::make_faces_tags< Dim, @@ -275,8 +275,7 @@ struct InitializeFacesAndMortars { tmpl::conditional_t< std::is_same_v, tmpl::list<>, tmpl::list>>; - template + template static void apply( const gsl::not_null>*> face_directions, const gsl::not_null>*> @@ -306,8 +305,46 @@ struct InitializeFacesAndMortars { const InverseJacobian& inv_jacobian, const domain::FunctionsOfTimeMap& functions_of_time, - const Background& background = nullptr, - const Metavariables& /*meta*/ = nullptr) { + const AmrData&... amr_data) { + apply(face_directions, faces_inertial_coords, face_normals, + face_normal_vectors, face_normal_magnitudes, face_jacobians, + face_jacobian_times_inv_jacobian, deriv_unnormalized_face_normals, + mortar_meshes, mortar_sizes, mortar_jacobians, mesh, element, + neighbor_meshes, element_map, inv_jacobian, functions_of_time, + nullptr, nullptr, amr_data...); + } + template + static void apply( + const gsl::not_null>*> face_directions, + const gsl::not_null>*> + faces_inertial_coords, + const gsl::not_null>*> + face_normals, + const gsl::not_null>*> + face_normal_vectors, + const gsl::not_null>*> + face_normal_magnitudes, + const gsl::not_null>*> + face_jacobians, + const gsl::not_null>*> + face_jacobian_times_inv_jacobian, + const gsl::not_null>*> + deriv_unnormalized_face_normals, + const gsl::not_null<::dg::MortarMap>*> mortar_meshes, + const gsl::not_null<::dg::MortarMap>*> + mortar_sizes, + const gsl::not_null<::dg::MortarMap>*> + mortar_jacobians, + const Mesh& mesh, const Element& element, + const DirectionalIdMap>& neighbor_meshes, + const ElementMap& element_map, + const InverseJacobian& inv_jacobian, + const domain::FunctionsOfTimeMap& functions_of_time, + const Background& background, const Metavariables& /*meta*/, + const AmrData&... /*amr_data*/) { static_assert(std::is_same_v or not(std::is_same_v), "Supply an analytic background from which the 'InvMetricTag' " @@ -383,6 +420,9 @@ struct InitializeFacesAndMortars { detail::deriv_unnormalized_face_normals_impl( deriv_unnormalized_face_normals, mesh, element, inv_jacobian); // Mortars (internal directions) + mortar_meshes->clear(); + mortar_sizes->clear(); + mortar_jacobians->clear(); const auto& element_id = element.id(); for (const auto& [direction, neighbors] : element.neighbors()) { const auto face_mesh = mesh.slice_away(direction.dimension()); From 7b921b697bfa85a38e3c2f5d49192bb8462b94d4 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Sun, 28 Jan 2024 16:20:38 -0800 Subject: [PATCH 7/9] Add AMR projector for elliptic background --- .../DiscontinuousGalerkin/Initialization.hpp | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp index d6c8ff9be1ba..fd602faffd97 100644 --- a/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp +++ b/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp @@ -488,7 +488,7 @@ struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> { /// Initialize background quantities for the elliptic DG operator, possibly /// including the metric necessary for normalizing face normals template -struct InitializeBackground { +struct InitializeBackground : tt::ConformsTo<::amr::protocols::Projector> { using return_tags = tmpl::list<::Tags::Variables, domain::Tags::Faces>>; @@ -499,7 +499,8 @@ struct InitializeBackground { Frame::Inertial>, BackgroundTag, Parallel::Tags::Metavariables>; - template + template static void apply( const gsl::not_null*> background_fields, const gsl::not_null>*> @@ -507,7 +508,18 @@ struct InitializeBackground { const tnsr::I& inertial_coords, const Mesh& mesh, const InverseJacobian& inv_jacobian, - const BackgroundBase& background, const Metavariables& /*meta*/) { + const BackgroundBase& background, const Metavariables& /*meta*/, + const AmrData&... amr_data) { + if constexpr (sizeof...(AmrData) == 1) { + if constexpr (std::is_same_v, Element>>) { + if (((mesh == amr_data.first) and ...)) { + // This element hasn't changed during AMR. Nothing to do. + return; + } + } + } + // Background fields in the volume using factory_classes = typename std::decay_t::factory_creation::factory_classes; From 7ae1f5ed183cce089a63c3b84c565dd762b63d48 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Fri, 26 Jan 2024 11:44:47 -0800 Subject: [PATCH 8/9] Add AMR projector for elliptic DG operator --- .../Actions/ApplyOperator.hpp | 68 ++++++++++++++----- .../DiscontinuousGalerkin/CMakeLists.txt | 1 + src/Elliptic/Executables/Solver.hpp | 19 +++--- .../Test_SubdomainOperator.cpp | 6 +- .../DiscontinuousGalerkin/Test_DgOperator.cpp | 5 +- 5 files changed, 67 insertions(+), 32 deletions(-) diff --git a/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp index ac3cfc591b1e..4f7365da03a2 100644 --- a/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp +++ b/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp @@ -37,6 +37,7 @@ #include "Parallel/GlobalCache.hpp" #include "Parallel/InboxInserters.hpp" #include "Parallel/Invoke.hpp" +#include "ParallelAlgorithms/Amr/Projectors/DefaultInitialize.hpp" #include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/GetOutput.hpp" #include "Utilities/Gsl.hpp" @@ -207,13 +208,16 @@ struct PrepareAndSendMortarData< // Can't `db::get` the arguments for the boundary conditions within // `db::mutate`, so we extract the data to mutate and move it back in when // we're done. + // Possible optimization: also keep memory for the mortar data around in the + // DataBox. Currently the mortar data is extracted and erased by + // `apply_operator` anyway, so we can just create a new map here to avoid + // dealing with AMR resizing the mortars. When we keep the memory around, + // its size has to be adjusted when the mesh changes during AMR. typename PrimalFluxesTag::type primal_fluxes; - typename all_mortar_data_tag::type all_mortar_data; - db::mutate( - [&primal_fluxes, &all_mortar_data](const auto local_primal_fluxes, - const auto local_all_mortar_data) { + typename all_mortar_data_tag::type all_mortar_data{}; + db::mutate( + [&primal_fluxes](const auto local_primal_fluxes) { primal_fluxes = std::move(*local_primal_fluxes); - all_mortar_data = std::move(*local_all_mortar_data); }, make_not_null(&box)); @@ -446,13 +450,27 @@ template using initialize_operator = tmpl::list< detail::InitializeFacesMortarsAndBackground>; +/*! + * \brief AMR projectors for the tags added by `initialize_operator` + */ +template +using amr_projectors = tmpl::append< + tmpl::list>, + tmpl::conditional_t< + std::is_same_v>, + tmpl::list<>, + tmpl::list>>>; + /*! * \brief Apply the DG operator to the `PrimalFieldsTag` and write the result to * the `OperatorAppliedToFieldsTag` * - * Add this list to the action list of a parallel component to compute the - * elliptic DG operator or its linearization. The operator involves a - * communication between nearest-neighbor elements. See `elliptic::dg` for + * Add the `apply_actions` list to the action list of a parallel component to + * compute the elliptic DG operator or its linearization. The operator involves + * a communication between nearest-neighbor elements. See `elliptic::dg` for * details on the elliptic DG operator. Make sure to add * `elliptic::dg::Actions::initialize_operator` to the initialization phase of * your parallel component so the required DataBox tags are set up before @@ -468,21 +486,37 @@ using initialize_operator = tmpl::list< * example when applying the nonlinear and linearized operator. They default to * the `PrimalFieldsTag` and the `PrimalFluxesTag`, meaning memory buffers * corresponding to these tags are set up in the DataBox. + * + * \par AMR + * Also add the `amr_projectors` to the list of AMR projectors to support AMR. */ template -using apply_operator = - tmpl::list, - detail::ReceiveMortarDataAndApplyOperator< - System, Linearized, TemporalIdTag, PrimalFieldsTag, - PrimalFluxesTag, OperatorAppliedToFieldsTag, - PrimalMortarFieldsTag, PrimalMortarFluxesTag>>; +struct DgOperator { + private: + static constexpr size_t Dim = System::volume_dim; + + public: + using apply_actions = + tmpl::list, + detail::ReceiveMortarDataAndApplyOperator< + System, Linearized, TemporalIdTag, PrimalFieldsTag, + PrimalFluxesTag, OperatorAppliedToFieldsTag, + PrimalMortarFieldsTag, PrimalMortarFluxesTag>>; + using amr_projectors = tmpl::list<::amr::projectors::DefaultInitialize< + PrimalFluxesTag, OperatorAppliedToFieldsTag, + ::Tags::Mortars, + Dim>>>; +}; /*! * \brief For linear systems, impose inhomogeneous boundary conditions as diff --git a/src/Elliptic/DiscontinuousGalerkin/CMakeLists.txt b/src/Elliptic/DiscontinuousGalerkin/CMakeLists.txt index 45326f61ea1f..4bfbb59294d3 100644 --- a/src/Elliptic/DiscontinuousGalerkin/CMakeLists.txt +++ b/src/Elliptic/DiscontinuousGalerkin/CMakeLists.txt @@ -33,6 +33,7 @@ target_link_libraries( ErrorHandling FunctionsOfTime LinearOperators + ParallelAmr Spectral Utilities INTERFACE diff --git a/src/Elliptic/Executables/Solver.hpp b/src/Elliptic/Executables/Solver.hpp index 985385de867a..7034141b2582 100644 --- a/src/Elliptic/Executables/Solver.hpp +++ b/src/Elliptic/Executables/Solver.hpp @@ -232,7 +232,7 @@ struct Solver { typename schwarz_smoother::register_element>>; template - using build_operator_actions = elliptic::dg::Actions::apply_operator< + using dg_operator = elliptic::dg::Actions::DgOperator< system, Linearized, tmpl::conditional_t, @@ -243,14 +243,13 @@ struct Solver { using build_matrix_actions = LinearSolver::Actions::build_matrix_actions< linear_solver_iteration_id, vars_tag, operator_applied_to_vars_tag, - build_operator_actions, + typename dg_operator::apply_actions, domain::Tags::Coordinates, LinearSolver::multigrid::Tags::IsFinestGrid>; template - using smooth_actions = - typename schwarz_smoother::template solve, - Label>; + using smooth_actions = typename schwarz_smoother::template solve< + typename dg_operator::apply_actions, Label>; /// This data needs to be communicated on subdomain overlap regions using communicated_overlap_tags = tmpl::flatten, + typename dg_operator::apply_actions, // Schwarz smoothing on each multigrid level smooth_actions, smooth_actions>, // Support disabling the preconditioner ::LinearSolver::Actions::make_identity_if_skipped< - multigrid, build_operator_actions>>, + multigrid, typename dg_operator::apply_actions>>, StepActions>; template using nonlinear_solve_actions = typename nonlinear_solver::template solve< - build_operator_actions, + typename dg_operator::apply_actions, tmpl::list< // Transfer data down the multigrid hierachy LinearSolver::multigrid::Actions::ReceiveFieldsFromFinerGrid< @@ -307,10 +306,10 @@ struct Solver { // Linear solve tmpl::list< // Apply the DG operator to the initial guess - elliptic::dg::Actions::apply_operator< + 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>, + fluxes_vars_tag>::apply_actions, // Modify fixed sources with boundary conditions elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource< system, fixed_sources_tag>, diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp b/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp index 507969062b77..7e2ae25519c8 100644 --- a/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp +++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp @@ -374,9 +374,9 @@ struct ElementArray { SubdomainOperatorAppliedToDataTag; using apply_full_dg_operator_actions = - ::elliptic::dg::Actions::apply_operator; + typename ::elliptic::dg::Actions::DgOperator< + System, true, TemporalIdTag, fields_tag, fluxes_tag, + operator_applied_to_fields_tag>::apply_actions; using background_tag = elliptic::Tags::Background; diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp index 5e6cd94243de..2db773d408c3 100644 --- a/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp +++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp @@ -145,9 +145,10 @@ struct ElementArray { ImposeInhomogeneousBoundaryConditionsOnSource< System, fixed_sources_tag>, ::Actions::Label, - ::elliptic::dg::Actions::apply_operator< + typename ::elliptic::dg::Actions::DgOperator< System, Linearized, TemporalIdTag, vars_tag, - primal_fluxes_vars_tag, operator_applied_to_vars_tag>, + primal_fluxes_vars_tag, + operator_applied_to_vars_tag>::apply_actions, IncrementTemporalId, Parallel::Actions::TerminatePhase>>>; }; From 6e95e1f0334e80d4dd72a1c675a3f56950d5571d Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Fri, 16 Feb 2024 17:12:31 -0800 Subject: [PATCH 9/9] Add test for elliptic AMR projectors --- .../DiscontinuousGalerkin/CMakeLists.txt | 2 + .../DiscontinuousGalerkin/Test_DgOperator.cpp | 249 +++++++++++++++--- 2 files changed, 208 insertions(+), 43 deletions(-) diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/CMakeLists.txt b/tests/Unit/Elliptic/DiscontinuousGalerkin/CMakeLists.txt index 204400e538c5..3efe6fde9355 100644 --- a/tests/Unit/Elliptic/DiscontinuousGalerkin/CMakeLists.txt +++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/CMakeLists.txt @@ -14,6 +14,7 @@ add_test_library(${LIBRARY} "${LIBRARY_SOURCES}") target_link_libraries( ${LIBRARY} PRIVATE + AmrCriteria AnalyticSolutions DataStructures DiscontinuousGalerkin @@ -24,6 +25,7 @@ target_link_libraries( EllipticDg ErrorHandling Parallel + ParallelAmr Poisson PoissonSolutions Spectral diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp index 2db773d408c3..2ef4749fc5c2 100644 --- a/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp +++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp @@ -4,6 +4,7 @@ #include "Framework/TestingFramework.hpp" #include +#include #include #include #include @@ -43,8 +44,16 @@ #include "Parallel/AlgorithmExecution.hpp" #include "Parallel/Phase.hpp" #include "ParallelAlgorithms/Actions/Goto.hpp" +#include "ParallelAlgorithms/Actions/InitializeItems.hpp" #include "ParallelAlgorithms/Actions/SetData.hpp" #include "ParallelAlgorithms/Actions/TerminatePhase.hpp" +#include "ParallelAlgorithms/Amr/Actions/Component.hpp" +#include "ParallelAlgorithms/Amr/Actions/EvaluateRefinementCriteria.hpp" +#include "ParallelAlgorithms/Amr/Actions/Initialize.hpp" +#include "ParallelAlgorithms/Amr/Criteria/Random.hpp" +#include "ParallelAlgorithms/Amr/Policies/Isotropy.hpp" +#include "ParallelAlgorithms/Amr/Policies/Tags.hpp" +#include "ParallelAlgorithms/Amr/Protocols/AmrMetavariables.hpp" #include "PointwiseFunctions/AnalyticSolutions/Poisson/ProductOfSinusoids.hpp" #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp" #include "Utilities/CartesianProduct.hpp" @@ -115,6 +124,10 @@ struct ElementArray { using primal_fluxes_vars_tag = ::Tags::Variables; using operator_applied_to_vars_tag = ::Tags::Variables>; + using dg_operator = + ::elliptic::dg::Actions::DgOperator; // Don't wrap the fixed sources in the `Var` prefix because typically we want // to impose inhomogeneous boundary conditions on the un-prefixed vars, i.e. // not necessarily the vars we apply the linearized operator to @@ -138,6 +151,8 @@ struct ElementArray { ::elliptic::Actions::InitializeFixedSources< System, analytic_solution_tag>, ::elliptic::dg::Actions::initialize_operator, + Initialization::Actions::InitializeItems< + ::amr::Initialization::Initialize>, Parallel::Actions::TerminatePhase>>, Parallel::PhaseActions< Parallel::Phase::Testing, @@ -145,18 +160,65 @@ struct ElementArray { ImposeInhomogeneousBoundaryConditionsOnSource< System, fixed_sources_tag>, ::Actions::Label, - typename ::elliptic::dg::Actions::DgOperator< - System, Linearized, TemporalIdTag, vars_tag, - primal_fluxes_vars_tag, - operator_applied_to_vars_tag>::apply_actions, - IncrementTemporalId, Parallel::Actions::TerminatePhase>>>; + typename dg_operator::apply_actions, IncrementTemporalId, + Parallel::Actions::TerminatePhase>>>; +}; + +template +struct AmrComponent { + using metavariables = Metavariables; + using chare_type = ActionTesting::MockSingletonChare; + using array_index = int; + using component_being_mocked = ::amr::Component; + using const_global_cache_tags = + tmpl::list<::amr::Criteria::Tags::Criteria, ::amr::Tags::Policies, + logging::Tags::Verbosity>; + + using phase_dependent_action_list = tmpl::list>>>>; +}; + +template +struct ProjectTemporalId : tt::ConformsTo<::amr::protocols::Projector> { + using return_tags = + tmpl::list>; + using argument_tags = tmpl::list<>; + // p-refinement + template + static void apply( + const gsl::not_null /*temporal_id*/, + const gsl::not_null**> /*cache*/, + const std::pair, Element>& /*old_mesh_and_element*/) {} + // h-refinement + template + static void apply( + const gsl::not_null temporal_id, + const gsl::not_null**> /*cache*/, + const tuples::TaggedTuple& parent_items) { + *temporal_id = get(parent_items); + } + // h-coarsening + template + static void apply( + const gsl::not_null temporal_id, + const gsl::not_null**> /*cache*/, + const std::unordered_map< + ElementId, tuples::TaggedTuple>& children_items) { + *temporal_id = get(children_items.begin()->second); + } }; template struct Metavariables { + static constexpr size_t volume_dim = System::volume_dim; using analytic_solution = AnalyticSolution; using element_array = ElementArray; - using component_list = tmpl::list; + using amr_component = AmrComponent; + using component_list = tmpl::list; using const_global_cache_tags = tmpl::list<::Tags::AnalyticSolution>; struct factory_creation @@ -167,8 +229,31 @@ struct Metavariables { tmpl::list>>, tmpl::pair>, - tmpl::pair>>; + tmpl::pair>, + tmpl::pair<::amr::Criterion, tmpl::list<::amr::Criteria::Random>>>; + }; + struct amr : tt::ConformsTo<::amr::protocols::AmrMetavariables> { + using projectors = tmpl::flatten, + ::amr::projectors::DefaultInitialize< + tmpl::list, + domain::Tags::InitialRefinementLevels, + typename element_array::vars_tag>>, + elliptic::dg::ProjectGeometry, + elliptic::Actions::InitializeFixedSources< + System, typename element_array::analytic_solution_tag>, + ::elliptic::Actions::InitializeOptionalAnalyticSolution< + volume_dim, typename element_array::analytic_solution_tag, + tmpl::append, + analytic_solution>, + elliptic::dg::Actions::amr_projectors< + System, typename element_array::analytic_solution_tag>, + typename element_array::dg_operator::amr_projectors>>; }; + + // NOLINTNEXTLINE(google-runtime-references) + void pup(PUP::er& /*p*/) {} }; template < @@ -192,7 +277,8 @@ void test_dg_operator( std::unordered_map< ElementId, typename ElementArray::operator_applied_to_vars_tag::type>>>& - tests_data) { + tests_data, + const bool test_amr = false) { CAPTURE(penalty_parameter); CAPTURE(use_massive_dg_operator); @@ -222,20 +308,33 @@ void test_dg_operator( } } + // Configure AMR criteria so we _increase_ (never decrease) resolution + std::vector> amr_criteria{}; + amr_criteria.push_back(std::make_unique<::amr::Criteria::Random>( + std::unordered_map<::amr::Flag, size_t>{ + // h-refinement is not supported yet in the action testing framework + {::amr::Flag::Split, 0}, + {::amr::Flag::IncreaseResolution, 1}, + {::amr::Flag::DoNothing, 1}})); + ActionTesting::MockRuntimeSystem runner{tuples::TaggedTuple< domain::Tags::Domain, domain::Tags::FunctionsOfTimeInitialize, domain::Tags::ExternalBoundaryConditions, ::elliptic::dg::Tags::PenaltyParameter, ::elliptic::dg::Tags::Massive, ::elliptic::dg::Tags::Quadrature, ::elliptic::dg::Tags::Formulation, - ::Tags::AnalyticSolution>{ + ::Tags::AnalyticSolution, + ::amr::Criteria::Tags::Criteria, ::amr::Tags::Policies, + logging::Tags::Verbosity<::amr::OptionTags::AmrGroup>>{ std::move(domain), domain_creator.functions_of_time(), std::move(boundary_conditions), penalty_parameter, use_massive_dg_operator, quadrature, ::dg::Formulation::StrongInertial, - analytic_solution}}; + analytic_solution, std::move(amr_criteria), + ::amr::Policies{::amr::Isotropy::Anisotropic}, ::Verbosity::Debug}}; // DataBox shortcuts - const auto get_tag = [&runner]( - auto tag_v, const ElementId& local_element_id) -> const auto& { + const auto get_tag = + [&runner](auto tag_v, + const ElementId& local_element_id) -> const auto& { using tag = std::decay_t; return ActionTesting::get_databox_tag(runner, local_element_id); @@ -379,36 +478,100 @@ void test_dg_operator( apply_operator_and_check_result({}, all_zero_primal_fluxes, all_zero_operator_vars); } - { - INFO("Test A(x) = b with analytic solution"); - std::unordered_map, Vars> analytic_primal_vars{}; - std::unordered_map, PrimalFluxesVars> - analytic_primal_fluxes{}; - std::unordered_map, OperatorAppliedToVars> - analytic_fixed_source_with_inhom_bc{}; - for (const auto& element_id : all_element_ids) { - const auto& inertial_coords = get_tag( - domain::Tags::Coordinates{}, element_id); - analytic_primal_vars[element_id] = - variables_from_tagged_tuple(analytic_solution.variables( - inertial_coords, typename System::primal_fields{})); - analytic_primal_fluxes[element_id] = - variables_from_tagged_tuple(analytic_solution.variables( - inertial_coords, typename System::primal_fluxes{})); - analytic_fixed_source_with_inhom_bc[element_id] = - get_tag(fixed_sources_tag{}, element_id); - } - apply_operator_and_check_result( - analytic_primal_vars, analytic_primal_fluxes, - analytic_fixed_source_with_inhom_bc, analytic_solution_aux_approx, - analytic_solution_operator_approx); - } + const auto test_analytic_solution = + [&analytic_solution, &all_element_ids, &get_tag, + &apply_operator_and_check_result, &analytic_solution_aux_approx, + &analytic_solution_operator_approx]() { + INFO("Test A(x) = b with analytic solution"); + std::unordered_map, Vars> analytic_primal_vars{}; + std::unordered_map, PrimalFluxesVars> + analytic_primal_fluxes{}; + std::unordered_map, OperatorAppliedToVars> + analytic_fixed_source_with_inhom_bc{}; + for (const auto& element_id : all_element_ids) { + const auto& inertial_coords = get_tag( + domain::Tags::Coordinates{}, element_id); + analytic_primal_vars[element_id] = + variables_from_tagged_tuple(analytic_solution.variables( + inertial_coords, typename System::primal_fields{})); + analytic_primal_fluxes[element_id] = + variables_from_tagged_tuple(analytic_solution.variables( + inertial_coords, typename System::primal_fluxes{})); + analytic_fixed_source_with_inhom_bc[element_id] = + get_tag(fixed_sources_tag{}, element_id); + } + apply_operator_and_check_result( + analytic_primal_vars, analytic_primal_fluxes, + analytic_fixed_source_with_inhom_bc, analytic_solution_aux_approx, + analytic_solution_operator_approx); + }; + test_analytic_solution(); { INFO("Test A(x) = b with custom x and b"); for (const auto& test_data : tests_data) { std::apply(apply_operator_and_check_result, test_data); } } + + // Stop here if we're not testing AMR + if (not test_amr) { + return; + } + + INFO("Run AMR!"); + const auto invoke_all_simple_actions = [&runner, &all_element_ids]() { + bool quiescence = false; + while (not quiescence) { + quiescence = true; + for (const auto& element_id : all_element_ids) { + while (not ActionTesting::is_simple_action_queue_empty( + runner, element_id)) { + ActionTesting::invoke_queued_simple_action( + make_not_null(&runner), element_id); + quiescence = false; + } + } + } + }; + using amr_component = typename Metavars::amr_component; + ActionTesting::emplace_singleton_component( + make_not_null(&runner), ActionTesting::NodeId{0}, + ActionTesting::LocalCoreId{0}); + auto& cache = ActionTesting::cache(runner, 0); + Parallel::simple_action<::amr::Actions::EvaluateRefinementCriteria>( + Parallel::get_parallel_component(cache)); + invoke_all_simple_actions(); + Parallel::simple_action<::amr::Actions::AdjustDomain>( + Parallel::get_parallel_component(cache)); + invoke_all_simple_actions(); + // h-refinement is not supported yet in the action testing framework + // // Invoke `CreateChild` on AMR component + // while (not ActionTesting::is_simple_action_queue_empty( + // runner, 0)) { + // ActionTesting::invoke_queued_simple_action( + // make_not_null(&runner), 0); + // } + // // Update list of element IDs + // all_element_ids.clear(); + // for (const auto& [element_id, element] : + // runner.template mock_distributed_objects()) { + // (void)element; + // all_element_ids.insert(element_id); + // } + // CAPTURE(all_element_ids); + // Go to the start of the testing phase to run the + // `ImposeInhomogeneousBoundaryConditionsOnSource` action again + ActionTesting::set_phase(make_not_null(&runner), Parallel::Phase::Testing); + if constexpr (Linearized) { + for (const auto& element_id : all_element_ids) { + ActionTesting::next_action(make_not_null(&runner), + element_id); + } + } + + // Test analytic solution again after AMR. It should still pass because we + // only refined, never coarsened. + test_analytic_solution(); } } // namespace @@ -561,7 +724,7 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") { test_dg_operator( domain_creator, penalty_parameter, use_massive_dg_operator, quadrature, analytic_solution, analytic_solution_aux_approx, - analytic_solution_operator_approx, {}); + analytic_solution_operator_approx, {}, true); } } } @@ -668,7 +831,7 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") { test_dg_operator( domain_creator, penalty_parameter, use_massive_dg_operator, quadrature, analytic_solution, analytic_solution_aux_approx, - analytic_solution_operator_approx, {}); + analytic_solution_operator_approx, {}, true); } } } @@ -808,7 +971,7 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") { test_dg_operator( domain_creator, penalty_parameter, use_massive_dg_operator, quadrature, analytic_solution, analytic_solution_aux_approx, - analytic_solution_operator_approx, {}); + analytic_solution_operator_approx, {}, true); } } } @@ -934,10 +1097,10 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") { cartesian_product(make_array(true, false), make_array(Spectral::Quadrature::GaussLobatto, Spectral::Quadrature::Gauss))) { - test_dg_operator(domain_creator, penalty_parameter, - massive, quadrature, analytic_solution, - analytic_solution_aux_approx, - analytic_solution_operator_approx, {}); + test_dg_operator( + domain_creator, penalty_parameter, massive, quadrature, + analytic_solution, analytic_solution_aux_approx, + analytic_solution_operator_approx, {}, true); } } }