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

Generated by: LCOV version 1.14