SolvePoisson.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"
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/RemoveOptionsAndTerminatePhase.hpp"
41 #include "ParallelAlgorithms/LinearSolver/Gmres/Gmres.hpp"
43 #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
44 #include "PointwiseFunctions/AnalyticSolutions/Poisson/AnalyticSolution.hpp"
45 #include "PointwiseFunctions/AnalyticSolutions/Poisson/Lorentzian.hpp"
46 #include "PointwiseFunctions/AnalyticSolutions/Poisson/Moustache.hpp"
47 #include "PointwiseFunctions/AnalyticSolutions/Poisson/ProductOfSinusoids.hpp"
48 #include "PointwiseFunctions/AnalyticSolutions/Poisson/Zero.hpp"
49 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
50 #include "Utilities/Blas.hpp"
52 #include "Utilities/Functional.hpp"
53 #include "Utilities/MemoryHelpers.hpp"
54 #include "Utilities/ProtocolHelpers.hpp"
55 #include "Utilities/TMPL.hpp"
56 
57 /// \cond
58 namespace PUP {
59 class er;
60 } // namespace PUP
61 /// \endcond
62 
63 namespace SolvePoisson::OptionTags {
65  static std::string name() noexcept { return "LinearSolver"; }
66  static constexpr Options::String help =
67  "The iterative Krylov-subspace linear solver";
68 };
69 struct GmresGroup {
70  static std::string name() noexcept { return "GMRES"; }
71  static constexpr Options::String help = "Options for the GMRES linear solver";
72  using group = LinearSolverGroup;
73 };
74 } // namespace SolvePoisson::OptionTags
75 
76 /// \cond
77 template <size_t Dim>
78 struct Metavariables {
79  static constexpr size_t volume_dim = Dim;
80  using system =
82 
83  // List the possible backgrounds, i.e. the variable-independent part of the
84  // equations that define the problem to solve (along with the boundary
85  // conditions). We'll probably always have an analytic solution for Poisson
86  // problems, so we don't bother supporting non-solution backgrounds for now.
87  using analytic_solution_registrars = tmpl::flatten<tmpl::list<
88  Poisson::Solutions::Registrars::ProductOfSinusoids<Dim>,
89  tmpl::conditional_t<Dim == 1 or Dim == 2,
90  Poisson::Solutions::Registrars::Moustache<Dim>,
91  tmpl::list<>>,
92  tmpl::conditional_t<Dim == 3,
93  Poisson::Solutions::Registrars::Lorentzian<Dim>,
94  tmpl::list<>>>>;
95  using analytic_solution_tag = elliptic::Tags::Background<
97 
98  // List the possible initial guesses
99  using initial_guess_registrars =
100  tmpl::append<tmpl::list<Poisson::Solutions::Registrars::Zero<Dim>>,
101  analytic_solution_registrars>;
102  using initial_guess_tag = elliptic::Tags::InitialGuess<
104 
105  static constexpr Options::String help{
106  "Find the solution to a Poisson problem."};
107 
108  // These are the fields we solve for
110  // These are the fixed sources, i.e. the RHS of the equations
111  using fixed_sources_tag = db::add_tag_prefix<::Tags::FixedSource, fields_tag>;
112  // This is the linear operator applied to the fields. We'll only use it to
113  // apply the operator to the initial guess, so an optimization would be to
114  // re-use the `operator_applied_to_vars_tag` below. This optimization needs a
115  // few minor changes to the parallel linear solver algorithm.
116  using operator_applied_to_fields_tag =
118 
119  // The linear solver algorithm. We must use GMRES since the operator is
120  // not guaranteed to be symmetric. It can be made symmetric by multiplying by
121  // the DG mass matrix.
122  using linear_solver =
123  LinearSolver::gmres::Gmres<Metavariables, fields_tag,
125  false>;
126  using linear_solver_iteration_id =
128  // For the GMRES linear solver we need to apply the DG operator to its
129  // internal "operand" in every iteration of the algorithm.
130  using vars_tag = typename linear_solver::operand_tag;
131  using operator_applied_to_vars_tag =
133  // We'll buffer the corresponding fluxes in this tag, but won't actually need
134  // to access them outside applying the operator
135  using fluxes_vars_tag =
137  typename system::primal_fluxes>>;
138 
139  using analytic_solution_fields = typename system::primal_fields;
140  using observe_fields = analytic_solution_fields;
141 
142  // Collect all items to store in the cache.
143  using const_global_cache_tags =
144  tmpl::list<analytic_solution_tag, initial_guess_tag,
146 
147  struct factory_creation
148  : tt::ConformsTo<Options::protocols::FactoryCreation> {
149  using factory_classes = tmpl::map<
150  tmpl::pair<Event, tmpl::flatten<tmpl::list<
152  dg::Events::field_observations<
153  volume_dim, linear_solver_iteration_id,
154  observe_fields, analytic_solution_fields>>>>,
155  tmpl::pair<Trigger,
156  tmpl::push_back<Triggers::logical_triggers,
158  linear_solver_iteration_id>>>>;
159  };
160 
161  // Collect all reduction tags for observers
162  using observed_reduction_data_tags =
163  observers::collect_reduction_data_tags<tmpl::flatten<tmpl::list<
164  tmpl::at<typename factory_creation::factory_classes, Event>,
165  linear_solver>>>;
166 
167  // Specify all global synchronization points.
168  enum class Phase { Initialization, RegisterWithObserver, Solve, Exit };
169 
170  using initialization_actions = tmpl::list<
173  typename linear_solver::initialize_element,
177  analytic_solution_tag, tmpl::append<typename system::primal_fields,
178  typename system::primal_fluxes>>,
181  system, fixed_sources_tag>,
182  // Apply the DG operator to the initial guess
184  system, true, linear_solver_iteration_id, fields_tag, fluxes_vars_tag,
185  operator_applied_to_fields_tag, vars_tag, fluxes_vars_tag>,
187 
188  using build_linear_operator_actions = elliptic::dg::Actions::apply_operator<
189  system, true, linear_solver_iteration_id, vars_tag, fluxes_vars_tag,
190  operator_applied_to_vars_tag>;
191 
192  using register_actions =
195 
196  using solve_actions = tmpl::list<
197  typename linear_solver::template solve<tmpl::list<
198  Actions::RunEventsAndTriggers, build_linear_operator_actions>>,
200 
201  using dg_element_array = elliptic::DgElementArray<
202  Metavariables,
203  tmpl::list<Parallel::PhaseActions<Phase, Phase::Initialization,
204  initialization_actions>,
205  Parallel::PhaseActions<Phase, Phase::RegisterWithObserver,
206  register_actions>,
208 
209  // Specify all parallel components that will execute actions at some point.
210  using component_list = tmpl::flatten<
211  tmpl::list<dg_element_array, typename linear_solver::component_list,
214 
215  // Specify the transitions between phases.
216  template <typename... Tags>
217  static Phase determine_next_phase(
218  const gsl::not_null<
219  tuples::TaggedTuple<Tags...>*> /*phase_change_decision_data*/,
220  const Phase& current_phase,
221  const Parallel::CProxy_GlobalCache<
222  Metavariables>& /*cache_proxy*/) noexcept {
223  switch (current_phase) {
224  case Phase::Initialization:
225  return Phase::RegisterWithObserver;
226  case Phase::RegisterWithObserver:
227  return Phase::Solve;
228  case Phase::Solve:
229  return Phase::Exit;
230  case Phase::Exit:
231  ERROR(
232  "Should never call determine_next_phase with the current phase "
233  "being 'Exit'");
234  default:
235  ERROR(
236  "Unknown type of phase. Did you static_cast<Phase> an integral "
237  "value?");
238  }
239  }
240 
241  // NOLINTNEXTLINE(google-runtime-references)
242  void pup(PUP::er& /*p*/) noexcept {}
243 };
244 
245 static const std::vector<void (*)()> charm_init_node_funcs{
246  &setup_error_handling, &setup_memory_allocation_failure_reporting,
248  &domain::creators::register_derived_with_charm,
250  metavariables::analytic_solution_tag::type::element_type>,
252  metavariables::initial_guess_tag::type::element_type>,
254  metavariables::system::boundary_conditions_base>,
255  &Parallel::register_factory_classes_with_charm<metavariables>};
256 static const std::vector<void (*)()> charm_init_proc_funcs{
258 /// \endcond
FloatingPointExceptions.hpp
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
std::string
SolvePoisson::OptionTags::GmresGroup
Definition: SolvePoisson.hpp:69
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
elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource
For linear systems, impose inhomogeneous boundary conditions as contributions to the fixed sources (i...
Definition: ApplyOperator.hpp:523
GlobalCache.hpp
Options.hpp
LinearSolver::Tags::Operand
The operand that the local linear operator is applied to.
Definition: Tags.hpp:40
Tags.hpp
std::vector
std::system
T system(T... args)
elliptic::Triggers::EveryNIterations
Definition: EveryNIterations.hpp:19
elliptic::Actions::InitializeAnalyticSolution
Place the analytic solution of the system fields in the DataBox.
Definition: InitializeAnalyticSolution.hpp:58
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
SolvePoisson::OptionTags::LinearSolverGroup
Definition: SolvePoisson.hpp:64
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()
Poisson::FirstOrderSystem
The Poisson equation formulated as a set of coupled first-order PDEs.
Definition: FirstOrderSystem.hpp:77
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
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
FirstOrderSystem.hpp
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
cstddef
Actions::RunEventsAndTriggers
Run the events and triggers.
Definition: RunEventsAndTriggers.hpp:27
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
Poisson::Solutions::AnalyticSolution
Base class for analytic solutions of the Poisson equation.
Definition: AnalyticSolution.hpp:32
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
Tags::EventsAndTriggers
Definition: Tags.hpp:51
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
db::wrap_tags_in
tmpl::transform< TagList, tmpl::bind< Wrapper, tmpl::_1, tmpl::pin< Args >... > > wrap_tags_in
Create a new tmpl::list of tags by wrapping each tag in TagList in Wrapper<_, Args....
Definition: PrefixHelpers.hpp:30
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
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
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13