From b74d7338810cf8141dd3f704979fca8f07c92b91 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Mon, 29 Jan 2024 00:23:03 -0800 Subject: [PATCH] Re-initialize Schwarz subdomain before each solve To account for changes in domain geometry (by AMR). --- .../DiscontinuousGalerkin/Initialization.hpp | 3 + .../SubdomainOperator/InitializeSubdomain.hpp | 37 +- src/Elliptic/Executables/Solver.hpp | 64 ++- .../Actions/CommunicateOverlapFields.hpp | 35 +- .../Test_SubdomainOperator.cpp | 437 ++++++++++++------ .../Actions/Test_CommunicateOverlapFields.cpp | 3 +- 6 files changed, 399 insertions(+), 180 deletions(-) diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp index fd602faffd97..68f6068592db 100644 --- a/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp +++ b/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp @@ -103,6 +103,7 @@ struct InitializeGeometry { domain::Tags::InitialRefinementLevels, domain::Tags::Domain, domain::Tags::FunctionsOfTime, elliptic::dg::Tags::Quadrature>; + using volume_tags = argument_tags; static void apply( gsl::not_null*> mesh, gsl::not_null*> element, gsl::not_null>*> neighbor_meshes, @@ -140,6 +141,8 @@ struct ProjectGeometry : tt::ConformsTo<::amr::protocols::Projector> { using argument_tags = tmpl::list, domain::Tags::Element, domain::Tags::Domain, domain::Tags::FunctionsOfTime>; + using volume_tags = + tmpl::list, domain::Tags::FunctionsOfTime>; // p-refinement static void apply( diff --git a/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/InitializeSubdomain.hpp b/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/InitializeSubdomain.hpp index 46daa9ac6bb3..52daba375eec 100644 --- a/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/InitializeSubdomain.hpp +++ b/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/InitializeSubdomain.hpp @@ -107,19 +107,36 @@ struct InitializeOverlapGeometry { * Note that the geometry depends on the system and on the choice of background * through the normalization of face normals, which involves a metric. * + * The `FromInitialDomain` template parameter controls how the + * next-to-nearest-neighbor information is acquired: + * - If `FromInitialDomain = true`, uses the initial parameters from the input + * file. This can be used at the beginning of a simulation or if the domain + * never changes (no AMR). + * - If `FromInitialDomain = false`, uses the mesh, element, and neighbor meshes + * of each neighbor. This typically requires communication to exchange this + * data between neighbors each time the domain changes, e.g. with AMR. + * * DataBox: * - Uses: * - `BackgroundTag` * - `domain::Tags::Element` - * - `domain::Tags::InitialExtents` - * - `domain::Tags::InitialRefinementLevels` + * - If `FromInitialDomain = true`: + * - `domain::Tags::InitialExtents` + * - `domain::Tags::InitialRefinementLevels` + * - `elliptic::dg::Tags::Quadrature` + * - If `FromInitialDomain = false`: + * - `Overlaps>` + * - `Overlaps>` + * - `Overlaps>` * - `domain::Tags::Domain` + * - `domain::Tags::FunctionsOfTime` * - `LinearSolver::Schwarz::Tags::MaxOverlap` * - Adds: Many tags prefixed with `LinearSolver::Schwarz::Tags::Overlaps`. See * `elliptic::dg::Actions::InitializeDomain` and * `elliptic::dg::Actions::initialize_operator` for a complete list. */ -template +template struct InitializeSubdomain { private: static constexpr size_t Dim = System::volume_dim; @@ -128,7 +145,10 @@ struct InitializeSubdomain { static constexpr bool has_background_fields = not std::is_same_v>; - using InitializeGeometry = elliptic::dg::InitializeGeometry; + using InitializeGeometry = + tmpl::conditional_t, + elliptic::dg::ProjectGeometry>; using InitializeOverlapGeometry = detail::InitializeOverlapGeometry; using InitializeFacesAndMortars = elliptic::dg::InitializeFacesAndMortars< Dim, typename System::inv_metric_tag, BackgroundTag>; @@ -186,8 +206,13 @@ struct InitializeSubdomain { elliptic::util::mutate_apply_at< db::wrap_tags_in, - typename InitializeGeometry::argument_tags, - typename InitializeGeometry::argument_tags>( + tmpl::append< + db::wrap_tags_in>, + typename InitializeGeometry::volume_tags>, + typename InitializeGeometry::volume_tags>( InitializeGeometry{}, make_not_null(&box), overlap_id, neighbor_id); // Initialize subdomain-specific tags on overlaps elliptic::util::mutate_apply_at< diff --git a/src/Elliptic/Executables/Solver.hpp b/src/Elliptic/Executables/Solver.hpp index 7034141b2582..5ede602b73b9 100644 --- a/src/Elliptic/Executables/Solver.hpp +++ b/src/Elliptic/Executables/Solver.hpp @@ -213,8 +213,6 @@ struct Solver { typename system::primal_fluxes>, elliptic::analytic_data::AnalyticSolution>, elliptic::dg::Actions::initialize_operator, - elliptic::dg::subdomain_operator::Actions::InitializeSubdomain< - system, background_tag, typename schwarz_smoother::options_group>, tmpl::conditional_t< is_linear, tmpl::list<>, ::Initialization::Actions::AddComputeTags::apply_actions, Label>; + // These tags are communicated on subdomain overlaps to initialize the + // subdomain geometry. AMR updates these tags, so we have to communicate them + // after each AMR step. + using subdomain_init_tags = + tmpl::list, + domain::Tags::Element, + domain::Tags::NeighborMesh>; + /// This data needs to be communicated on subdomain overlap regions using communicated_overlap_tags = tmpl::flatten>>>; + using init_subdomain_action = + elliptic::dg::subdomain_operator::Actions::InitializeSubdomain< + system, background_tag, typename schwarz_smoother::options_group, + false>; + template using linear_solve_actions = typename linear_solver::template solve< tmpl::list< @@ -289,10 +300,12 @@ struct Solver { // Communicate data on subdomain overlap regions LinearSolver::Schwarz::Actions::SendOverlapFields< communicated_overlap_tags, - typename schwarz_smoother::options_group, false>, + typename schwarz_smoother::options_group, false, + nonlinear_solver_iteration_id>, LinearSolver::Schwarz::Actions::ReceiveOverlapFields< volume_dim, communicated_overlap_tags, - typename schwarz_smoother::options_group>, + typename schwarz_smoother::options_group, false, + nonlinear_solver_iteration_id>, // Reset Schwarz subdomain solver LinearSolver::Schwarz::Actions::ResetSubdomainSolver< typename schwarz_smoother::options_group>, @@ -301,22 +314,33 @@ struct Solver { StepActions>; template - using solve_actions = tmpl::conditional_t< - is_linear, - // Linear solve - tmpl::list< - // Apply the DG operator to the initial guess - 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>::apply_actions, - // Modify fixed sources with boundary conditions - elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource< - system, fixed_sources_tag>, - // Krylov solve - linear_solve_actions>, - // Nonlinear solve - nonlinear_solve_actions>; + using solve_actions = tmpl::flatten, + LinearSolver::Schwarz::Actions::ReceiveOverlapFields< + volume_dim, subdomain_init_tags, + typename schwarz_smoother::options_group, false>, + init_subdomain_action, + // Linear or nonlinear solve + tmpl::conditional_t< + is_linear, + // Linear solve + tmpl::list< + // Apply the DG operator to the initial guess + 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>::apply_actions, + // Modify fixed sources with boundary conditions + elliptic::dg::Actions:: + ImposeInhomogeneousBoundaryConditionsOnSource< + system, fixed_sources_tag>, + // Krylov solve + linear_solve_actions>, + // Nonlinear solve + nonlinear_solve_actions>>>; using component_list = tmpl::flatten< tmpl::list, diff --git a/src/ParallelAlgorithms/LinearSolver/Schwarz/Actions/CommunicateOverlapFields.hpp b/src/ParallelAlgorithms/LinearSolver/Schwarz/Actions/CommunicateOverlapFields.hpp index 1019e7af1d7d..b5a601c23bc5 100644 --- a/src/ParallelAlgorithms/LinearSolver/Schwarz/Actions/CommunicateOverlapFields.hpp +++ b/src/ParallelAlgorithms/LinearSolver/Schwarz/Actions/CommunicateOverlapFields.hpp @@ -70,14 +70,15 @@ struct OverlapFieldsTag * This actions should be followed by * `LinearSolver::Schwarz::Actions::ReceiveOverlapFields` in the action list. */ -template +template > struct SendOverlapFields; /// \cond template + bool RestrictToOverlap, typename TemporalIdTag> struct SendOverlapFields, OptionsGroup, - RestrictToOverlap> { + RestrictToOverlap, TemporalIdTag> { using const_global_cache_tags = tmpl::list, logging::Tags::Verbosity>; @@ -92,15 +93,16 @@ struct SendOverlapFields, OptionsGroup, const ParallelComponent* const /*meta*/) { const auto& element = get>(box); - // Skip communicating if the overlap is empty - if (UNLIKELY(db::get>(box) == 0 or + // Skip communicating if the overlap is empty and we RestrictToOverlap, + // because we would end up with empty tensor data. + if (UNLIKELY((RestrictToOverlap and + db::get>(box) == 0) or element.number_of_neighbors() == 0)) { return {Parallel::AlgorithmExecution::Continue, std::nullopt}; } // Do some logging - const auto& iteration_id = - get>(box); + const auto& iteration_id = get(box); if (UNLIKELY(get>(box) >= ::Verbosity::Debug)) { Parallel::printf("%s %s(%zu): Send overlap fields\n", element_id, @@ -172,12 +174,16 @@ struct SendOverlapFields, OptionsGroup, * This actions should be preceded by * `LinearSolver::Schwarz::Actions::SendOverlapFields` in the action list. */ -template +template > struct ReceiveOverlapFields; /// \cond -template -struct ReceiveOverlapFields, OptionsGroup> { +template +struct ReceiveOverlapFields, OptionsGroup, + RestrictToOverlap, TemporalIdTag> { private: using overlap_fields_tag = detail::OverlapFieldsTag, OptionsGroup>; @@ -197,12 +203,12 @@ struct ReceiveOverlapFields, OptionsGroup> { const Parallel::GlobalCache& /*cache*/, const ElementId& element_id, const ActionList /*meta*/, const ParallelComponent* const /*meta*/) { - const auto& iteration_id = - get>(box); + const auto& iteration_id = get(box); const auto& element = get>(box); - // Nothing to receive if overlap is empty - if (UNLIKELY(db::get>(box) == 0 or + // Nothing to receive if nothing was sent in the action above + if (UNLIKELY((RestrictToOverlap and + db::get>(box) == 0) or element.number_of_neighbors() == 0)) { return {Parallel::AlgorithmExecution::Continue, std::nullopt}; } @@ -227,6 +233,7 @@ struct ReceiveOverlapFields, OptionsGroup> { db::mutate...>( [&received_overlap_fields](const auto... local_overlap_fields) { if constexpr (sizeof...(OverlapFields) > 1) { + (local_overlap_fields->clear(), ...); for (auto& [overlap_id, overlap_fields] : received_overlap_fields) { expand_pack((*local_overlap_fields)[overlap_id] = std::move(get(overlap_fields))...); diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp b/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp index 7e2ae25519c8..fc3c0907c6a8 100644 --- a/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp +++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/SubdomainOperator/Test_SubdomainOperator.cpp @@ -60,8 +60,18 @@ #include "Parallel/AlgorithmExecution.hpp" #include "Parallel/Phase.hpp" #include "Parallel/PhaseDependentActionList.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/Projectors/CopyFromCreatorOrLeaveAsIs.hpp" +#include "ParallelAlgorithms/Amr/Protocols/AmrMetavariables.hpp" +#include "ParallelAlgorithms/LinearSolver/Schwarz/Actions/CommunicateOverlapFields.hpp" #include "ParallelAlgorithms/LinearSolver/Schwarz/ElementCenteredSubdomainData.hpp" #include "ParallelAlgorithms/LinearSolver/Schwarz/OverlapHelpers.hpp" #include "ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp" @@ -373,10 +383,10 @@ struct ElementArray { using subdomain_operator_applied_to_fields_tag = SubdomainOperatorAppliedToDataTag; - using apply_full_dg_operator_actions = - typename ::elliptic::dg::Actions::DgOperator< - System, true, TemporalIdTag, fields_tag, fluxes_tag, - operator_applied_to_fields_tag>::apply_actions; + using dg_operator = + ::elliptic::dg::Actions::DgOperator; using background_tag = elliptic::Tags::Background; @@ -386,38 +396,70 @@ struct ElementArray { using array_index = ElementId; using const_global_cache_tags = tmpl::list, background_tag>; + + // These tags are communicated on subdomain overlaps to initialize the + // subdomain geometry. AMR updates these tags, so we have to communicate them + // after each AMR step. + using subdomain_init_tags = + tmpl::list, domain::Tags::Element, + domain::Tags::NeighborMesh>; + using init_subdomain_action = + ::elliptic::dg::subdomain_operator::Actions::InitializeSubdomain< + System, background_tag, DummyOptionsGroup, false>; + using init_random_subdomain_data_action = + InitializeRandomSubdomainData; + using phase_dependent_action_list = tmpl::list< Parallel::PhaseActions< Parallel::Phase::Initialization, - tmpl::list< - ActionTesting::InitializeDataBox< - tmpl::list, - domain::Tags::InitialExtents, - SubdomainOperatorTag, - subdomain_operator_applied_to_fields_tag, - OverrideBoundaryConditionsTag>>, - ::elliptic::dg::Actions::InitializeDomain, - ::elliptic::dg::Actions::initialize_operator, - ::elliptic::dg::subdomain_operator::Actions::InitializeSubdomain< - System, - tmpl::conditional_t, - DummyOptionsGroup>, - InitializeRandomSubdomainData, - ExtraInitActions, Parallel::Actions::TerminatePhase>>, + tmpl::list, + domain::Tags::InitialExtents, + SubdomainOperatorTag, + subdomain_operator_applied_to_fields_tag, + OverrideBoundaryConditionsTag>>, + ::elliptic::dg::Actions::InitializeDomain, + ::elliptic::dg::Actions::initialize_operator< + System, background_tag>, + Initialization::Actions::InitializeItems< + ::amr::Initialization::Initialize>, + ExtraInitActions, Parallel::Actions::TerminatePhase>>, Parallel::PhaseActions< Parallel::Phase::Testing, - tmpl::list, - TestSubdomainOperatorMatrix< - SubdomainOperator, typename fields_tag::tags_list>, - Parallel::Actions::TerminatePhase>>>; + tmpl::list< + // Init subdomain + LinearSolver::Schwarz::Actions::SendOverlapFields< + subdomain_init_tags, DummyOptionsGroup, false, TemporalIdTag>, + LinearSolver::Schwarz::Actions::ReceiveOverlapFields< + Dim, subdomain_init_tags, DummyOptionsGroup, false, + TemporalIdTag>, + init_subdomain_action, init_random_subdomain_data_action, + Parallel::Actions::TerminatePhase, + // Full DG operator + typename dg_operator::apply_actions, + Parallel::Actions::TerminatePhase, + // Subdomain operator + ApplySubdomainOperator, + TestSubdomainOperatorMatrix, + 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<::amr::OptionTags::AmrGroup>>; + + using phase_dependent_action_list = tmpl::list>>>>; }; template ; - using component_list = tmpl::list; + using amr_component = AmrComponent; + using component_list = tmpl::list; using const_global_cache_tags = tmpl::list<>; + static constexpr size_t volume_dim = System::volume_dim; struct factory_creation : tt::ConformsTo { using factory_classes = tmpl::map< tmpl::pair< - elliptic::BoundaryConditions::BoundaryCondition, + elliptic::BoundaryConditions::BoundaryCondition, tmpl::list>>, tmpl::pair>>>; + tmpl::list>>, + tmpl::pair<::amr::Criterion, tmpl::list<::amr::Criteria::Random>>>; }; + template + using overlaps_tag = + LinearSolver::Schwarz::Tags::Overlaps; + struct amr : tt::ConformsTo<::amr::protocols::AmrMetavariables> { + using projectors = tmpl::flatten, + domain::Tags::InitialRefinementLevels, + SubdomainOperatorTag, + typename element_array:: + subdomain_operator_applied_to_fields_tag, + typename element_array::fields_tag>, + // Tags communicated on subdomain overlaps. No need to project + // these during AMR because they will be communicated. + db::wrap_tags_in, + // Tags initialized on subdomains. No need to project these during + // AMR because they will get re-initialized after communication. + typename element_array::init_subdomain_action::simple_tags, + typename element_array::init_random_subdomain_data_action:: + simple_tags>>, + ::amr::projectors::CopyFromCreatorOrLeaveAsIs< + OverrideBoundaryConditionsTag, TemporalIdTag, + // Work around a segfault because this tag isn't handled + // correctly by the testing framework + Parallel::Tags::GlobalCacheImpl>, + elliptic::dg::ProjectGeometry, + elliptic::dg::Actions::amr_projectors< + System, typename element_array::background_tag>, + typename element_array::dg_operator::amr_projectors, ExtraInitActions>>; + }; + + // NOLINTNEXTLINE(google-runtime-references) + void pup(PUP::er& /*p*/) {} }; // The test should work for any elliptic system. For systems with fluxes or @@ -490,18 +569,30 @@ void test_subdomain_operator( const auto element_ids = ::initial_element_ids(initial_ref_levs); CAPTURE(element_ids.size()); + // Configure AMR criteria + 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::IncreaseResolution, 1}, + {::amr::Flag::DoNothing, 1}})); + ActionTesting::MockRuntimeSystem runner{tuples::TaggedTuple< domain::Tags::Domain, domain::Tags::FunctionsOfTimeInitialize, domain::Tags::ExternalBoundaryConditions, elliptic::Tags::Background, LinearSolver::Schwarz::Tags::MaxOverlap, + logging::Tags::Verbosity, elliptic::dg::Tags::PenaltyParameter, elliptic::dg::Tags::Massive, - elliptic::dg::Tags::Quadrature, elliptic::dg::Tags::Formulation>{ + elliptic::dg::Tags::Quadrature, elliptic::dg::Tags::Formulation, + ::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), - std::make_unique>(), overlap, penalty_parameter, - use_massive_dg_operator, quadrature, - ::dg::Formulation::StrongInertial}}; + std::make_unique>(), overlap, + ::Verbosity::Verbose, penalty_parameter, use_massive_dg_operator, + quadrature, ::dg::Formulation::StrongInertial, std::move(amr_criteria), + ::amr::Policies{::amr::Isotropy::Anisotropic}, ::Verbosity::Debug}}; // Initialize all elements, generating random subdomain data for (const auto& element_id : element_ids) { @@ -517,7 +608,6 @@ void test_subdomain_operator( element_id); } } - ActionTesting::set_phase(make_not_null(&runner), Parallel::Phase::Testing); // DataBox shortcuts const auto get_tag = [&runner](const ElementId& local_element_id, auto tag_v) -> decltype(auto) { @@ -533,125 +623,181 @@ void test_subdomain_operator( make_not_null(&runner), local_element_id, value); }; - // For selection of expensive tests - std::uniform_int_distribution dist_select_subdomain_center( - 0, element_ids.size() - 1); - const size_t rnd_subdomain_center = dist_select_subdomain_center(gen); - size_t subdomain_center_id = 0; - - // Take each element as the subdomain-center in turn - for (const auto& subdomain_center : element_ids) { - CAPTURE(subdomain_center); - - // First, reset the data on all elements to zero + const auto test_subdomain_operator_equals_dg_operator = [&runner, + &element_ids, + &overlap, + &rnd_overlap, &gen, + &get_tag, + &set_tag]() { + ActionTesting::set_phase(make_not_null(&runner), + Parallel::Phase::Testing); for (const auto& element_id : element_ids) { - set_tag(element_id, fields_tag{}, - typename fields_tag::type{ - get_tag(element_id, domain::Tags::Mesh{}) - .number_of_grid_points(), - 0.}); - } - - // Set data on the central element and its neighbors to the subdomain data - const auto& subdomain_data = - get_tag(subdomain_center, subdomain_data_tag{}); - const auto& all_overlap_extents = - get_tag(subdomain_center, - LinearSolver::Schwarz::Tags::Overlaps< - elliptic::dg::subdomain_operator::Tags::ExtrudingExtent, - Dim, DummyOptionsGroup>{}); - const auto& central_element = - get_tag(subdomain_center, domain::Tags::Element{}); - set_tag(subdomain_center, fields_tag{}, subdomain_data.element_data); - for (const auto& [overlap_id, overlap_data] : - subdomain_data.overlap_data) { - const auto& [direction, neighbor_id] = overlap_id; - const auto direction_from_neighbor = - central_element.neighbors().at(direction).orientation()( - direction.opposite()); - set_tag( - neighbor_id, fields_tag{}, - LinearSolver::Schwarz::extended_overlap_data( - overlap_data, - get_tag(neighbor_id, domain::Tags::Mesh{}).extents(), - all_overlap_extents.at(overlap_id), direction_from_neighbor)); + ActionTesting::next_action_if_ready( + make_not_null(&runner), element_id); } - - // Run actions to compute the full DG-operator for (const auto& element_id : element_ids) { - CAPTURE(element_id); - runner.template mock_distributed_objects() - .at(element_id) - .force_next_action_to_be(0); - runner.template mock_distributed_objects() - .at(element_id) - .set_terminate(false); while (not ActionTesting::get_terminate(runner, element_id) and ActionTesting::next_action_if_ready( make_not_null(&runner), element_id)) { } } - // Break here so all elements have sent mortar data before receiving it - for (const auto& element_id : element_ids) { - CAPTURE(element_id); - while (not ActionTesting::get_terminate(runner, - element_id) and - ActionTesting::next_action_if_ready( - make_not_null(&runner), element_id)) { + + // For selection of expensive tests + std::uniform_int_distribution dist_select_subdomain_center( + 0, element_ids.size() - 1); + const auto& rnd_subdomain_center = dist_select_subdomain_center(gen); + size_t subdomain_center_id = 0; + + // Take each element as the subdomain-center in turn + for (const auto& subdomain_center : element_ids) { + CAPTURE(subdomain_center); + + // First, reset the data on all elements to zero + for (const auto& element_id : element_ids) { + set_tag(element_id, fields_tag{}, + typename fields_tag::type{ + get_tag(element_id, domain::Tags::Mesh{}) + .number_of_grid_points(), + 0.}); } - } - // Invoke ApplySubdomainOperator action only on the subdomain center - ActionTesting::next_action(make_not_null(&runner), - subdomain_center); - - // Test that the subdomain operator and the full DG-operator computed the - // same result within the subdomain - const auto& subdomain_result = - get_tag(subdomain_center, subdomain_operator_applied_to_fields_tag{}); - Approx custom_approx = Approx::custom().epsilon(1.e-12).scale(1.0); - CHECK_VARIABLES_CUSTOM_APPROX( - subdomain_result.element_data, - get_tag(subdomain_center, operator_applied_to_fields_tag{}), - custom_approx); - REQUIRE(subdomain_result.overlap_data.size() == - subdomain_data.overlap_data.size()); - for (const auto& [overlap_id, overlap_result] : - subdomain_result.overlap_data) { - CAPTURE(overlap_id); - const auto& [direction, neighbor_id] = overlap_id; - const auto direction_from_neighbor = - central_element.neighbors().at(direction).orientation()( - direction.opposite()); - const auto expected_overlap_result = - LinearSolver::Schwarz::data_on_overlap( - get_tag(neighbor_id, operator_applied_to_fields_tag{}), - get_tag(neighbor_id, domain::Tags::Mesh{}).extents(), - all_overlap_extents.at(overlap_id), direction_from_neighbor); - CHECK_VARIABLES_CUSTOM_APPROX(overlap_result, expected_overlap_result, - custom_approx); - } + // Set data on the central element and its neighbors to the subdomain + // data + const auto& subdomain_data = + get_tag(subdomain_center, subdomain_data_tag{}); + const auto& all_overlap_extents = + get_tag(subdomain_center, + LinearSolver::Schwarz::Tags::Overlaps< + elliptic::dg::subdomain_operator::Tags::ExtrudingExtent, + Dim, DummyOptionsGroup>{}); + const auto& central_element = + get_tag(subdomain_center, domain::Tags::Element{}); + set_tag(subdomain_center, fields_tag{}, subdomain_data.element_data); + for (const auto& [overlap_id, overlap_data] : + subdomain_data.overlap_data) { + const auto& [direction, neighbor_id] = overlap_id; + const auto direction_from_neighbor = + central_element.neighbors().at(direction).orientation()( + direction.opposite()); + set_tag( + neighbor_id, fields_tag{}, + LinearSolver::Schwarz::extended_overlap_data( + overlap_data, + get_tag(neighbor_id, domain::Tags::Mesh{}).extents(), + all_overlap_extents.at(overlap_id), direction_from_neighbor)); + } + + // Run actions to compute the full DG-operator + for (const auto& element_id : element_ids) { + CAPTURE(element_id); + runner.template mock_distributed_objects() + .at(element_id) + .force_next_action_to_be(5); + runner.template mock_distributed_objects() + .at(element_id) + .set_terminate(false); + while (not ActionTesting::get_terminate(runner, + element_id) and + ActionTesting::next_action_if_ready( + make_not_null(&runner), element_id)) { + } + } + // Break here so all elements have sent mortar data before receiving it + for (const auto& element_id : element_ids) { + CAPTURE(element_id); + while (not ActionTesting::get_terminate(runner, + element_id) and + ActionTesting::next_action_if_ready( + make_not_null(&runner), element_id)) { + } + } - // Now build the matrix representation of the subdomain operator - // explicitly, and apply it to the data to make sure the matrix is - // equivalent to the matrix-free operator. This is important to test - // because the subdomain operator includes optimizations for when it is - // invoked on sparse data, i.e. data that is mostly zero, which is the - // case when building it explicitly column-by-column. - if (overlap == rnd_overlap and - subdomain_center_id == rnd_subdomain_center) { + // Invoke ApplySubdomainOperator action only on the subdomain center ActionTesting::next_action(make_not_null(&runner), subdomain_center); + + // Test that the subdomain operator and the full DG-operator computed + // the same result within the subdomain + const auto& subdomain_result = get_tag( + subdomain_center, subdomain_operator_applied_to_fields_tag{}); + Approx custom_approx = Approx::custom().epsilon(1.e-12).scale(1.0); + CHECK_VARIABLES_CUSTOM_APPROX( + subdomain_result.element_data, + get_tag(subdomain_center, operator_applied_to_fields_tag{}), + custom_approx); + REQUIRE(subdomain_result.overlap_data.size() == + subdomain_data.overlap_data.size()); + for (const auto& [overlap_id, overlap_result] : + subdomain_result.overlap_data) { + CAPTURE(overlap_id); + const auto& [direction, neighbor_id] = overlap_id; + const auto direction_from_neighbor = + central_element.neighbors().at(direction).orientation()( + direction.opposite()); + const auto expected_overlap_result = + LinearSolver::Schwarz::data_on_overlap( + get_tag(neighbor_id, operator_applied_to_fields_tag{}), + get_tag(neighbor_id, domain::Tags::Mesh{}).extents(), + all_overlap_extents.at(overlap_id), direction_from_neighbor); + CHECK_VARIABLES_CUSTOM_APPROX(overlap_result, expected_overlap_result, + custom_approx); + } + + // Now build the matrix representation of the subdomain operator + // explicitly, and apply it to the data to make sure the matrix is + // equivalent to the matrix-free operator. This is important to test + // because the subdomain operator includes optimizations for when it is + // invoked on sparse data, i.e. data that is mostly zero, which is the + // case when building it explicitly column-by-column. + if (overlap == rnd_overlap and + subdomain_center_id == rnd_subdomain_center) { + ActionTesting::next_action(make_not_null(&runner), + subdomain_center); + } + ++subdomain_center_id; + } // loop over subdomain centers + }; + + test_subdomain_operator_equals_dg_operator(); + + INFO("Run AMR!"); + const auto invoke_all_simple_actions = [&runner, &element_ids]() { + bool quiescence = false; + while (not quiescence) { + quiescence = true; + for (const auto& element_id : 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; + } + } } - ++subdomain_center_id; - } // loop over subdomain centers - } // loop over overlaps + }; + using amr_component = typename metavariables::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(); + + // Test again after AMR + test_subdomain_operator_equals_dg_operator(); + } // loop over overlaps } // Add a constitutive relation for elasticity systems to the DataBox template -struct InitializeConstitutiveRelation { +struct InitializeConstitutiveRelation + : tt::ConformsTo<::amr::protocols::Projector> { + public: // Iterative action using simple_tags = tmpl::list>; template >; + using argument_tags = tmpl::list<>; + template + static void apply( + const gsl::not_null>*> + material, + const AmrData&... /*amr_data*/) { + *material = std::make_unique< + Elasticity::ConstitutiveRelations::IsotropicHomogeneous>(1., 2.); + } }; } // namespace diff --git a/tests/Unit/ParallelAlgorithms/LinearSolver/Schwarz/Actions/Test_CommunicateOverlapFields.cpp b/tests/Unit/ParallelAlgorithms/LinearSolver/Schwarz/Actions/Test_CommunicateOverlapFields.cpp index 0540827f273d..87d4400ec733 100644 --- a/tests/Unit/ParallelAlgorithms/LinearSolver/Schwarz/Actions/Test_CommunicateOverlapFields.cpp +++ b/tests/Unit/ParallelAlgorithms/LinearSolver/Schwarz/Actions/Test_CommunicateOverlapFields.cpp @@ -64,7 +64,8 @@ struct ElementArray { LinearSolver::Schwarz::Actions::SendOverlapFields< tmpl::list, DummyOptionsGroup, RestrictToOverlap>, LinearSolver::Schwarz::Actions::ReceiveOverlapFields< - Dim, tmpl::list, DummyOptionsGroup>, + Dim, tmpl::list, DummyOptionsGroup, + RestrictToOverlap>, Parallel::Actions::TerminatePhase>>>; };