SolveXcts.hpp
1 // 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/RegisterDerivedWithCharm.hpp"
10 #include "Domain/Tags.hpp"
11 #include "Elliptic/Actions/InitializeAnalyticSolution.hpp"
12 #include "Elliptic/Actions/InitializeFields.hpp"
13 #include "Elliptic/Actions/InitializeFixedSources.hpp"
14 #include "Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp"
15 #include "Elliptic/DiscontinuousGalerkin/Actions/InitializeDomain.hpp"
16 #include "Elliptic/DiscontinuousGalerkin/DgElementArray.hpp"
17 #include "Elliptic/Systems/Xcts/FirstOrderSystem.hpp"
18 #include "Elliptic/Tags.hpp"
19 #include "Elliptic/Triggers/EveryNIterations.hpp"
20 #include "IO/Observer/Actions/RegisterEvents.hpp"
21 #include "IO/Observer/Helpers.hpp"
22 #include "IO/Observer/ObserverComponent.hpp"
23 #include "NumericalAlgorithms/Convergence/Tags.hpp"
24 #include "Options/Options.hpp"
25 #include "Options/Protocols/FactoryCreation.hpp"
26 #include "Parallel/Actions/SetupDataBox.hpp"
27 #include "Parallel/Actions/TerminatePhase.hpp"
28 #include "Parallel/GlobalCache.hpp"
29 #include "Parallel/InitializationFunctions.hpp"
30 #include "Parallel/PhaseDependentActionList.hpp"
31 #include "Parallel/Reduction.hpp"
33 #include "ParallelAlgorithms/Events/Factory.hpp"
34 #include "ParallelAlgorithms/EventsAndTriggers/Actions/RunEventsAndTriggers.hpp"
35 #include "ParallelAlgorithms/EventsAndTriggers/Completion.hpp"
36 #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
37 #include "ParallelAlgorithms/EventsAndTriggers/LogicalTriggers.hpp"
39 #include "ParallelAlgorithms/EventsAndTriggers/Trigger.hpp"
40 #include "ParallelAlgorithms/Initialization/Actions/AddComputeTags.hpp"
41 #include "ParallelAlgorithms/Initialization/Actions/RemoveOptionsAndTerminatePhase.hpp"
42 #include "ParallelAlgorithms/LinearSolver/Gmres/Gmres.hpp"
44 #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/NewtonRaphson.hpp"
45 #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
46 #include "PointwiseFunctions/AnalyticData/Xcts/AnalyticData.hpp"
47 #include "PointwiseFunctions/AnalyticData/Xcts/Binary.hpp"
48 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
49 #include "PointwiseFunctions/AnalyticSolutions/Xcts/AnalyticSolution.hpp"
50 #include "PointwiseFunctions/AnalyticSolutions/Xcts/Flatness.hpp"
51 #include "PointwiseFunctions/AnalyticSolutions/Xcts/Kerr.hpp"
52 #include "PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.hpp"
53 #include "Utilities/Blas.hpp"
55 #include "Utilities/Functional.hpp"
56 #include "Utilities/MemoryHelpers.hpp"
57 #include "Utilities/ProtocolHelpers.hpp"
58 #include "Utilities/TMPL.hpp"
59 
60 namespace SolveXcts::OptionTags {
62  static std::string name() noexcept { return "NonlinearSolver"; }
63  static constexpr Options::String help = "The iterative nonlinear solver";
64 };
66  static std::string name() noexcept { return "NewtonRaphson"; }
67  static constexpr Options::String help =
68  "Options for the Newton-Raphson nonlinear solver";
70 };
72  static std::string name() noexcept { return "LinearSolver"; }
73  static constexpr Options::String help =
74  "The iterative Krylov-subspace linear solver";
75 };
76 struct GmresGroup {
77  static std::string name() noexcept { return "Gmres"; }
78  static constexpr Options::String help = "Options for the GMRES linear solver";
79  using group = LinearSolverGroup;
80 };
81 } // namespace SolveXcts::OptionTags
82 
83 /// \cond
84 struct Metavariables {
85  static constexpr size_t volume_dim = 3;
86  static constexpr int conformal_matter_scale = 0;
87  using system =
89  Xcts::Geometry::Curved, conformal_matter_scale>;
90 
91  // List the possible backgrounds, i.e. the variable-independent part of the
92  // equations that define the problem to solve (along with the boundary
93  // conditions)
94  using analytic_solution_registrars =
95  tmpl::list<Xcts::Solutions::Registrars::Schwarzschild,
96  Xcts::Solutions::Registrars::Kerr>;
97  using analytic_data_registrars = tmpl::list<
98  Xcts::AnalyticData::Registrars::Binary<analytic_solution_registrars>>;
99  using background_tag = elliptic::Tags::Background<
100  ::AnalyticData<volume_dim, tmpl::append<analytic_solution_registrars,
101  analytic_data_registrars>>>;
102 
103  // List the possible initial guesses
104  using initial_guess_registrars =
105  tmpl::append<tmpl::list<Xcts::Solutions::Registrars::Flatness>,
106  analytic_solution_registrars, analytic_data_registrars>;
107  using initial_guess_tag = elliptic::Tags::InitialGuess<
109 
110  static constexpr Options::String help{
111  "Find the solution to an XCTS problem."};
112 
113  // These are the fields we solve for
115  // These are the fluxes corresponding to the fields, i.e. essentially their
116  // first derivatives. These are background fields for the linearized sources.
118  // These are the fixed sources, i.e. the RHS of the equations
119  using fixed_sources_tag = db::add_tag_prefix<::Tags::FixedSource, fields_tag>;
120  using operator_applied_to_fields_tag =
122 
123  using nonlinear_solver = NonlinearSolver::newton_raphson::NewtonRaphson<
124  Metavariables, fields_tag, SolveXcts::OptionTags::NewtonRaphsonGroup,
125  fixed_sources_tag>;
126  using nonlinear_solver_iteration_id =
128 
129  // The linear solver algorithm. We must use GMRES since the operator is
130  // not guaranteed to be symmetric.
131  using linear_solver = LinearSolver::gmres::Gmres<
132  Metavariables, typename nonlinear_solver::linear_solver_fields_tag,
134  typename nonlinear_solver::linear_solver_source_tag>;
135  using linear_solver_iteration_id =
137  // For the GMRES linear solver we need to apply the DG operator to its
138  // internal "operand" in every iteration of the algorithm.
139  using correction_vars_tag = typename linear_solver::operand_tag;
140  using operator_applied_to_correction_vars_tag =
142  correction_vars_tag>;
143  // The correction fluxes can be stored in an arbitrary tag. We don't need to
144  // access them anywhere, they're just a memory buffer for the linearized
145  // operator.
146  using correction_fluxes_tag =
148 
149  using analytic_solution_fields = tmpl::append<typename system::primal_fields,
150  typename system::primal_fluxes>;
151  using observe_fields = tmpl::append<analytic_solution_fields,
152  typename system::background_fields>;
153 
154  // Collect all items to store in the cache.
155  using const_global_cache_tags =
156  tmpl::list<background_tag, initial_guess_tag, Tags::EventsAndTriggers>;
157 
158  struct factory_creation
159  : tt::ConformsTo<Options::protocols::FactoryCreation> {
160  using factory_classes = tmpl::map<
161  tmpl::pair<Event, tmpl::flatten<tmpl::list<
163  dg::Events::field_observations<
164  volume_dim, nonlinear_solver_iteration_id,
165  observe_fields, analytic_solution_fields>>>>,
166  tmpl::pair<Trigger,
167  tmpl::push_back<Triggers::logical_triggers,
169  nonlinear_solver_iteration_id>>>>;
170  };
171 
172  // Collect all reduction tags for observers
173  using observed_reduction_data_tags =
174  observers::collect_reduction_data_tags<tmpl::flatten<
175  tmpl::list<tmpl::at<factory_creation::factory_classes, Event>,
176  nonlinear_solver, linear_solver>>>;
177 
178  // Specify all global synchronization points.
179  enum class Phase { Initialization, RegisterWithObserver, Solve, Exit };
180 
181  using initialization_actions = tmpl::list<
184  typename nonlinear_solver::initialize_element,
185  typename linear_solver::initialize_element,
189  background_tag, analytic_solution_fields,
191  analytic_solution_registrars, analytic_data_registrars>>>,
194 
195  template <bool Linearized>
196  using build_operator_actions = elliptic::dg::Actions::apply_operator<
197  system, Linearized,
198  tmpl::conditional_t<Linearized, linear_solver_iteration_id,
199  nonlinear_solver_iteration_id>,
200  tmpl::conditional_t<Linearized, correction_vars_tag, fields_tag>,
201  tmpl::conditional_t<Linearized, correction_fluxes_tag, fluxes_tag>,
202  tmpl::conditional_t<Linearized, operator_applied_to_correction_vars_tag,
203  operator_applied_to_fields_tag>>;
204 
205  using register_actions =
207  typename nonlinear_solver::register_element,
209 
210  using solve_actions = tmpl::list<
211  typename nonlinear_solver::template solve<
212  build_operator_actions<false>,
213  typename linear_solver::template solve<build_operator_actions<true>>,
216 
217  using dg_element_array = elliptic::DgElementArray<
218  Metavariables,
219  tmpl::list<Parallel::PhaseActions<Phase, Phase::Initialization,
220  initialization_actions>,
221  Parallel::PhaseActions<Phase, Phase::RegisterWithObserver,
222  register_actions>,
224 
225  // Specify all parallel components that will execute actions at some point.
226  using component_list = tmpl::flatten<
227  tmpl::list<dg_element_array, typename nonlinear_solver::component_list,
228  typename linear_solver::component_list,
231 
232  // Specify the transitions between phases.
233  template <typename... Tags>
234  static Phase determine_next_phase(
235  const gsl::not_null<
236  tuples::TaggedTuple<Tags...>*> /*phase_change_decision_data*/,
237  const Phase& current_phase,
238  const Parallel::CProxy_GlobalCache<
239  Metavariables>& /*cache_proxy*/) noexcept {
240  switch (current_phase) {
241  case Phase::Initialization:
242  return Phase::RegisterWithObserver;
243  case Phase::RegisterWithObserver:
244  return Phase::Solve;
245  case Phase::Solve:
246  return Phase::Exit;
247  case Phase::Exit:
248  ERROR(
249  "Should never call determine_next_phase with the current phase "
250  "being 'Exit'");
251  default:
252  ERROR(
253  "Unknown type of phase. Did you static_cast<Phase> an integral "
254  "value?");
255  }
256  }
257 
258  // NOLINTNEXTLINE(google-runtime-references)
259  void pup(PUP::er& /*p*/) noexcept {}
260 };
261 
262 static const std::vector<void (*)()> charm_init_node_funcs{
263  &setup_error_handling, &setup_memory_allocation_failure_reporting,
265  &domain::creators::register_derived_with_charm,
267  metavariables::background_tag::type::element_type>,
269  metavariables::initial_guess_tag::type::element_type>,
272  typename metavariables::analytic_solution_registrars>>,
274  metavariables::system::boundary_conditions_base>,
275  &Parallel::register_factory_classes_with_charm<metavariables>};
276 static const std::vector<void (*)()> charm_init_proc_funcs{
278 /// \endcond
FloatingPointExceptions.hpp
Xcts::FirstOrderSystem
The Extended Conformal Thin Sandwich (XCTS) decomposition of the Einstein constraint equations,...
Definition: FirstOrderSystem.hpp:142
RegisterDerivedClassesWithCharm.hpp
Actions::SetupDataBox
Add into the DataBox default constructed items for the collection of tags requested by any of the act...
Definition: SetupDataBox.hpp:102
LinearSolver::gmres::Gmres
A GMRES solver for nonsymmetric linear systems of equations .
Definition: Gmres.hpp:87
SolveXcts::OptionTags::LinearSolverGroup
Definition: SolveXcts.hpp:71
std::string
observers::Actions::RegisterEventsWithObservers
Registers this element of a parallel component with the local Observer parallel component for each tr...
Definition: RegisterEvents.hpp:108
Tags.hpp
LinearSolver::Tags::OperatorAppliedTo
The linear operator applied to the data in Tag
Definition: Tags.hpp:53
GlobalCache.hpp
Options.hpp
Tags.hpp
std::vector
std::system
T system(T... args)
elliptic::Triggers::EveryNIterations
Definition: EveryNIterations.hpp:19
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
Tags::Variables
Definition: VariablesTag.hpp:21
Initialization::Actions::RemoveOptionsAndTerminatePhase
Definition: RemoveOptionsAndTerminatePhase.hpp:27
Parallel::Actions::TerminatePhase
Terminate the algorithm to proceed to the next phase.
Definition: TerminatePhase.hpp:26
elliptic::Tags::Background
The variable-independent part of the elliptic equations, e.g. the fixed-sources in a Poisson equatio...
Definition: Tags.hpp:44
disable_openblas_multithreading
void disable_openblas_multithreading() noexcept
Disable OpenBLAS multithreading since it conflicts with Charm++ parallelism.
Tags.hpp
domain::push_back
CoordinateMap< SourceFrame, TargetFrame, Maps..., NewMap > push_back(CoordinateMap< SourceFrame, TargetFrame, Maps... > old_map, NewMap new_map) noexcept
Creates a CoordinateMap by appending the new map to the end of the old maps.
enable_floating_point_exceptions
void enable_floating_point_exceptions()
Xcts::Geometry::Curved
@ Curved
The conformal geometry is either curved or employs curved coordinates, so non-vanishing Christoffel s...
Events::Completion
Definition: Completion.hpp:15
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:37
Xcts::Equations::HamiltonianLapseAndShift
@ HamiltonianLapseAndShift
The full XCTS equations, solved for , and .
Parallel::PhaseActions
List of all the actions to be executed in the specified phase.
Definition: PhaseDependentActionList.hpp:16
observers::ObserverWriter
The nodegroup parallel component that is responsible for writing data to disk.
Definition: ObserverComponent.hpp:51
Event
Definition: Event.hpp:19
elliptic::DgElementArray
The parallel component responsible for managing the DG elements that compose the computational domain...
Definition: DgElementArray.hpp:37
elliptic::dg::Actions::apply_operator
tmpl::list< detail::PrepareAndSendMortarData< System, Linearized, TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag, OperatorAppliedToFieldsTag, PrimalMortarFieldsTag, PrimalMortarFluxesTag >, detail::ReceiveMortarDataAndApplyOperator< System, Linearized, TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag, OperatorAppliedToFieldsTag, PrimalMortarFieldsTag, PrimalMortarFluxesTag > > apply_operator
Apply the DG operator to the PrimalFieldsTag and write the result to the OperatorAppliedToFieldsTag
Definition: ApplyOperator.hpp:511
elliptic::Actions::InitializeOptionalAnalyticSolution
Definition: InitializeAnalyticSolution.hpp:87
cstddef
Actions::RunEventsAndTriggers
Run the events and triggers.
Definition: RunEventsAndTriggers.hpp:27
SolveXcts::OptionTags::NonlinearSolverGroup
Definition: SolveXcts.hpp:61
Xcts::Solutions::AnalyticSolution
Base class for analytic solutions of the XCTS equations.
Definition: AnalyticSolution.hpp:34
SolveXcts::OptionTags::NewtonRaphsonGroup
Definition: SolveXcts.hpp:65
tuples::TaggedTuple< Tags... >
observers::Observer
The group parallel component that is responsible for reducing data to be observed.
Definition: ObserverComponent.hpp:29
Trigger
Definition: Trigger.hpp:17
elliptic::dg::Actions::initialize_operator
tmpl::list< detail::InitializeFacesMortarsAndBackground< System, BackgroundTag > > initialize_operator
Initialize geometric and background quantities for the elliptic DG operator.
Definition: ApplyOperator.hpp:473
AnalyticData
Provides analytic tensor data as a function of the spatial coordinates.
Definition: AnalyticData.hpp:34
elliptic::Actions::InitializeFields
Initialize the dynamic fields of the elliptic system, i.e. those we solve for.
Definition: InitializeFields.hpp:43
Convergence::Tags::IterationId
Identifies a step in an iterative algorithm.
Definition: Tags.hpp:74
elliptic::dg::Actions::InitializeDomain
Initialize items related to the basic structure of the element.
Definition: InitializeDomain.hpp:53
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
elliptic::Tags::InitialGuess
The initial guess for the elliptic solve.
Definition: Tags.hpp:56
Parallel::register_derived_classes_with_charm
void register_derived_classes_with_charm() noexcept
Register derived classes of the Base class.
Definition: RegisterDerivedClassesWithCharm.hpp:35
Blas.hpp
NonlinearSolver::newton_raphson::NewtonRaphson
A Newton-Raphson correction scheme for nonlinear systems of equations .
Definition: NewtonRaphson.hpp:73
tt::ConformsTo
Indicate a class conforms to the Protocol.
Definition: ProtocolHelpers.hpp:22
elliptic::Actions::InitializeFixedSources
Initialize the "fixed sources" of the elliptic equations, i.e. their variable-independent source term...
Definition: InitializeFixedSources.hpp:54
TMPL.hpp
SolveXcts::OptionTags::GmresGroup
Definition: SolveXcts.hpp:76
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13