From 0e38ae42be9e9c050a295d44ec2a70f9102d4f92 Mon Sep 17 00:00:00 2001 From: William Throwe Date: Thu, 23 Feb 2023 15:06:35 -0500 Subject: [PATCH] Add IMEX SolveImplicitSector mutator --- .clang-tidy | 8 +- src/DataStructures/Tensor/TypeAliases.hpp | 6 + src/Evolution/Imex/CMakeLists.txt | 4 + .../Imex/Protocols/ImplicitSector.hpp | 8 + src/Evolution/Imex/SolveImplicitSector.hpp | 118 +++ src/Evolution/Imex/SolveImplicitSector.tpp | 520 ++++++++++++ tests/Unit/Evolution/Imex/CMakeLists.txt | 2 + .../Imex/Test_SolveImplicitSector.cpp | 749 ++++++++++++++++++ tests/Unit/Helpers/Evolution/CMakeLists.txt | 1 + .../Helpers/Evolution/Imex/CMakeLists.txt | 21 + .../Helpers/Evolution/Imex/TestSector.hpp | 209 +++++ 11 files changed, 1643 insertions(+), 3 deletions(-) create mode 100644 src/Evolution/Imex/SolveImplicitSector.hpp create mode 100644 src/Evolution/Imex/SolveImplicitSector.tpp create mode 100644 tests/Unit/Evolution/Imex/Test_SolveImplicitSector.cpp create mode 100644 tests/Unit/Helpers/Evolution/Imex/CMakeLists.txt create mode 100644 tests/Unit/Helpers/Evolution/Imex/TestSector.hpp diff --git a/.clang-tidy b/.clang-tidy index 78c08a289f2d5..f8a9a14fa2c95 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -84,12 +84,14 @@ Checks: '*, # we are okay with lower case, -readability-uppercase-literal-suffix,' CheckOptions: - - key: performance-move-const-arg.CheckTriviallyCopyableMove - value: false - - key: cppcoreguidelines-special-member-functions.AllowSoleDefaultDtor + - key: cppcoreguidelines-avoid-do-while.IgnoreMacros value: true - key: cppcoreguidelines-special-member-functions.AllowMissingMoveFunctions value: true + - key: cppcoreguidelines-special-member-functions.AllowSoleDefaultDtor + value: true + - key: performance-move-const-arg.CheckTriviallyCopyableMove + value: false WarningsAsErrors: '*' # It is unclear if the header filter actually works or how to use it so # just include all headers diff --git a/src/DataStructures/Tensor/TypeAliases.hpp b/src/DataStructures/Tensor/TypeAliases.hpp index f319c99a726c0..647bd601c4ad8 100644 --- a/src/DataStructures/Tensor/TypeAliases.hpp +++ b/src/DataStructures/Tensor/TypeAliases.hpp @@ -399,6 +399,12 @@ using iJkk = Tensor, SpatialIndex, SpatialIndex>>; +template +using iiJJ = Tensor, + index_list, + SpatialIndex, + SpatialIndex, + SpatialIndex>>; } // namespace tnsr /*! diff --git a/src/Evolution/Imex/CMakeLists.txt b/src/Evolution/Imex/CMakeLists.txt index bdc5e5290dc55..15b023f69bff6 100644 --- a/src/Evolution/Imex/CMakeLists.txt +++ b/src/Evolution/Imex/CMakeLists.txt @@ -19,6 +19,8 @@ spectre_target_headers( GuessResult.hpp Mode.hpp NamespaceDocs.hpp + SolveImplicitSector.hpp + SolveImplicitSector.tpp ) target_link_libraries( @@ -28,6 +30,8 @@ target_link_libraries( Options INTERFACE DataStructures + LinearSolver + RootFinding Time Utilities ) diff --git a/src/Evolution/Imex/Protocols/ImplicitSector.hpp b/src/Evolution/Imex/Protocols/ImplicitSector.hpp index d6d279864886d..8623d365b8b41 100644 --- a/src/Evolution/Imex/Protocols/ImplicitSector.hpp +++ b/src/Evolution/Imex/Protocols/ImplicitSector.hpp @@ -86,6 +86,14 @@ namespace imex::protocols { /// /// All `Variables` in the DataBox, including the sources and source /// jacobian, will be initialized to zero with a single grid point. +/// +/// \snippet Test_SolveImplicitSector.cpp ImplicitSector +/// \snippet Test_SolveImplicitSector.cpp initial_guess +/// +/// Examples of definitions of the implicit source and jacobian: +/// +/// \snippet Test_SolveImplicitSector.cpp source +/// \snippet Test_SolveImplicitSector.cpp jacobian struct ImplicitSector { template struct test { diff --git a/src/Evolution/Imex/SolveImplicitSector.hpp b/src/Evolution/Imex/SolveImplicitSector.hpp new file mode 100644 index 0000000000000..4e423d66b6cec --- /dev/null +++ b/src/Evolution/Imex/SolveImplicitSector.hpp @@ -0,0 +1,118 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include + +#include "DataStructures/Tensor/TypeAliases.hpp" +#include "Evolution/Imex/Mode.hpp" +#include "Evolution/Imex/Protocols/ImplicitSector.hpp" +#include "Time/Tags/TimeStepper.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/ProtocolHelpers.hpp" +#include "Utilities/TMPL.hpp" + +/// \cond +class DataVector; +class ImexTimeStepper; +class TimeDelta; +template +class Variables; +namespace Tags { +struct TimeStep; +} // namespace Tags +namespace TimeSteppers { +template +class History; +} // namespace TimeSteppers +namespace imex::Tags { +template +struct ImplicitHistory; +struct Mode; +template +struct SolveFailures; +struct SolveTolerance; +} // namespace imex::Tags +/// \endcond + +namespace imex { +namespace solve_implicit_sector_detail { +template +using ForwardTuple = tmpl::wrap< + tmpl::transform>>>, + std::tuple>; +} // namespace solve_implicit_sector_detail + +/// Perform the implicit solve for one implicit sector. +/// +/// This will update the tensors in the implicit sector and clean up +/// the corresponding time stepper history. A new history entry is +/// not added, because that should be done with the same values of the +/// variables used for the explicit portion of the time derivative, +/// which may still undergo variable-fixing-like corrections. +/// +/// \warning +/// This will use the value of `::Tags::Time` from the DataBox. Most +/// of the time, the value appropriate for evaluating the explicit RHS +/// is stored there, so it will likely need to be set to the +/// appropriate value for the implicit RHS for the duration of this +/// mutation. +template +struct SolveImplicitSector { + static_assert( + tt::assert_conforms_to_v); + + public: + using SystemVariables = typename SystemVariablesTag::type; + using SectorVariables = Variables; + + private: + template + struct get_tags_from_evolution { + using type = typename Attempt::tags_from_evolution; + }; + + using attempt_tags = tmpl::transform>; + using evolution_data_tags = + tmpl::push_front; + + using EvolutionDataTuple = solve_implicit_sector_detail::ForwardTuple< + tmpl::join>; + + static void apply_impl( + gsl::not_null system_variables, + gsl::not_null*> implicit_history, + gsl::not_null*> solve_failures, + const ImexTimeStepper& time_stepper, const TimeDelta& time_step, + Mode implicit_solve_mode, double implicit_solve_tolerance, + const EvolutionDataTuple& joined_evolution_data); + + public: + using return_tags = tmpl::list, + Tags::SolveFailures>; + using argument_tags = + tmpl::append, ::Tags::TimeStep, + Tags::Mode, Tags::SolveTolerance>, + tmpl::join>; + + template + static void apply(const gsl::not_null system_variables, + const gsl::not_null*> + implicit_history, + const gsl::not_null*> solve_failures, + const ImexTimeStepper& time_stepper, + const TimeDelta& time_step, const Mode implicit_solve_mode, + const double implicit_solve_tolerance, + const ForwardArgs&... forward_args) { + apply_impl(system_variables, implicit_history, solve_failures, time_stepper, + time_step, implicit_solve_mode, implicit_solve_tolerance, + std::forward_as_tuple(forward_args...)); + } +}; +} // namespace imex diff --git a/src/Evolution/Imex/SolveImplicitSector.tpp b/src/Evolution/Imex/SolveImplicitSector.tpp new file mode 100644 index 0000000000000..4e1de1680efa1 --- /dev/null +++ b/src/Evolution/Imex/SolveImplicitSector.tpp @@ -0,0 +1,520 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include "Evolution/Imex/SolveImplicitSector.hpp" + +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataBox/PrefixHelpers.hpp" +#include "DataStructures/DataBox/Prefixes.hpp" +#include "DataStructures/DataBox/Tag.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/ExtractPoint.hpp" +#include "DataStructures/Matrix.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "DataStructures/VariablesTag.hpp" +#include "Evolution/Imex/GuessResult.hpp" +#include "Evolution/Imex/Mode.hpp" +#include "Evolution/Imex/Protocols/ImplicitSector.hpp" +#include "Evolution/Imex/Tags/Jacobian.hpp" +#include "NumericalAlgorithms/LinearSolver/Lapack.hpp" +#include "NumericalAlgorithms/RootFinding/GslMultiRoot.hpp" +#include "Time/History.hpp" +#include "Time/TimeSteppers/ImexTimeStepper.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/ErrorHandling/Error.hpp" +#include "Utilities/ErrorHandling/Exceptions.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/ProtocolHelpers.hpp" +#include "Utilities/SplitTuple.hpp" +#include "Utilities/StdArrayHelpers.hpp" +#include "Utilities/TMPL.hpp" +#include "Utilities/TaggedTuple.hpp" +#include "Utilities/TypeTraits/IsA.hpp" + +/// \cond +class TimeDelta; +/// \endcond + +namespace imex { +namespace solve_implicit_sector_detail { +template +class ImplicitEquation { + public: + template + ImplicitEquation( + const gsl::not_null system_variables, + const gsl::not_null*> + implicit_history, + const ImexTimeStepper& time_stepper, const TimeDelta& time_step, + const ForwardTuple& + initial_guess_arguments, + tmpl::type_ /*meta*/) + : implicit_weight_( + time_stepper.implicit_weight(implicit_history, time_step)) { + static_assert(std::is_same_v); + + if (implicit_weight_ == 0.0) { + // Explicit substep. No solves. Just write the correct answer + // into the evolution box. + auto sector_subset = system_variables->template reference_subset< + typename SectorVariables::tags_list>(); + time_stepper.add_inhomogeneous_implicit_terms( + make_not_null(§or_subset), implicit_history, time_step); + return; + } + + inhomogeneous_terms_ = + system_variables + ->template extract_subset(); + time_stepper.add_inhomogeneous_implicit_terms( + make_not_null(&inhomogeneous_terms_), implicit_history, time_step); + + { + const auto initial_guess_return = + tmpl::as_pack( + [&](auto... tags) { + return std::make_tuple( + make_not_null(&get>( + *system_variables))...); + }); + initial_guess_types_ = std::apply( + [&](const auto&... args) { + return ImplicitSector::initial_guess::apply( + args..., inhomogeneous_terms_, implicit_weight_); + }, + std::tuple_cat(initial_guess_return, initial_guess_arguments)); + } + + ASSERT(initial_guess_types_.size() == + system_variables->number_of_grid_points() or + initial_guess_types_.empty(), + "initial_guess must return one GuessResult per point or an empty " + "vector for GuessResult::InitialGuess everywhere."); + } + + double implicit_weight() const { return implicit_weight_; } + const SectorVariables& inhomogeneous_terms() const { + return inhomogeneous_terms_; + } + + bool do_solve() const { return implicit_weight_ != 0.0; } + GuessResult initial_guess_result(const size_t point_index) const { + return initial_guess_types_.empty() ? GuessResult::InitialGuess + : initial_guess_types_[point_index]; + } + + private: + double implicit_weight_; + SectorVariables inhomogeneous_terms_{}; + std::vector initial_guess_types_{}; +}; + +template +class ImplicitSolver { + static_assert( + tt::assert_conforms_to_v); + + using sector_variables_tag = + ::Tags::Variables; + using SectorVariables = typename sector_variables_tag::type; + static constexpr size_t solve_dimension = + SectorVariables::number_of_independent_components; + + using tags_from_evolution = typename SolveAttempt::tags_from_evolution; + using EvolutionData = ForwardTuple; + + struct EvolutionDataTag : db::SimpleTag { + using type = const EvolutionData*; + }; + + struct SolverPointIndex : db::SimpleTag { + using type = size_t; + }; + + template + struct FromEvolution : Tag, db::ReferenceTag { + using base = Tag; + using parent_tag = EvolutionDataTag; + using argument_tags = tmpl::list; + static const typename base::type& get( + const EvolutionData* const evolution_data) { + return std::get::value>( + *evolution_data); + } + }; + + template + struct FromEvolution> : Tag, db::ComputeTag { + using base = Tag; + using argument_tags = tmpl::list; + static constexpr auto function( + const gsl::not_null result, + const EvolutionData* const evolution_data, const size_t index) { + result->initialize(1); + extract_point( + result, + get::value>(*evolution_data), + index); + } + }; + + // Tensor always allocates, so instead of creating one + // we create a one-tensor Variables, and the DataBox will allow + // access to the tensor transparently. We could instead manually + // set all the DataVectors as non-owning, pointing at individual + // doubles, but that's more work and it's not clear it gains us + // anything. + template + struct FromEvolution> + : ::Tags::Variables>, db::ComputeTag { + using base = ::Tags::Variables>; + using argument_tags = tmpl::list; + static constexpr auto function( + const gsl::not_null result, + const EvolutionData* const evolution_data, const size_t index) { + result->initialize(1); + extract_point( + make_not_null(&get(*result)), + get::value>(*evolution_data), + index); + } + }; + + using all_mutators = tmpl::remove_duplicates< + tmpl::append, + typename SolveAttempt::source_prep, + typename SolveAttempt::jacobian_prep>>; + + using source_tag = db::add_tag_prefix<::Tags::Source, sector_variables_tag>; + using jacobian_tag = + ::Tags::Variables>; + + using internal_simple_tags = + tmpl::list; + using wrapped_tags_from_evolution = + tmpl::transform>; + + using simple_tags = + tmpl::append; + using compute_tags = tmpl::append; + + using SolveBox = + db::compute_databox_type>; + + public: + template + ImplicitSolver(const ImplicitEquation& implicit_equation, + PassedEvolutionData&& data_from_evolution) + : solve_box_(db::create()), + implicit_equation_(&implicit_equation) { + static_assert(std::is_same_v, + "ImplicitSolver was passed a temporary."); + db::mutate_apply< + tmpl::push_front< + tmpl::filter< + simple_tags, + tt::is_a>>, + EvolutionDataTag>, + tmpl::list<>>( + [&data_from_evolution]( + const gsl::not_null evolution_data_pointer, + const auto... vars) { + *evolution_data_pointer = &data_from_evolution; + expand_pack((vars->initialize(1, 0.0), 0)...); + }, + make_not_null(&solve_box_)); + } + + void set_index(const size_t index) { + db::mutate( + [&index](const gsl::not_null box_index) { + *box_index = index; + }, + make_not_null(&solve_box_)); + + extract_point(make_not_null(&inhomogeneous_terms_), + implicit_equation_->inhomogeneous_terms(), index); + + completed_mutators_ = decltype(completed_mutators_){}; + } + + std::array operator()( + const std::array& sector_variables_array) const { + ASSERT(implicit_equation_->implicit_weight() != 0.0, + "Should not be performing solves on explicit substeps"); + set_sector_variables(sector_variables_array); + run_mutators>(); + std::array residual_array{}; + SectorVariables residual(residual_array.data(), residual_array.size()); + residual = + inhomogeneous_terms_ - db::get(solve_box_) + + implicit_equation_->implicit_weight() * + db::get>( + solve_box_); + return residual_array; + } + + std::array, solve_dimension> jacobian( + const std::array& sector_variables_array) const { + ASSERT(implicit_equation_->implicit_weight() != 0.0, + "Should not be performing solves on explicit substeps"); + set_sector_variables(sector_variables_array); + run_mutators>(); + + std::array, solve_dimension> + jacobian_array{}; + // The storage order for the tensors does not match the required + // order for the returned array, so we have to copy components + // individually. + // + // Despite repeated references to then, the result of this is + // independent of the *_for_offsets variables. They are only used + // for calculating offsets into the returned array. + const auto& variables_for_offsets = + db::get(solve_box_); + tmpl::for_each( + [&](auto dependent_tag_v) { + using dependent_tag = tmpl::type_from; + const auto& dependent_for_offsets = + get(variables_for_offsets); + for (size_t dependent_component = 0; + dependent_component < dependent_for_offsets.size(); + ++dependent_component) { + const auto dependent_index = + dependent_for_offsets.get_tensor_index(dependent_component); + auto& result_row = jacobian_array[static_cast( + dependent_for_offsets[dependent_component].data() - + variables_for_offsets.data())]; + tmpl::for_each( + [&](auto independent_tag_v) { + using independent_tag = + tmpl::type_from; + using jacobian_component_tag = + imex::Tags::Jacobian>; + const auto& independent_for_offsets = + get(variables_for_offsets); + + for (size_t independent_component = 0; + independent_component < independent_for_offsets.size(); + ++independent_component) { + const auto independent_index = + independent_for_offsets.get_tensor_index( + independent_component); + result_row[static_cast( + independent_for_offsets[independent_component].data() - + variables_for_offsets.data())] = + get(solve_box_) + .get(concatenate(independent_index, + dependent_index))[0]; + } + }); + } + }); + + jacobian_array *= implicit_equation_->implicit_weight(); + + for (size_t i = 0; i < solve_dimension; ++i) { + jacobian_array[i][i] -= 1.0; + } + return jacobian_array; + } + + private: + void set_sector_variables( + std::array sector_variables_array) const { + const SectorVariables sector_variables(sector_variables_array.data(), + sector_variables_array.size()); + set_sector_variables(sector_variables); + } + + void set_sector_variables(const SectorVariables& sector_variables) const { + if (sector_variables == most_recent_sector_variables_) { + return; + } + most_recent_sector_variables_ = sector_variables; + db::mutate( + [§or_variables](const gsl::not_null vars) { + *vars = sector_variables; + }, + make_not_null(&solve_box_)); + completed_mutators_ = decltype(completed_mutators_){}; + } + + template + void run_mutators() const { + tmpl::for_each([this](auto mutator_v) { + using mutator = tmpl::type_from; + if (not get>(completed_mutators_)) { + db::mutate_apply(make_not_null(&solve_box_)); + get>(completed_mutators_) = true; + } + }); + } + + template + struct RanMutator { + using type = bool; + }; + + // Re mutables: This struct is only used locally in serial + // single-threaded implicit solves. The gsl_multiroot interface + // takes a const solver object, but we want to be able to share + // calculations between the source and jacobian calculations. + // NOLINTNEXTLINE(spectre-mutable) + mutable SolveBox solve_box_; + SectorVariables inhomogeneous_terms_{1}; + gsl::not_null*> implicit_equation_; + // NOLINTNEXTLINE(spectre-mutable) + mutable SectorVariables most_recent_sector_variables_{}; + // NOLINTNEXTLINE(spectre-mutable) + mutable tuples::tagged_tuple_from_typelist< + tmpl::transform>> + completed_mutators_{}; +}; +} // namespace solve_implicit_sector_detail + +template +void SolveImplicitSector::apply_impl( + const gsl::not_null system_variables, + const gsl::not_null*> + implicit_history, + const gsl::not_null*> solve_failures, + const ImexTimeStepper& time_stepper, const TimeDelta& time_step, + const Mode implicit_solve_mode, const double implicit_solve_tolerance, + const EvolutionDataTuple& joined_evolution_data) { + get(*solve_failures) = 0.0; + + const auto evolution_data = split_tuple< + tmpl::transform>>( + joined_evolution_data); + const auto& initial_guess_arguments = std::get<0>(evolution_data); + + const solve_implicit_sector_detail::ImplicitEquation + equation(system_variables, implicit_history, time_stepper, time_step, + initial_guess_arguments, tmpl::type_{}); + if (not equation.do_solve()) { + return; + } + + const size_t number_of_grid_points = get(*solve_failures).size(); + // Only allocated if used. + Matrix semi_implicit_jacobian{}; + + bool had_failure = true; + tmpl::for_each< + typename ImplicitSector::solve_attempts>([&](auto solve_attempt_v) { + using solve_attempt = tmpl::type_from; + if (not had_failure) { + return; + } + had_failure = false; + constexpr bool have_fallback = + not std::is_same_v>; + constexpr auto attempt_number = + tmpl::index_of::value; + // Entry 0 is for the initial guess. + const auto& attempt_evolution_data = + std::get(evolution_data); + solve_implicit_sector_detail::ImplicitSolver + solver(equation, attempt_evolution_data); + + for (size_t point = 0; point < number_of_grid_points; ++point) { + if (get(*solve_failures)[point] < attempt_number) { + continue; + } + if (equation.initial_guess_result(point) == GuessResult::ExactSolution) { + // Initial guess was written into the evolution DataBox + // when it was computed. + continue; + } + + std::array + pointwise_vars_array; + SectorVariables pointwise_vars(pointwise_vars_array.data(), + pointwise_vars_array.size()); + solver.set_index(point); + std::array + initial_guess; + { + SectorVariables guess_vars(initial_guess.data(), initial_guess.size()); + extract_point(make_not_null(&guess_vars), + system_variables->template reference_subset< + typename SectorVariables::tags_list>(), + point); + } + switch (implicit_solve_mode) { + case Mode::Implicit: { + const size_t max_iterations = 100; + try { + pointwise_vars_array = RootFinder::gsl_multiroot( + solver, initial_guess, + RootFinder::StoppingConditions::Residual( + implicit_solve_tolerance), + max_iterations); + } catch (const convergence_error&) { + if constexpr (have_fallback) { + ++get(*solve_failures)[point]; + had_failure = true; + continue; + } else { + throw; + } + } + break; + } + case Mode::SemiImplicit: { + auto correction_array = solver(initial_guess); + DataVector correction(correction_array.data(), + correction_array.size()); + correction *= -1.0; + semi_implicit_jacobian = solver.jacobian(initial_guess); + const int lapack_info = lapack::general_matrix_linear_solve( + &correction, &semi_implicit_jacobian); + if (lapack_info != 0) { + if (lapack_info < 0) { + ERROR("LAPACK invalid argument: " << -lapack_info); + } else { + if constexpr (have_fallback) { + ++get(*solve_failures)[point]; + had_failure = true; + continue; + } else { + ERROR("Semi-implicit inversion was singular at\n" + << pointwise_vars); + } + } + } + pointwise_vars_array = initial_guess + correction_array; + break; + } + default: + ERROR("Invalid implicit mode"); + } + + // Write the result into the evolution variables. + auto sector_reference = system_variables->template reference_subset< + typename SectorVariables::tags_list>(); + overwrite_point(make_not_null(§or_reference), pointwise_vars, point); + } + }); +} +} // namespace imex diff --git a/tests/Unit/Evolution/Imex/CMakeLists.txt b/tests/Unit/Evolution/Imex/CMakeLists.txt index 0951fcbaa1103..d12d2754a203b 100644 --- a/tests/Unit/Evolution/Imex/CMakeLists.txt +++ b/tests/Unit/Evolution/Imex/CMakeLists.txt @@ -6,6 +6,7 @@ set(LIBRARY "Test_Imex") set(LIBRARY_SOURCES Test_GuessResult.cpp Test_Mode.cpp + Test_SolveImplicitSector.cpp ) add_subdirectory(Tags) @@ -19,6 +20,7 @@ target_link_libraries( DataStructuresHelpers Framework Imex + ImexHelpers Utilities ) diff --git a/tests/Unit/Evolution/Imex/Test_SolveImplicitSector.cpp b/tests/Unit/Evolution/Imex/Test_SolveImplicitSector.cpp new file mode 100644 index 0000000000000..5d56fa548130a --- /dev/null +++ b/tests/Unit/Evolution/Imex/Test_SolveImplicitSector.cpp @@ -0,0 +1,749 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataBox/PrefixHelpers.hpp" +#include "DataStructures/DataBox/Prefixes.hpp" +#include "DataStructures/DataBox/Tag.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "DataStructures/VariablesTag.hpp" +#include "Evolution/Imex/GuessResult.hpp" +#include "Evolution/Imex/Mode.hpp" +#include "Evolution/Imex/Protocols/ImplicitSector.hpp" +#include "Evolution/Imex/SolveImplicitSector.hpp" +#include "Evolution/Imex/SolveImplicitSector.tpp" +#include "Evolution/Imex/Tags/ImplicitHistory.hpp" +#include "Evolution/Imex/Tags/Jacobian.hpp" +#include "Evolution/Imex/Tags/Mode.hpp" +#include "Evolution/Imex/Tags/SolveFailures.hpp" +#include "Evolution/Imex/Tags/SolveTolerance.hpp" +#include "Framework/TestHelpers.hpp" +#include "Helpers/DataStructures/MakeWithRandomValues.hpp" +#include "Helpers/Evolution/Imex/TestSector.hpp" +#include "Time/History.hpp" +#include "Time/Slab.hpp" +#include "Time/Tags/TimeStep.hpp" +#include "Time/Tags/TimeStepper.hpp" +#include "Time/TimeStepId.hpp" +#include "Time/TimeSteppers/Heun2.hpp" +#include "Utilities/ErrorHandling/Error.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/Literals.hpp" +#include "Utilities/MakeWithValue.hpp" +#include "Utilities/ProtocolHelpers.hpp" +#include "Utilities/TMPL.hpp" +#include "Utilities/TaggedTuple.hpp" + +namespace { +// Set temporarily to verify that the solver correctly skips most of +// the work when the step is explicit. +// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) +bool performing_step_with_no_implicit_term = false; + +struct Var1 : db::SimpleTag { + using type = Scalar; +}; + +struct Var2 : db::SimpleTag { + using type = tnsr::II; +}; + +struct Var3 : db::SimpleTag { + using type = tnsr::I; +}; + +struct NonTensor : db::SimpleTag { + using type = double; +}; + +// These next several tags aren't used in the calculation, just for +// testing DataBox handling. +struct TensorFromEvolution : db::SimpleTag { + using type = Scalar; +}; + +using VariablesFromEvolution = Tags::Variables>; + +struct TensorTemporary : db::SimpleTag { + using type = Scalar; +}; + +struct VariablesTemporary : db::SimpleTag { + using type = Variables>; +}; + +struct SomeComputeTagBase : db::SimpleTag { + using type = double; +}; + +struct SomeComputeTag : SomeComputeTagBase, db::ComputeTag { + using base = SomeComputeTagBase; + using argument_tags = + tmpl::list; + static void function( + const gsl::not_null result, const Scalar& var1, + const Scalar& from_evolution, + const Variables>& temporary) { + // Check the initialization of the temporary Variables in the + // solver DataBox. None of the mutators modify the object, so it + // should always have that state. + CHECK(temporary.number_of_grid_points() == 1); + CHECK(get(get(temporary))[0] == 0.0); + + // Check slicing + CHECK(get(from_evolution).size() == 1); + CHECK(get(from_evolution)[0] == 2.0 * get(var1)[0]); + + *result = get(var1)[0] + 1.0; + } +}; + +enum class PrepId { Source, Jacobian, Shared }; + +struct RecordPreparersForTest : db::SimpleTag { + using type = std::array, 4>; +}; + +template +struct Preparer { + using return_tags = tmpl::list; + using argument_tags = tmpl::list; + + static void apply( + const gsl::not_null prep_run_values, + const Var2::type& var2, const Var3::type& var3) { + CHECK(not performing_step_with_no_implicit_term); + + std::pair current_values{var2, var3}; + CHECK((*prep_run_values)[static_cast(Prep)] != current_values); + (*prep_run_values)[static_cast(Prep)] = std::move(current_values); + } +}; +// End stuff only used for DataBox handling + +// [initial_guess] +struct AnalyticSolution { + using return_tags = tmpl::list; + using argument_tags = tmpl::list; + static std::vector apply( + const gsl::not_null*> var2, + const gsl::not_null*> var3, + const Scalar& var1, const double non_tensor, + const Variables>& inhomogeneous_terms, + const double implicit_weight) { + // Solution for source terms + // S[v2^ij] = v3^i v3^j - nt v2^ij + // S[v3^i] = -v1 v3^i + + // Solving v3^i = X - w v1 v3^i gives v3^i = X / (1 + w v1) + tenex::evaluate(var3, get(inhomogeneous_terms)(ti::I) / + (1.0 + implicit_weight * var1())); + tenex::evaluate( + var2, (get(inhomogeneous_terms)(ti::I, ti::J) + + implicit_weight * (*var3)(ti::I) * (*var3)(ti::J)) / + (1.0 + implicit_weight * non_tensor)); + return {get(var1).size(), imex::GuessResult::ExactSolution}; + } +}; +// [initial_guess] + +struct InitialGuess { + using return_tags = tmpl::list; + using argument_tags = tmpl::list<>; + static std::vector apply( + const gsl::not_null*> var2, + const gsl::not_null*> var3, + const Variables>& /*inhomogeneous_terms*/, + const double /*implicit_weight*/) { + CHECK(not performing_step_with_no_implicit_term); + + for (auto& component : *var2) { + component *= 2.0; + } + for (auto& component : *var3) { + component *= 3.0; + } + return {}; + } +}; + +// [source] +struct Source { + using return_tags = tmpl::list<::Tags::Source, ::Tags::Source>; + using argument_tags = tmpl::list; + + static void apply(const gsl::not_null*> source_var2, + const gsl::not_null*> source_var3, + const Scalar& var1, + const tnsr::II& var2, + const tnsr::I& var3, const double non_tensor, + const RecordPreparersForTest::type& prep_run_values, + const double compute_tag_value) { + // [source] + CHECK(not performing_step_with_no_implicit_term); + + const std::pair current_values{var2, var3}; + CHECK(prep_run_values[static_cast(PrepId::Shared)] == + current_values); + CHECK(prep_run_values[static_cast(PrepId::Source)] == + current_values); + + CHECK(compute_tag_value == get(var1)[0] + 1.0); + + work(source_var2, source_var3, var1, var2, var3, non_tensor); + } + + // Used in the test below. Not part of the IMEX interface. + static void work(const gsl::not_null*> source_var2, + const gsl::not_null*> source_var3, + const Scalar& var1, + const tnsr::II& var2, + const tnsr::I& var3, + const double non_tensor) { + tenex::evaluate( + source_var2, + var3(ti::I) * var3(ti::J) - non_tensor * var2(ti::I, ti::J)); + tenex::evaluate(source_var3, -var1() * var3(ti::I)); + } +}; + +// [jacobian] +struct Jacobian { + using return_tags = + tmpl::list>, + imex::Tags::Jacobian>, + imex::Tags::Jacobian>>; + using argument_tags = + tmpl::list; + + static void apply(const gsl::not_null*> dvar2_dvar2, + const gsl::not_null*> dvar2_dvar3, + const gsl::not_null*> dvar3_dvar3, + const Scalar& var1, + const tnsr::I& var3, const double non_tensor, + const RecordPreparersForTest::type& prep_run_values) { + // [jacobian] + CHECK(not performing_step_with_no_implicit_term); + + // We don't need var2 for anything else in this function. Hard to + // imagine a way not taking one of the variables as an argument + // could break anything, but easy to test so we do it. + CHECK(prep_run_values[static_cast(PrepId::Shared)].second == var3); + CHECK(prep_run_values[static_cast(PrepId::Jacobian)].second == + var3); + + std::fill(dvar2_dvar3->begin(), dvar2_dvar3->end(), 0.0); + for (size_t i = 0; i < 2; ++i) { + dvar2_dvar2->get(i, i, i, i) = -non_tensor; + dvar2_dvar3->get(i, i, i) = 2.0 * var3.get(i); + dvar3_dvar3->get(i, i) = -get(var1); + for (size_t j = 0; j < i; ++j) { + dvar2_dvar2->get(i, j, i, j) = -non_tensor; + dvar2_dvar3->get(i, i, j) += var3.get(j); + dvar2_dvar3->get(j, i, j) += var3.get(i); + } + } + } +}; + +// [ImplicitSector] +template +struct ImplicitSector : tt::ConformsTo { + using tensors = tmpl::list; + using initial_guess = tmpl::conditional_t; + + struct SolveAttempt { + using tags_from_evolution = + tmpl::list; + using simple_tags = tmpl::list; + using compute_tags = tmpl::list; + + using source_prep = + tmpl::list, Preparer>; + using jacobian_prep = + tmpl::list, Preparer>; + + using source = Source; + using jacobian = + tmpl::conditional_t; + }; + + using solve_attempts = tmpl::list; +}; +// [ImplicitSector] + +// ::tensors doesn't depend on the template parameter +using sector_variables_tag = Tags::Variables::tensors>; +using SectorVariables = sector_variables_tag::type; + +tuples::TaggedTuple +arbitrary_test_values() { + SectorVariables explicit_values(1); + tnsr::II& var2 = get(explicit_values); + get<0, 0>(var2) = 3.0; + get<0, 1>(var2) = 4.0; + get<1, 1>(var2) = 5.0; + tnsr::I& var3 = get(explicit_values); + get<0>(var3) = 6.0; + get<1>(var3) = 7.0; + + Scalar var1{}; + get(var1) = DataVector{8.0}; + const double non_tensor = 9.0; + Variables> test_variables(1); + get(test_variables)[0] = 2.0 * get(var1)[0]; + return {std::move(explicit_values), std::move(var1), non_tensor, + std::move(test_variables)}; +} + +template +void test_test_sector() { + using sector = ImplicitSector; + auto values = arbitrary_test_values(); + TestHelpers::imex::test_sector( + 1.0e-1, 1.0e-12, std::move(get(values)), + {std::move(get(values)), get(values), + std::move(get(values))}); +} + +void test_internal_jacobian_ordering() { + // This test doesn't make sense on the analytic solution version. + using sector = ImplicitSector; + using solve_attempt = sector::SolveAttempt; + + const Slab slab(0.0, 1.0); + TimeSteppers::History history(2); + history.insert(TimeStepId(true, 0, slab.start()), decltype(history)::no_value, + db::prefix_variables(1, 3.0)); + + auto values = arbitrary_test_values(); + + SectorVariables variables = std::move(get(values)); + + const imex::solve_implicit_sector_detail::ImplicitEquation + equation(make_not_null(&variables), &history, TimeSteppers::Heun2{}, + slab.duration(), std::tuple<>{}, tmpl::type_{}); + const auto evolution_data = std::forward_as_tuple( + std::as_const(get(values)), std::as_const(get(values)), + std::as_const(get(values))); + imex::solve_implicit_sector_detail::ImplicitSolver + solver(equation, evolution_data); + solver.set_index(0); + std::array + initial_guess{}; + SectorVariables(initial_guess.data(), initial_guess.size()) = variables; + const auto jacobian = solver.jacobian(initial_guess); + + auto deriv_approx = Approx::custom().epsilon(1.0e-12); + // gsl_multiroot wants jacobian[i][j] = dfi/dxj + for (size_t j = 0; j < initial_guess.size(); ++j) { + const auto derivative = + numerical_derivative(solver, initial_guess, j, 1.0e-1); + for (size_t i = 0; i < initial_guess.size(); ++i) { + CHECK(gsl::at(gsl::at(jacobian, i), j) == deriv_approx(derivative[i])); + } + } +} + +template +void test_solve_implicit_sector(const imex::Mode solve_mode) { + using sector = ImplicitSector; + // No solve is done with an analytic solution. + const bool doing_semi_implicit_solve = + solve_mode == imex::Mode::SemiImplicit and not TestWithAnalyticSolution; + // We handle v1 entirely explicitly and v2, v3 entirely implicitly. + // The evolution equations for the latter two (coded in `Source` + // above) are + // d/dt[v2^ij] = v3^i v3^j - nt v2^ij + // d/dt[v3^i] = -v1 v3^i + + using variables_tag = Tags::Variables>; + using implicit_variables_source_tag = + Tags::Variables, ::Tags::Source>>; + using DtImplicitVariables = + Variables, ::Tags::dt>>; + using history_tag = imex::Tags::ImplicitHistory; + + const size_t number_of_grid_points = 5; + const Slab slab(3.0, 5.0); + const TimeStepId initial_time_step_id(true, 0, slab.start()); + const auto time_step = Slab(3.0, 5.0).duration() / 3; + + MAKE_GENERATOR(gen); + // Keep values positive to prevent the denominators in the analytic + // solution from becoming small. + std::uniform_real_distribution dist(0.0, 5.0); + const auto non_tensor = make_with_random_values(make_not_null(&gen), + make_not_null(&dist)); + const auto initial_vars = make_with_random_values( + make_not_null(&gen), make_not_null(&dist), number_of_grid_points); + auto box = db::create, Tags::TimeStep, history_tag, + imex::Tags::Mode, imex::Tags::SolveFailures, + imex::Tags::SolveTolerance>>( + initial_vars, non_tensor, VariablesFromEvolution::type{}, + std::make_unique(), time_step, + typename history_tag::type{2}, solve_mode, + Scalar(DataVector(number_of_grid_points, 0.0)), 1.0e-10); + + // Perform updates as if taking an explicit step. + const auto simulate_explicit_step = [&dist, &gen, &initial_vars]( + const auto box, + const TimeStepId& time_step_id) { + db::mutate( + [&dist, &gen, &initial_vars, &time_step_id]( + const gsl::not_null history, + const gsl::not_null var1, + const gsl::not_null var2, + const gsl::not_null var3, + const gsl::not_null test_variables, + const NonTensor::type& non_tensor) { + implicit_variables_source_tag::type source_vars(number_of_grid_points, + 0.0); + Source::work(&get>(source_vars), + &get>(source_vars), *var1, *var2, + *var3, non_tensor); + + history->insert( + time_step_id, history_tag::type::no_value, + source_vars + .reference_with_different_prefixes()); + // Update the explicitly evolved variable. + fill_with_random_values(var1, make_not_null(&gen), + make_not_null(&dist)); + // The explicit time derivative for var2 and var3 is + // zero, so the explicit integration will consider them + // constant and reset them to the initial value. + *var2 = get(initial_vars); + *var3 = get(initial_vars); + // This isn't evolved but we test obtaining it from the + // evolution box by checking for this value. + test_variables->initialize(get(*var1).size()); + get(get(*test_variables)) = 2.0 * get(*var1); + }, + box, db::get(*box)); + }; + + simulate_explicit_step(make_not_null(&box), initial_time_step_id); + + auto guess = initial_vars; + InitialGuess::apply(&get(guess), &get(guess), {}, {}); + + db::mutate_apply>( + make_not_null(&box)); + + const double dt = time_step.value(); + const auto final_vars = db::get(box); + Var2::type expected_var2{}; + Var3::type expected_var3{}; + Var2::type expected_var2_final{}; + if (not doing_semi_implicit_solve) { + // The first implicit substep for the Heun stepper is + // y(dt) = y(0) + dt/2 (d/d[y(0)] + d/dt[y(dt)]) + + // The analytic solution for the result of the first substep is + // v3^i(dt) = v3^i(0) (1 - dt/2 v1(0)) / (1 + dt/2 v1(dt)) + // v2^ij(dt) = (v2^ij(0) (1 - dt/2 nt) + + // + dt/2 (v3^i(0) v3^j(0) + v3^i(dt) v3^j(dt))) + // / (1 + dt/2 nt) + + tenex::evaluate(make_not_null(&expected_var3), + (1.0 - 0.5 * dt * get(initial_vars)()) / + (1.0 + 0.5 * dt * get(final_vars)()) * + get(initial_vars)(ti::I)); + tenex::evaluate( + make_not_null(&expected_var2), + ((1.0 - 0.5 * dt * non_tensor) * get(initial_vars)(ti::I, ti::J) + + 0.5 * dt * + (get(initial_vars)(ti::I) * get(initial_vars)(ti::J) + + expected_var3(ti::I) * expected_var3(ti::J))) / + (1.0 + 0.5 * dt * non_tensor)); + + // The second implicit substep is simpler, since it isn't actually + // implicit, and is in fact the same as the first substep. + expected_var2_final = expected_var2; + } else { + // For a semi-implicit method, we expand the source term around + // the initial guess: + // + // S(y(dt)) = S(guess + (y(dt) - guess)) + // ~= S(guess) + J(guess) (y(dt) - guess) + // + // For Heun's method + // + // y(dt) = y(0) + dt/2 (d/d[y(0)] + d/dt[y(dt)]) + // + // this gives the equation + // + // [1 - dt/2 J(guess)] (y(dt) - guess) + // = y(0) - guess + dt/2 (d/d[y(0)] + S(guess)) + // + // with + // + // S[v2^ij] = v3^i v3^j - nt v2^ij + // S[v3^i] = -v1 v3^i + // J[v2^ij] : dS[v2^ij]/dv2^kl = -nt delta^i_k delta^j_l + // dS[v2^ij]/dv3^k = v3^i delta^j_k + v3^j delta^i_k + // J[v3^i] : dS[v3^i]/dv2^jk = 0 + // dS[v3^i]/dv3^j = -v1 delta^i_j + // + // This is fairly easy to solve as (denoting the guess as g2 and g3): + // + // v3^i(dt) = (1 - dt/2 v1(0)) / (1 + dt/2 v1(dt)) v3^i(0) + // (solved exactly as the source is linear) + // + // v2^ij(dt) = + // (1 - dt/2 nt) / (1 + dt/2 nt) v2^ij(0) + // + dt/2 / (1 + dt/2 nt) [ + // v3^i(0) v3^j(0) + g3^i v3^j(dt) + v3^i(dt) g3^j - g3^i g3^j] + + tenex::evaluate(make_not_null(&expected_var3), + (1.0 - 0.5 * dt * get(initial_vars)()) / + (1.0 + 0.5 * dt * get(final_vars)()) * + get(initial_vars)(ti::I)); + tenex::evaluate( + make_not_null(&expected_var2), + ((1.0 - 0.5 * dt * non_tensor) * get(initial_vars)(ti::I, ti::J) + + 0.5 * dt * + (get(initial_vars)(ti::I) * get(initial_vars)(ti::J) + + get(guess)(ti::I) * expected_var3(ti::J) + + expected_var3(ti::I) * get(guess)(ti::J) - + get(guess)(ti::I) * get(guess)(ti::J))) / + (1.0 + 0.5 * dt * non_tensor)); + + // The second substep is explicit, so doesn't do a semi-implicit + // solve, but still uses the results from the first substep. Var3 + // is again exact. + tenex::evaluate( + make_not_null(&expected_var2_final), + (1.0 - 0.5 * dt * non_tensor) * get(initial_vars)(ti::I, ti::J) + + 0.5 * dt * + (get(initial_vars)(ti::I) * + get(initial_vars)(ti::J) + + expected_var3(ti::I) * expected_var3(ti::J) - + non_tensor * expected_var2(ti::I, ti::J))); + } + CHECK_ITERABLE_APPROX(get(final_vars), expected_var2); + CHECK_ITERABLE_APPROX(get(final_vars), expected_var3); + + CHECK(db::get(box).size() == 1); + CHECK(db::get(box).substeps().empty()); + + simulate_explicit_step(make_not_null(&box), + initial_time_step_id.next_substep(time_step, 1.0)); + performing_step_with_no_implicit_term = true; + db::mutate_apply>( + make_not_null(&box)); + performing_step_with_no_implicit_term = false; + CHECK_ITERABLE_APPROX(get(db::get(box)), + expected_var2_final); + CHECK_ITERABLE_APPROX(get(db::get(box)), expected_var3); + + CHECK(db::get(box).size() == 1); + CHECK(db::get(box).substeps().size() == 1); + + // Take another substep just to test the history cleanup. + simulate_explicit_step(make_not_null(&box), + initial_time_step_id.next_step(time_step)); + db::mutate_apply>( + make_not_null(&box)); + + CHECK(db::get(box).size() == 1); + CHECK(db::get(box).substeps().empty()); +} + +struct ResettingTestSector : tt::ConformsTo { + using tensors = tmpl::list; + using initial_guess = imex::GuessExplicitResult; + + struct SolveAttempt { + using tags_from_evolution = tmpl::list; + using simple_tags = tmpl::list<>; + using compute_tags = tmpl::list<>; + + using source_prep = tmpl::list<>; + using jacobian_prep = tmpl::list<>; + + struct source { + using return_tags = tmpl::list<::Tags::Source>; + using argument_tags = tmpl::list; + + static void apply(const gsl::not_null*> source_var1, + const tnsr::II& var2) { + get(*source_var1) = get<0, 0>(var2); + } + }; + + struct jacobian { + using return_tags = tmpl::list<>; + using argument_tags = tmpl::list<>; + + static void apply() {} + }; + }; + + using solve_attempts = tmpl::list; +}; + +// There was a bug where internal cached values did not clear properly +// between points. +void test_point_reseting() { + using variables_tag = ::Tags::Variables>; + + const Slab slab(0.0, 2.0); + const auto time_step = slab.duration(); + + auto var2 = make_with_value>(2_st, 0.0); + get<0, 0>(var2) = DataVector{1.0, 2.0}; + + // Set the initial value and derivative to zero so we can ignore + // those terms in the time stepper equation. + variables_tag::type initial_value(2, 0.0); // NOLINT(misc-const-correctness) + TimeSteppers::History history(2); + history.insert(TimeStepId(true, 0, slab.start()), decltype(history)::no_value, + db::prefix_variables(2, 0.0)); + + auto box = db::create, + imex::Tags::Mode, Tags::TimeStepper, Tags::TimeStep, + imex::Tags::SolveFailures, + imex::Tags::SolveTolerance>>( + var2, std::move(initial_value), std::move(history), + imex::Mode::SemiImplicit, std::make_unique(), + time_step, Scalar(DataVector(2, 0.0)), 1.0e-10); + db::mutate_apply< + imex::SolveImplicitSector>( + make_not_null(&box)); + + // The equation being solved is: y(dt) = dt/2 source(y(dt)) + // where: dt = 2, source(y) = var2(0, 0) + CHECK_ITERABLE_APPROX(get(get(box)), (get<0, 0>(var2))); +} + +struct DesiredLevel : db::SimpleTag { + using type = Scalar; +}; + +struct SectorWithFallback : tt::ConformsTo { + using tensors = tmpl::list; + using initial_guess = imex::GuessExplicitResult; + + template + struct SolveAttempt { + using tags_from_evolution = tmpl::list; + using simple_tags = tmpl::list<>; + using compute_tags = tmpl::list<>; + + using source_prep = tmpl::list<>; + using jacobian_prep = tmpl::list<>; + + struct source { + using return_tags = tmpl::list<::Tags::Source>; + using argument_tags = tmpl::list; + + static void apply(const gsl::not_null*> source_var1, + const Scalar& desired_level) { + CHECK(get(desired_level)[0] <= Level); + get(*source_var1) = Level; + } + }; + + struct jacobian { + using return_tags = + tmpl::list>>; + using argument_tags = tmpl::list; + + static void apply(const gsl::not_null*> dvar1_dvar1, + const Scalar& desired_level) { + if (get(desired_level)[0] == Level) { + get(*dvar1_dvar1) = 0.0; + } else { + // Cause a failure from a non-invertable y = y + constant + get(*dvar1_dvar1) = 1.0; + } + } + }; + }; + + using solve_attempts = + tmpl::list, SolveAttempt<3>, SolveAttempt<2>, + SolveAttempt<1>, SolveAttempt<0>>; +}; + +void test_fallback() { + using sector = SectorWithFallback; + using variables_tag = ::Tags::Variables>; + using history_tag = imex::Tags::ImplicitHistory; + + const Slab slab(0.0, 2.0); + const auto time_step = slab.duration(); + const Scalar desired_level{{{{1.0, 3.0, 4.0, 4.0, 3.0}}}}; + const size_t number_of_grid_points = get(desired_level).size(); + + // Set the initial value and derivative to zero so we can ignore + // those terms in the time stepper equation. + // NOLINTNEXTLINE(misc-const-correctness) + variables_tag::type initial_value(number_of_grid_points, 0.0); + TimeSteppers::History history(2); + history.insert(TimeStepId(true, 0, slab.start()), decltype(history)::no_value, + db::prefix_variables( + number_of_grid_points, 0.0)); + + auto box = db::create, Tags::TimeStep, + imex::Tags::SolveFailures, imex::Tags::SolveTolerance>>( + desired_level, std::move(initial_value), std::move(history), + imex::Mode::SemiImplicit, std::make_unique(), + time_step, Scalar(DataVector(number_of_grid_points, 0.0)), + 1.0e-10); + + // Jacobian is only right when the template parameter and DataVector + // agree. (This is also a test of test_sector will fallbacks.) + TestHelpers::imex::test_sector>( + 1.0e-1, 1.0e-14, db::get(box), + {Scalar{{{{2.0}}}}}); + + db::mutate_apply>( + make_not_null(&box)); + + // The equation being solved is: y(dt) = dt/2 source(y(dt)) + // where: dt = 2, source(y) = desired_level + CHECK_ITERABLE_APPROX(get(get(box)), get(desired_level)); + CHECK(get(get>(box)) == + 4.0 - get(desired_level)); +} +} // namespace + +SPECTRE_TEST_CASE("Unit.Evolution.Imex.SolveImplicitSector", + "[Unit][Evolution]") { + test_test_sector(); + test_test_sector(); + test_internal_jacobian_ordering(); + test_solve_implicit_sector(imex::Mode::Implicit); + test_solve_implicit_sector(imex::Mode::Implicit); + test_solve_implicit_sector(imex::Mode::SemiImplicit); + test_solve_implicit_sector(imex::Mode::SemiImplicit); + test_point_reseting(); + test_fallback(); +} diff --git a/tests/Unit/Helpers/Evolution/CMakeLists.txt b/tests/Unit/Helpers/Evolution/CMakeLists.txt index 059860cb365db..fa8d156f74d60 100644 --- a/tests/Unit/Helpers/Evolution/CMakeLists.txt +++ b/tests/Unit/Helpers/Evolution/CMakeLists.txt @@ -4,4 +4,5 @@ add_subdirectory(EventsAndDenseTriggers) add_subdirectory(DgSubcell) add_subdirectory(DiscontinuousGalerkin) +add_subdirectory(Imex) add_subdirectory(Systems) diff --git a/tests/Unit/Helpers/Evolution/Imex/CMakeLists.txt b/tests/Unit/Helpers/Evolution/Imex/CMakeLists.txt new file mode 100644 index 0000000000000..5df5d9313a02a --- /dev/null +++ b/tests/Unit/Helpers/Evolution/Imex/CMakeLists.txt @@ -0,0 +1,21 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +set(LIBRARY "ImexHelpers") + +add_spectre_library(${LIBRARY} INTERFACE) + +spectre_target_headers( + ${LIBRARY} + INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/tests/Unit + HEADERS + TestSector.hpp + ) + +target_link_libraries( + ${LIBRARY} + INTERFACE + DataStructures + Imex + Utilities + ) diff --git a/tests/Unit/Helpers/Evolution/Imex/TestSector.hpp b/tests/Unit/Helpers/Evolution/Imex/TestSector.hpp new file mode 100644 index 0000000000000..64140a0c1b556 --- /dev/null +++ b/tests/Unit/Helpers/Evolution/Imex/TestSector.hpp @@ -0,0 +1,209 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include "Framework/TestingFramework.hpp" + +#include +#include +#include +#include + +#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataBox/PrefixHelpers.hpp" +#include "DataStructures/DataBox/Prefixes.hpp" +#include "DataStructures/DataBox/Tag.hpp" +#include "DataStructures/Tensor/ContractFirstNIndices.hpp" +#include "DataStructures/Variables.hpp" +#include "DataStructures/VariablesTag.hpp" +#include "Evolution/Imex/GuessResult.hpp" +#include "Evolution/Imex/Protocols/ImplicitSector.hpp" +#include "Evolution/Imex/Tags/Jacobian.hpp" +#include "Framework/TestHelpers.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/Literals.hpp" +#include "Utilities/PrettyType.hpp" +#include "Utilities/ProtocolHelpers.hpp" +#include "Utilities/TMPL.hpp" +#include "Utilities/TaggedTuple.hpp" +#include "Utilities/TypeTraits/IsA.hpp" + +namespace TestHelpers::imex { +namespace test_sector_detail { +template +struct ReadOnlySource : ::db::SimpleTag { + using type = typename Tag::type; +}; + +template +struct ReadOnly : Tag, ::db::ComputeTag { + using base = Tag; + using argument_tags = tmpl::list>; + static void function(gsl::not_null dest, + const typename ReadOnly::type& source) { + *dest = source; + } +}; + +template +struct OnlyEntry; + +template +struct OnlyEntry> { + using type = T; +}; +} // namespace test_sector_detail + +/// Check that an implicit sector conforms to +/// ::imex::protocols::ImplicitSector and that its `source` and +/// `jacobian` are consistent with each other. +template ::type> +void test_sector(const double stencil_size, const double tolerance, + const Variables& explicit_values, + tuples::tagged_tuple_from_typelist< + typename SolveAttempt::tags_from_evolution> + evolution_data) { + static_assert( + tt::assert_conforms_to_v); + static_assert( + tmpl::list_contains_v); + + using sector_variables_tag = ::Tags::Variables; + using SectorVariables = typename sector_variables_tag::type; + using sector_source_tag = + ::db::add_tag_prefix<::Tags::Source, sector_variables_tag>; + + using sector_jacobian_tag = ::Tags::Variables<::imex::jacobian_tags< + typename Sector::tensors, typename sector_source_tag::type::tags_list>>; + + // Make tags that are immutable in real use immutable here too. + using read_only_tags_from_evolution_source = + tmpl::transform>; + using read_only_tags_from_evolution = + tmpl::transform>; + + using simple_tags = tmpl::append< + read_only_tags_from_evolution_source, typename SolveAttempt::simple_tags, + tmpl::list>; + using compute_tags = tmpl::append; + auto box = ::db::create(); + tmpl::for_each([&](auto tag) { + using Tag = tmpl::type_from; + ::db::mutate>( + [&](const gsl::not_null box_value) { + *box_value = std::move(get(evolution_data)); + }, + make_not_null(&box)); + }); + using variables_tags = tmpl::push_front< + tmpl::filter>>, + sector_source_tag, sector_jacobian_tag>; + ::db::mutate_apply>( + [](const auto... vars) { expand_pack((vars->initialize(1, 0.0), 0)...); }, + make_not_null(&box)); + ::db::mutate( + [&](const gsl::not_null vars) { + *vars = explicit_values; + }, + make_not_null(&box)); + + // Arbitrary values. These would depend on the choice of time + // stepper and step size and such. + const SectorVariables inhomogeneous_terms = 1.1 * explicit_values; + const double implicit_weight = 0.4; + + const std::vector<::imex::GuessResult> guess_type = + ::db::mutate_apply( + make_not_null(&box), inhomogeneous_terms, implicit_weight); + CHECK(guess_type.size() <= 1); + const SectorVariables initial_guess = ::db::get(box); + if (guess_type == std::vector{::imex::GuessResult::ExactSolution}) { + // If the sector can always be solved analytically, the jacobian + // need not be coded. Instead, we check that the supplied + // solution is correct. + tmpl::for_each([&](auto mutator) { + using Mutator = tmpl::type_from; + ::db::mutate_apply(make_not_null(&box)); + }); + ::db::mutate_apply(make_not_null(&box)); + const auto& source = ::db::get(box); + const SectorVariables step_result = + inhomogeneous_terms + implicit_weight * source; + CHECK_VARIABLES_CUSTOM_APPROX(initial_guess, step_result, + Approx::custom().epsilon(tolerance)); + } else { + // Check the jacobian numerically. + tmpl::for_each([&](auto mutator) { + using Mutator = tmpl::type_from; + ::db::mutate_apply(make_not_null(&box)); + }); + ::db::mutate_apply(make_not_null(&box)); + const auto jacobian = ::db::get(box); + tmpl::for_each([&](auto varying_tensor_v) { + using varying_tensor = tmpl::type_from; + CAPTURE(pretty_type::get_name()); + for (size_t varying_component = 0; + varying_component < varying_tensor::type::size(); + ++varying_component) { + ::db::mutate( + [&](const gsl::not_null vars) { + *vars = initial_guess; + }, + make_not_null(&box)); + const auto finite_difference_derivative = numerical_derivative( + [&](const std::array& component_value) { + ::db::mutate( + [&](const gsl::not_null var) { + (*var)[varying_component] = component_value[0]; + }, + make_not_null(&box)); + tmpl::for_each( + [&](auto mutator) { + using Mutator = tmpl::type_from; + ::db::mutate_apply(make_not_null(&box)); + }); + ::db::mutate_apply( + make_not_null(&box)); + return ::db::get(box); + }, + std::array{ + get(initial_guess)[varying_component][0]}, + 0, stencil_size); + typename varying_tensor::type variation(1_st, 0.0); + variation[varying_component] = 1.0; + CAPTURE(variation); + tmpl::for_each( + [&](auto source_tensor_v) { + using source_tensor = tmpl::type_from; + CAPTURE(pretty_type::get_name()); + auto analytic_derivative = contract_first_n_indices< + varying_tensor::type::num_tensor_indices>( + variation, + get<::imex::Tags::Jacobian>( + jacobian)); + // We want the derivative with respect to a single + // component, but with multiplicities all the equivalent + // components were actually set, so the result is larger + // than it should be. + for (auto& component : analytic_derivative) { + component /= variation.multiplicity(varying_component); + } + static_assert( + std::is_same_v, + typename source_tensor::type>); + CHECK_ITERABLE_CUSTOM_APPROX( + analytic_derivative, + get(finite_difference_derivative), + Approx::custom().epsilon(tolerance)); + }); + } + }); + } +} +} // namespace TestHelpers::imex