Merge pull request #5653 from wthrowe/imex/solve_implicit_sector

Add IMEX SolveImplicitSector mutator

.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

src/DataStructures/Tensor/TypeAliases.hpp

@@ -399,6 +399,12 @@ using iJkk = Tensor<DataType, tmpl::integral_list<std::int32_t, 3, 2, 1, 1>,
                                SpatialIndex<SpatialDim, UpLo::Lo, Fr>,
                                SpatialIndex<SpatialDim, UpLo::Lo, Fr>>>;
 
+template <typename DataType, size_t SpatialDim, typename Fr = Frame::Inertial>
+using iiJJ = Tensor<DataType, tmpl::integral_list<std::int32_t, 2, 2, 1, 1>,
+                    index_list<SpatialIndex<SpatialDim, UpLo::Lo, Fr>,
+                               SpatialIndex<SpatialDim, UpLo::Lo, Fr>,
+                               SpatialIndex<SpatialDim, UpLo::Up, Fr>,
+                               SpatialIndex<SpatialDim, UpLo::Up, Fr>>>;
 }  // namespace tnsr
 
 /*!

src/Evolution/Imex/CMakeLists.txt

@@ -9,22 +9,32 @@ spectre_target_sources(
   ${LIBRARY}
   PRIVATE
   GuessResult.cpp
+  Mode.cpp
   )
 
 spectre_target_headers(
   ${LIBRARY}
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
   GuessResult.hpp
+  Mode.hpp
   NamespaceDocs.hpp
+  SolveImplicitSector.hpp
+  SolveImplicitSector.tpp
   )
 
 target_link_libraries(
   ${LIBRARY}
   PUBLIC
   ErrorHandling
+  Options
   INTERFACE
+  DataStructures
+  LinearSolver
+  RootFinding
+  Time
   Utilities
   )
 
 add_subdirectory(Protocols)
+add_subdirectory(Tags)

src/Evolution/Imex/Mode.cpp

@@ -0,0 +1,22 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Imex/Mode.hpp"
+
+#include <string>
+
+#include "Options/ParseOptions.hpp"
+
+template <>
+imex::Mode Options::create_from_yaml<imex::Mode>::create<void>(
+    const Options::Option& options) {
+  const auto mode = options.parse_as<std::string>();
+  if (mode == "Implicit") {
+    return imex::Mode::Implicit;
+  } else if (mode == "SemiImplicit") {
+    return imex::Mode::SemiImplicit;
+  } else {
+    PARSE_ERROR(options.context(),
+                "Invalid IMEX mode.  Must be Implicit or SemiImplicit.");
+  }
+}

src/Evolution/Imex/Mode.hpp

@@ -0,0 +1,28 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include "Options/Options.hpp"
+
+namespace imex {
+/// IMEX implementations
+enum class Mode {
+  /// Solve the implicit equation using a nonlinear solver.
+  Implicit,
+  /// Solve a linearized version of the implicit equation.
+  SemiImplicit,
+};
+}  // namespace imex
+
+template <>
+struct Options::create_from_yaml<imex::Mode> {
+  template <typename Metavariables>
+  static imex::Mode create(const Options::Option& options) {
+    return create<void>(options);
+  }
+};
+
+template <>
+imex::Mode Options::create_from_yaml<imex::Mode>::create<void>(
+    const Options::Option& options);

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 <typename ConformingType>
   struct test {

src/Evolution/Imex/SolveImplicitSector.hpp

@@ -0,0 +1,122 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <tuple>
+#include <type_traits>
+
+#include "DataStructures/Tensor/TypeAliases.hpp"
+#include "Evolution/Imex/Mode.hpp"
+#include "Evolution/Imex/Protocols/ImplicitSector.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/ProtocolHelpers.hpp"
+#include "Utilities/TMPL.hpp"
+
+/// \cond
+class DataVector;
+class ImexTimeStepper;
+class TimeDelta;
+template <typename TagsList>
+class Variables;
+namespace Tags {
+struct TimeStep;
+template <typename StepperInterface>
+struct TimeStepper;
+}  // namespace Tags
+namespace TimeSteppers {
+template <typename Vars>
+class History;
+}  // namespace TimeSteppers
+namespace imex::Tags {
+template <typename ImplicitSector>
+struct ImplicitHistory;
+struct Mode;
+template <typename Sector>
+struct SolveFailures;
+struct SolveTolerance;
+}  // namespace imex::Tags
+/// \endcond
+
+namespace imex {
+namespace solve_implicit_sector_detail {
+template <typename Tags>
+using ForwardTuple = tmpl::wrap<
+    tmpl::transform<Tags, std::add_lvalue_reference<std::add_const<
+                              tmpl::bind<tmpl::type_from, tmpl::_1>>>>,
+    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 <typename SystemVariablesTag, typename ImplicitSector>
+struct SolveImplicitSector {
+  static_assert(
+      tt::assert_conforms_to_v<ImplicitSector, protocols::ImplicitSector>);
+
+ public:
+  using SystemVariables = typename SystemVariablesTag::type;
+  using SectorVariables = Variables<typename ImplicitSector::tensors>;
+
+ private:
+  template <typename Attempt>
+  struct get_tags_from_evolution {
+    using type = typename Attempt::tags_from_evolution;
+  };
+
+  using tags_for_each_attempt =
+      tmpl::transform<typename ImplicitSector::solve_attempts,
+                      get_tags_from_evolution<tmpl::_1>>;
+  // List of tags used for the initial guess followed by lists of tags
+  // used for each solve attempt.
+  using evolution_data_tags =
+      tmpl::push_front<tags_for_each_attempt,
+                       typename ImplicitSector::initial_guess::argument_tags>;
+
+  using EvolutionDataTuple = solve_implicit_sector_detail::ForwardTuple<
+      tmpl::join<evolution_data_tags>>;
+
+  static void apply_impl(
+      gsl::not_null<SystemVariables*> system_variables,
+      gsl::not_null<TimeSteppers::History<SectorVariables>*> implicit_history,
+      gsl::not_null<Scalar<DataVector>*> 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<SystemVariablesTag,
+                                 imex::Tags::ImplicitHistory<ImplicitSector>,
+                                 Tags::SolveFailures<ImplicitSector>>;
+  using argument_tags = tmpl::append<
+      tmpl::list<::Tags::TimeStepper<ImexTimeStepper>, ::Tags::TimeStep,
+                 Tags::Mode, Tags::SolveTolerance>,
+      tmpl::join<evolution_data_tags>>;
+
+  template <typename... ForwardArgs>
+  static void apply(const gsl::not_null<SystemVariables*> system_variables,
+                    const gsl::not_null<TimeSteppers::History<SectorVariables>*>
+                        implicit_history,
+                    const gsl::not_null<Scalar<DataVector>*> 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