SpECTRE Documentation Coverage Report
Current view: top level - Elliptic/Executables/Xcts - SolveXcts.hpp Hit Total Coverage
Commit: b5ffa4904430ccef0b226f73dcd25c74cb5188f6 Lines: 0 20 0.0 %
Date: 2021-07-28 22:05:18
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : 
       8             : #include "DataStructures/DataBox/PrefixHelpers.hpp"
       9             : #include "Domain/Creators/Factory3D.hpp"
      10             : #include "Domain/Creators/RegisterDerivedWithCharm.hpp"
      11             : #include "Domain/Tags.hpp"
      12             : #include "Elliptic/Actions/InitializeAnalyticSolution.hpp"
      13             : #include "Elliptic/Actions/InitializeFields.hpp"
      14             : #include "Elliptic/Actions/InitializeFixedSources.hpp"
      15             : #include "Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp"
      16             : #include "Elliptic/DiscontinuousGalerkin/Actions/InitializeDomain.hpp"
      17             : #include "Elliptic/DiscontinuousGalerkin/DgElementArray.hpp"
      18             : #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/InitializeSubdomain.hpp"
      19             : #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp"
      20             : #include "Elliptic/Systems/Xcts/FirstOrderSystem.hpp"
      21             : #include "Elliptic/Tags.hpp"
      22             : #include "Elliptic/Triggers/EveryNIterations.hpp"
      23             : #include "IO/Observer/Actions/RegisterEvents.hpp"
      24             : #include "IO/Observer/Helpers.hpp"
      25             : #include "IO/Observer/ObserverComponent.hpp"
      26             : #include "NumericalAlgorithms/Convergence/Tags.hpp"
      27             : #include "Options/Options.hpp"
      28             : #include "Options/Protocols/FactoryCreation.hpp"
      29             : #include "Parallel/Actions/SetupDataBox.hpp"
      30             : #include "Parallel/Actions/TerminatePhase.hpp"
      31             : #include "Parallel/GlobalCache.hpp"
      32             : #include "Parallel/InitializationFunctions.hpp"
      33             : #include "Parallel/PhaseDependentActionList.hpp"
      34             : #include "Parallel/Reduction.hpp"
      35             : #include "Parallel/RegisterDerivedClassesWithCharm.hpp"
      36             : #include "ParallelAlgorithms/Events/Factory.hpp"
      37             : #include "ParallelAlgorithms/EventsAndTriggers/Actions/RunEventsAndTriggers.hpp"
      38             : #include "ParallelAlgorithms/EventsAndTriggers/Completion.hpp"
      39             : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
      40             : #include "ParallelAlgorithms/EventsAndTriggers/LogicalTriggers.hpp"
      41             : #include "ParallelAlgorithms/EventsAndTriggers/Tags.hpp"
      42             : #include "ParallelAlgorithms/EventsAndTriggers/Trigger.hpp"
      43             : #include "ParallelAlgorithms/Initialization/Actions/AddComputeTags.hpp"
      44             : #include "ParallelAlgorithms/Initialization/Actions/RemoveOptionsAndTerminatePhase.hpp"
      45             : #include "ParallelAlgorithms/LinearSolver/Actions/MakeIdentityIfSkipped.hpp"
      46             : #include "ParallelAlgorithms/LinearSolver/Gmres/Gmres.hpp"
      47             : #include "ParallelAlgorithms/LinearSolver/Schwarz/Actions/ResetSubdomainSolver.hpp"
      48             : #include "ParallelAlgorithms/LinearSolver/Schwarz/Schwarz.hpp"
      49             : #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
      50             : #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/NewtonRaphson.hpp"
      51             : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
      52             : #include "PointwiseFunctions/AnalyticData/Xcts/AnalyticData.hpp"
      53             : #include "PointwiseFunctions/AnalyticData/Xcts/Binary.hpp"
      54             : #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
      55             : #include "PointwiseFunctions/AnalyticSolutions/Xcts/AnalyticSolution.hpp"
      56             : #include "PointwiseFunctions/AnalyticSolutions/Xcts/Flatness.hpp"
      57             : #include "PointwiseFunctions/AnalyticSolutions/Xcts/Kerr.hpp"
      58             : #include "PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.hpp"
      59             : #include "Utilities/Blas.hpp"
      60             : #include "Utilities/ErrorHandling/FloatingPointExceptions.hpp"
      61             : #include "Utilities/Functional.hpp"
      62             : #include "Utilities/MemoryHelpers.hpp"
      63             : #include "Utilities/ProtocolHelpers.hpp"
      64             : #include "Utilities/TMPL.hpp"
      65             : 
      66           0 : namespace SolveXcts::OptionTags {
      67           0 : struct NonlinearSolverGroup {
      68           0 :   static std::string name() noexcept { return "NonlinearSolver"; }
      69           0 :   static constexpr Options::String help = "The iterative nonlinear solver";
      70             : };
      71           0 : struct NewtonRaphsonGroup {
      72           0 :   static std::string name() noexcept { return "NewtonRaphson"; }
      73           0 :   static constexpr Options::String help =
      74             :       "Options for the Newton-Raphson nonlinear solver";
      75           0 :   using group = NonlinearSolverGroup;
      76             : };
      77           0 : struct LinearSolverGroup {
      78           0 :   static std::string name() noexcept { return "LinearSolver"; }
      79           0 :   static constexpr Options::String help =
      80             :       "The iterative Krylov-subspace linear solver";
      81             : };
      82           0 : struct GmresGroup {
      83           0 :   static std::string name() noexcept { return "Gmres"; }
      84           0 :   static constexpr Options::String help = "Options for the GMRES linear solver";
      85           0 :   using group = LinearSolverGroup;
      86             : };
      87           0 : struct SchwarzSmootherGroup {
      88           0 :   static std::string name() noexcept { return "SchwarzSmoother"; }
      89           0 :   static constexpr Options::String help = "Options for the Schwarz smoother";
      90           0 :   using group = LinearSolverGroup;
      91             : };
      92             : }  // namespace SolveXcts::OptionTags
      93             : 
      94             : /// \cond
      95             : struct Metavariables {
      96             :   static constexpr size_t volume_dim = 3;
      97             :   static constexpr int conformal_matter_scale = 0;
      98             :   using system =
      99             :       Xcts::FirstOrderSystem<Xcts::Equations::HamiltonianLapseAndShift,
     100             :                              Xcts::Geometry::Curved, conformal_matter_scale>;
     101             : 
     102             :   // List the possible backgrounds, i.e. the variable-independent part of the
     103             :   // equations that define the problem to solve (along with the boundary
     104             :   // conditions)
     105             :   using analytic_solution_registrars =
     106             :       tmpl::list<Xcts::Solutions::Registrars::Schwarzschild,
     107             :                  Xcts::Solutions::Registrars::Kerr>;
     108             :   using analytic_data_registrars = tmpl::list<
     109             :       Xcts::AnalyticData::Registrars::Binary<analytic_solution_registrars>>;
     110             :   using background_tag = elliptic::Tags::Background<
     111             :       ::AnalyticData<volume_dim, tmpl::append<analytic_solution_registrars,
     112             :                                               analytic_data_registrars>>>;
     113             : 
     114             :   // List the possible initial guesses
     115             :   using initial_guess_registrars =
     116             :       tmpl::append<tmpl::list<Xcts::Solutions::Registrars::Flatness>,
     117             :                    analytic_solution_registrars, analytic_data_registrars>;
     118             :   using initial_guess_tag = elliptic::Tags::InitialGuess<
     119             :       ::AnalyticData<volume_dim, initial_guess_registrars>>;
     120             : 
     121             :   static constexpr Options::String help{
     122             :       "Find the solution to an XCTS problem."};
     123             : 
     124             :   // These are the fields we solve for
     125             :   using fields_tag = ::Tags::Variables<typename system::primal_fields>;
     126             :   // These are the fluxes corresponding to the fields, i.e. essentially their
     127             :   // first derivatives. These are background fields for the linearized sources.
     128             :   using fluxes_tag = ::Tags::Variables<typename system::primal_fluxes>;
     129             :   // These are the fixed sources, i.e. the RHS of the equations
     130             :   using fixed_sources_tag = db::add_tag_prefix<::Tags::FixedSource, fields_tag>;
     131             :   using operator_applied_to_fields_tag =
     132             :       db::add_tag_prefix<NonlinearSolver::Tags::OperatorAppliedTo, fields_tag>;
     133             : 
     134             :   using nonlinear_solver = NonlinearSolver::newton_raphson::NewtonRaphson<
     135             :       Metavariables, fields_tag, SolveXcts::OptionTags::NewtonRaphsonGroup,
     136             :       fixed_sources_tag>;
     137             :   using nonlinear_solver_iteration_id =
     138             :       Convergence::Tags::IterationId<typename nonlinear_solver::options_group>;
     139             : 
     140             :   // The linear solver algorithm. We must use GMRES since the operator is
     141             :   // not guaranteed to be symmetric.
     142             :   using linear_solver = LinearSolver::gmres::Gmres<
     143             :       Metavariables, typename nonlinear_solver::linear_solver_fields_tag,
     144             :       SolveXcts::OptionTags::GmresGroup, true,
     145             :       typename nonlinear_solver::linear_solver_source_tag>;
     146             :   using linear_solver_iteration_id =
     147             :       Convergence::Tags::IterationId<typename linear_solver::options_group>;
     148             :   // Precondition each linear solver iteration with a number of Schwarz
     149             :   // smoothing steps
     150             :   using subdomain_operator =
     151             :       elliptic::dg::subdomain_operator::SubdomainOperator<
     152             :           system, SolveXcts::OptionTags::SchwarzSmootherGroup>;
     153             :   // This data needs to be communicated on subdomain overlap regions
     154             :   using communicated_overlap_tags = tmpl::list<
     155             :       // For linearized sources
     156             :       fields_tag, fluxes_tag>;
     157             :   using schwarz_smoother = LinearSolver::Schwarz::Schwarz<
     158             :       typename linear_solver::operand_tag,
     159             :       SolveXcts::OptionTags::SchwarzSmootherGroup, subdomain_operator,
     160             :       typename linear_solver::preconditioner_source_tag>;
     161             :   // For the GMRES linear solver we need to apply the DG operator to its
     162             :   // internal "operand" in every iteration of the algorithm.
     163             :   using correction_vars_tag = typename linear_solver::operand_tag;
     164             :   using operator_applied_to_correction_vars_tag =
     165             :       db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo,
     166             :                          correction_vars_tag>;
     167             :   // The correction fluxes can be stored in an arbitrary tag. We don't need to
     168             :   // access them anywhere, they're just a memory buffer for the linearized
     169             :   // operator.
     170             :   using correction_fluxes_tag =
     171             :       db::add_tag_prefix<NonlinearSolver::Tags::Correction, fluxes_tag>;
     172             : 
     173             :   using analytic_solution_fields = tmpl::append<typename system::primal_fields,
     174             :                                                 typename system::primal_fluxes>;
     175             :   using observe_fields = tmpl::append<analytic_solution_fields,
     176             :                                       typename system::background_fields>;
     177             : 
     178             :   // Collect all items to store in the cache.
     179             :   using const_global_cache_tags =
     180             :       tmpl::list<background_tag, initial_guess_tag, Tags::EventsAndTriggers>;
     181             : 
     182             :   struct factory_creation
     183             :       : tt::ConformsTo<Options::protocols::FactoryCreation> {
     184             :     using factory_classes = tmpl::map<
     185             :         tmpl::pair<DomainCreator<volume_dim>, domain_creators<volume_dim>>,
     186             :         tmpl::pair<Event, tmpl::flatten<tmpl::list<
     187             :                               Events::Completion,
     188             :                               dg::Events::field_observations<
     189             :                                   volume_dim, nonlinear_solver_iteration_id,
     190             :                                   observe_fields, analytic_solution_fields>>>>,
     191             :         tmpl::pair<Trigger,
     192             :                    tmpl::push_back<Triggers::logical_triggers,
     193             :                                    elliptic::Triggers::EveryNIterations<
     194             :                                        nonlinear_solver_iteration_id>>>>;
     195             :   };
     196             : 
     197             :   // Collect all reduction tags for observers
     198             :   using observed_reduction_data_tags =
     199             :       observers::collect_reduction_data_tags<tmpl::flatten<
     200             :           tmpl::list<tmpl::at<factory_creation::factory_classes, Event>,
     201             :                      nonlinear_solver, linear_solver, schwarz_smoother>>>;
     202             : 
     203             :   // Specify all global synchronization points.
     204             :   enum class Phase { Initialization, RegisterWithObserver, Solve, Exit };
     205             : 
     206             :   using initialization_actions = tmpl::list<
     207             :       Actions::SetupDataBox,
     208             :       elliptic::dg::Actions::InitializeDomain<volume_dim>,
     209             :       typename nonlinear_solver::initialize_element,
     210             :       typename linear_solver::initialize_element,
     211             :       typename schwarz_smoother::initialize_element,
     212             :       elliptic::Actions::InitializeFields<system, initial_guess_tag>,
     213             :       elliptic::Actions::InitializeFixedSources<system, background_tag>,
     214             :       elliptic::Actions::InitializeOptionalAnalyticSolution<
     215             :           background_tag, analytic_solution_fields,
     216             :           Xcts::Solutions::AnalyticSolution<tmpl::append<
     217             :               analytic_solution_registrars, analytic_data_registrars>>>,
     218             :       elliptic::dg::Actions::initialize_operator<system, background_tag>,
     219             :       elliptic::dg::subdomain_operator::Actions::InitializeSubdomain<
     220             :           system, background_tag, typename schwarz_smoother::options_group>,
     221             :       Initialization::Actions::RemoveOptionsAndTerminatePhase>;
     222             : 
     223             :   template <bool Linearized>
     224             :   using build_operator_actions = elliptic::dg::Actions::apply_operator<
     225             :       system, Linearized,
     226             :       tmpl::conditional_t<Linearized, linear_solver_iteration_id,
     227             :                           nonlinear_solver_iteration_id>,
     228             :       tmpl::conditional_t<Linearized, correction_vars_tag, fields_tag>,
     229             :       tmpl::conditional_t<Linearized, correction_fluxes_tag, fluxes_tag>,
     230             :       tmpl::conditional_t<Linearized, operator_applied_to_correction_vars_tag,
     231             :                           operator_applied_to_fields_tag>>;
     232             : 
     233             :   using register_actions =
     234             :       tmpl::list<observers::Actions::RegisterEventsWithObservers,
     235             :                  typename nonlinear_solver::register_element,
     236             :                  typename schwarz_smoother::register_element,
     237             :                  Parallel::Actions::TerminatePhase>;
     238             : 
     239             :   using solve_actions = tmpl::list<
     240             :       typename nonlinear_solver::template solve<
     241             :           build_operator_actions<false>,
     242             :           tmpl::list<LinearSolver::Schwarz::Actions::SendOverlapFields<
     243             :                          communicated_overlap_tags,
     244             :                          typename schwarz_smoother::options_group, false>,
     245             :                      LinearSolver::Schwarz::Actions::ReceiveOverlapFields<
     246             :                          volume_dim, communicated_overlap_tags,
     247             :                          typename schwarz_smoother::options_group>,
     248             :                      LinearSolver::Schwarz::Actions::ResetSubdomainSolver<
     249             :                          typename schwarz_smoother::options_group>,
     250             :                      typename linear_solver::template solve<tmpl::list<
     251             :                          typename schwarz_smoother::template solve<
     252             :                              build_operator_actions<true>>,
     253             :                          ::LinearSolver::Actions::make_identity_if_skipped<
     254             :                              schwarz_smoother, build_operator_actions<true>>>>>,
     255             :           Actions::RunEventsAndTriggers>,
     256             :       Parallel::Actions::TerminatePhase>;
     257             : 
     258             :   using dg_element_array = elliptic::DgElementArray<
     259             :       Metavariables,
     260             :       tmpl::list<Parallel::PhaseActions<Phase, Phase::Initialization,
     261             :                                         initialization_actions>,
     262             :                  Parallel::PhaseActions<Phase, Phase::RegisterWithObserver,
     263             :                                         register_actions>,
     264             :                  Parallel::PhaseActions<Phase, Phase::Solve, solve_actions>>>;
     265             : 
     266             :   // Specify all parallel components that will execute actions at some point.
     267             :   using component_list = tmpl::flatten<
     268             :       tmpl::list<dg_element_array, typename nonlinear_solver::component_list,
     269             :                  typename linear_solver::component_list,
     270             :                  typename schwarz_smoother::component_list,
     271             :                  observers::Observer<Metavariables>,
     272             :                  observers::ObserverWriter<Metavariables>>>;
     273             : 
     274             :   // Specify the transitions between phases.
     275             :   template <typename... Tags>
     276             :   static Phase determine_next_phase(
     277             :       const gsl::not_null<
     278             :           tuples::TaggedTuple<Tags...>*> /*phase_change_decision_data*/,
     279             :       const Phase& current_phase,
     280             :       const Parallel::CProxy_GlobalCache<
     281             :           Metavariables>& /*cache_proxy*/) noexcept {
     282             :     switch (current_phase) {
     283             :       case Phase::Initialization:
     284             :         return Phase::RegisterWithObserver;
     285             :       case Phase::RegisterWithObserver:
     286             :         return Phase::Solve;
     287             :       case Phase::Solve:
     288             :         return Phase::Exit;
     289             :       case Phase::Exit:
     290             :         ERROR(
     291             :             "Should never call determine_next_phase with the current phase "
     292             :             "being 'Exit'");
     293             :       default:
     294             :         ERROR(
     295             :             "Unknown type of phase. Did you static_cast<Phase> an integral "
     296             :             "value?");
     297             :     }
     298             :   }
     299             : 
     300             :   // NOLINTNEXTLINE(google-runtime-references)
     301             :   void pup(PUP::er& /*p*/) noexcept {}
     302             : };
     303             : 
     304             : static const std::vector<void (*)()> charm_init_node_funcs{
     305             :     &setup_error_handling,
     306             :     &setup_memory_allocation_failure_reporting,
     307             :     &disable_openblas_multithreading,
     308             :     &domain::creators::register_derived_with_charm,
     309             :     &Parallel::register_derived_classes_with_charm<
     310             :         metavariables::background_tag::type::element_type>,
     311             :     &Parallel::register_derived_classes_with_charm<
     312             :         metavariables::initial_guess_tag::type::element_type>,
     313             :     &Parallel::register_derived_classes_with_charm<
     314             :         Xcts::Solutions::AnalyticSolution<
     315             :             typename metavariables::analytic_solution_registrars>>,
     316             :     &Parallel::register_derived_classes_with_charm<
     317             :         metavariables::system::boundary_conditions_base>,
     318             :     &Parallel::register_derived_classes_with_charm<
     319             :         metavariables::schwarz_smoother::subdomain_solver>,
     320             :     &Parallel::register_factory_classes_with_charm<metavariables>};
     321             : static const std::vector<void (*)()> charm_init_proc_funcs{
     322             :     &enable_floating_point_exceptions};
     323             : /// \endcond

Generated by: LCOV version 1.14