SolveLinearEllipticProblem.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/InitializeSystem.hpp"
13 #include "Elliptic/DiscontinuousGalerkin/DgElementArray.hpp"
14 #include "Elliptic/DiscontinuousGalerkin/ImposeBoundaryConditions.hpp"
15 #include "Elliptic/DiscontinuousGalerkin/ImposeInhomogeneousBoundaryConditionsOnSource.hpp"
16 #include "Elliptic/DiscontinuousGalerkin/InitializeFirstOrderOperator.hpp"
17 #include "Elliptic/DiscontinuousGalerkin/NumericalFluxes/FirstOrderInternalPenalty.hpp"
18 #include "Elliptic/FirstOrderOperator.hpp"
20 #include "Elliptic/Tags.hpp"
21 #include "Elliptic/Triggers/EveryNIterations.hpp"
23 #include "IO/Observer/Actions.hpp"
24 #include "IO/Observer/Helpers.hpp"
25 #include "IO/Observer/ObserverComponent.hpp"
26 #include "IO/Observer/RegisterObservers.hpp"
27 #include "NumericalAlgorithms/DiscontinuousGalerkin/BoundarySchemes/FirstOrder/FirstOrderScheme.hpp"
28 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
29 #include "Options/Options.hpp"
30 #include "Parallel/Actions/TerminatePhase.hpp"
32 #include "Parallel/InitializationFunctions.hpp"
33 #include "Parallel/PhaseDependentActionList.hpp"
34 #include "Parallel/Reduction.hpp"
36 #include "ParallelAlgorithms/Actions/MutateApply.hpp"
37 #include "ParallelAlgorithms/DiscontinuousGalerkin/CollectDataForFluxes.hpp"
38 #include "ParallelAlgorithms/DiscontinuousGalerkin/FluxCommunication.hpp"
39 #include "ParallelAlgorithms/DiscontinuousGalerkin/InitializeDomain.hpp"
40 #include "ParallelAlgorithms/DiscontinuousGalerkin/InitializeInterfaces.hpp"
41 #include "ParallelAlgorithms/DiscontinuousGalerkin/InitializeMortars.hpp"
42 #include "ParallelAlgorithms/Events/ObserveErrorNorms.hpp"
43 #include "ParallelAlgorithms/Events/ObserveFields.hpp"
44 #include "ParallelAlgorithms/EventsAndTriggers/Actions/RunEventsAndTriggers.hpp"
45 #include "ParallelAlgorithms/Initialization/Actions/RemoveOptionsAndTerminatePhase.hpp"
46 #include "ParallelAlgorithms/LinearSolver/Actions/TerminateIfConverged.hpp"
47 #include "ParallelAlgorithms/LinearSolver/Gmres/Gmres.hpp"
49 #include "PointwiseFunctions/AnalyticSolutions/Poisson/Lorentzian.hpp"
50 #include "PointwiseFunctions/AnalyticSolutions/Poisson/Moustache.hpp"
51 #include "PointwiseFunctions/AnalyticSolutions/Poisson/ProductOfSinusoids.hpp"
52 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
53 #include "Utilities/Functional.hpp"
54 #include "Utilities/TMPL.hpp"
55 
56 namespace SolveLinearEllipticProblem {
57 namespace OptionTags {
59  static std::string name() noexcept { return "LinearSolver"; }
60  static constexpr OptionString help =
61  "Options for the iterative linear solver";
62 };
63 } // namespace OptionTags
64 } // namespace SolveLinearEllipticProblem
65 
66 /// \cond
67 template <typename System, typename InitialGuess, typename BoundaryConditions>
68 struct Metavariables {
69  using system = System;
70  static constexpr size_t volume_dim = system::volume_dim;
71  using initial_guess = InitialGuess;
72  using boundary_conditions = BoundaryConditions;
73 
74  static constexpr OptionString help{
75  "Find the solution to a linear elliptic problem.\n"
76  "Linear solver: GMRES\n"
77  "Numerical flux: FirstOrderInternalPenaltyFlux"};
78 
79  using fluxes_computer_tag =
81 
82  // Only Dirichlet boundary conditions are currently supported, and they are
83  // are all imposed by analytic solutions right now.
84  // This will be generalized ASAP. We will also support numeric initial guesses
85  // and analytic initial guesses that aren't solutions ("analytic data").
86  using analytic_solution_tag = Tags::AnalyticSolution<boundary_conditions>;
87 
88  // The linear solver algorithm. We must use GMRES since the operator is
89  // not positive-definite for the first-order system.
90  using linear_solver = LinearSolver::gmres::Gmres<
91  Metavariables, typename system::fields_tag,
93  using linear_solver_iteration_id =
95  // For the GMRES linear solver we need to apply the DG operator to its
96  // internal "operand" in every iteration of the algorithm.
97  using linear_operand_tag = db::add_tag_prefix<LinearSolver::Tags::Operand,
98  typename system::fields_tag>;
99  using primal_variables = db::wrap_tags_in<LinearSolver::Tags::Operand,
100  typename system::primal_fields>;
101  using auxiliary_variables =
103  typename system::auxiliary_fields>;
104 
105  // Parse numerical flux parameters from the input file to store in the cache.
106  using normal_dot_numerical_flux = Tags::NumericalFlux<
108  volume_dim, fluxes_computer_tag, primal_variables,
109  auxiliary_variables>>;
110  // Specify the DG boundary scheme. We use the strong first-order scheme here
111  // that only requires us to compute normals dotted into the first-order
112  // fluxes.
113  using boundary_scheme =
114  dg::FirstOrderScheme::FirstOrderScheme<volume_dim, linear_operand_tag,
115  normal_dot_numerical_flux,
116  linear_solver_iteration_id>;
117 
118  // Collect events and triggers
119  // (public for use by the Charm++ registration code)
120  using observe_fields =
121  db::get_variables_tags_list<typename system::fields_tag>;
122  using analytic_solution_fields = observe_fields;
123  using events =
125  volume_dim, linear_solver_iteration_id, observe_fields,
126  analytic_solution_fields>,
128  linear_solver_iteration_id, analytic_solution_fields>>;
129  using triggers = tmpl::list<elliptic::Triggers::Registrars::EveryNIterations<
130  linear_solver_iteration_id>>;
131 
132  // Collect all items to store in the cache.
133  using const_global_cache_tags =
134  tmpl::list<analytic_solution_tag, fluxes_computer_tag,
135  normal_dot_numerical_flux,
137 
138  // Collect all reduction tags for observers
139  struct element_observation_type {};
140  using observed_reduction_data_tags =
141  observers::collect_reduction_data_tags<tmpl::flatten<tmpl::list<
142  typename Event<events>::creatable_classes, linear_solver>>>;
143 
144  // Specify all global synchronization points.
145  enum class Phase { Initialization, RegisterWithObserver, Solve, Exit };
146 
147  using initialization_actions = tmpl::list<
154  typename linear_solver::initialize_element,
156  elliptic::Actions::InitializeAnalyticSolution<analytic_solution_tag,
157  analytic_solution_fields>,
159  Metavariables>,
162  volume_dim, typename system::fluxes, typename system::sources,
163  linear_operand_tag, primal_variables, auxiliary_variables>,
165 
166  using build_linear_operator_actions = tmpl::list<
172  linear_operand_tag>>,
174  linear_operand_tag, primal_variables>,
176  boundary_scheme,
180 
181  using dg_element_array = elliptic::DgElementArray<
182  Metavariables,
183  tmpl::list<Parallel::PhaseActions<Phase, Phase::Initialization,
184  initialization_actions>,
186  Phase, Phase::RegisterWithObserver,
189  linear_solver_iteration_id,
190  element_observation_type>>,
191  // We prepare the linear solve here to avoid
192  // adding an extra phase. We can't do that
193  // before registration because the
194  // `prepare_solve` action may contribute to
195  // observers.
196  typename linear_solver::prepare_solve,
199  Phase, Phase::Solve,
202  typename linear_solver::options_group>,
203  typename linear_solver::prepare_step,
204  build_linear_operator_actions,
205  typename linear_solver::perform_step>>>>;
206 
207  // Specify all parallel components that will execute actions at some point.
208  using component_list = tmpl::flatten<
209  tmpl::list<dg_element_array, typename linear_solver::component_list,
212 
213  // Specify the transitions between phases.
214  static Phase determine_next_phase(
215  const Phase& current_phase,
216  const Parallel::CProxy_ConstGlobalCache<
217  Metavariables>& /*cache_proxy*/) noexcept {
218  switch (current_phase) {
219  case Phase::Initialization:
220  return Phase::RegisterWithObserver;
221  case Phase::RegisterWithObserver:
222  return Phase::Solve;
223  case Phase::Solve:
224  return Phase::Exit;
225  case Phase::Exit:
226  ERROR(
227  "Should never call determine_next_phase with the current phase "
228  "being 'Exit'");
229  default:
230  ERROR(
231  "Unknown type of phase. Did you static_cast<Phase> an integral "
232  "value?");
233  }
234  }
235 };
236 
237 static const std::vector<void (*)()> charm_init_node_funcs{
238  &setup_error_handling, &domain::creators::register_derived_with_charm,
243 static const std::vector<void (*)()> charm_init_proc_funcs{
245 /// \endcond
FloatingPointExceptions.hpp
elliptic::Actions::InitializeSystem
Initializes the DataBox tags related to the elliptic system.
Definition: InitializeSystem.hpp:47
RegisterDerivedClassesWithCharm.hpp
LinearSolver::Actions::TerminateIfConverged
Terminate the algorithm if the linear solver has converged.
Definition: TerminateIfConverged.hpp:33
LinearSolver::gmres::Gmres
A GMRES solver for nonsymmetric linear systems of equations .
Definition: Gmres.hpp:80
std::string
elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource
Adds boundary contributions to the sources.
Definition: ImposeInhomogeneousBoundaryConditionsOnSource.hpp:72
LinearSolver::Tags::OperatorAppliedTo
The linear operator applied to the data in Tag
Definition: Tags.hpp:70
Options.hpp
LinearSolver::Tags::Operand
The operand that the local linear operator is applied to.
Definition: Tags.hpp:57
Tags.hpp
std::vector
std::system
T system(T... args)
elliptic::Actions::InitializeAnalyticSolution
Place the analytic solution of the system fields in the DataBox.
Definition: InitializeAnalyticSolution.hpp:46
LinearSolver::Tags::IterationId
Holds an IterationId that identifies a step in the linear solver algorithm.
Definition: Tags.hpp:84
dg::Actions::SendDataForFluxes
Send local boundary data needed for fluxes to neighbors.
Definition: FluxCommunication.hpp:79
Initialization::Actions::RemoveOptionsAndTerminatePhase
Definition: RemoveOptionsAndTerminatePhase.hpp:27
dg::Initialization::face_compute_tags
tmpl::list< Tags... > face_compute_tags
Definition: InitializeInterfaces.hpp:42
Parallel::Actions::TerminatePhase
Terminate the algorithm to proceed to the next phase.
Definition: TerminatePhase.hpp:26
Tags::NumericalFlux
The global cache tag for the numerical flux.
Definition: Tags.hpp:88
dg::Initialization::slice_tags_to_face
tmpl::list< Tags... > slice_tags_to_face
Definition: InitializeInterfaces.hpp:32
domain::Tags::BoundaryDirectionsInterior
Definition: Tags.hpp:250
db::add_tag_prefix
Prefix< DataBox_detail::dispatch_add_tag_prefix_impl< Prefix, Tag, Args... >, Args... > add_tag_prefix
Definition: PrefixHelpers.hpp:145
Tags.hpp
Registration::Registrar
A template for defining a registrar.
Definition: Registration.hpp:42
elliptic::FirstOrderOperator
Mutating DataBox invokable to compute the bulk contribution to the operator represented by the Operat...
Definition: FirstOrderOperator.hpp:201
enable_floating_point_exceptions
void enable_floating_point_exceptions()
Definition: FloatingPointExceptions.cpp:27
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
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:48
dg::Actions::InitializeInterfaces
Initialize items related to the interfaces between Elements and on external boundaries.
Definition: InitializeInterfaces.hpp:141
Event
Definition: Event.hpp:30
elliptic::DgElementArray
The parallel component responsible for managing the DG elements that compose the computational domain...
Definition: DgElementArray.hpp:37
dg::Actions::InitializeDomain
Initialize items related to the basic structure of the element.
Definition: InitializeDomain.hpp:64
FirstOrderSystem.hpp
cstddef
Actions::RunEventsAndTriggers
Run the events and triggers.
Definition: RunEventsAndTriggers.hpp:27
Actions::MutateApply
Apply the function Mutator::apply to the DataBox.
Definition: MutateApply.hpp:40
elliptic::Tags::FluxesComputer
Holds an object that computes the principal part of the elliptic PDEs.
Definition: Tags.hpp:27
Tags::AnalyticSolution
Definition: Tags.hpp:52
observers::Observer
The group parallel component that is responsible for reducing data to be observed.
Definition: ObserverComponent.hpp:27
Trigger
Definition: Trigger.hpp:34
dg::Actions::InitializeMortars
Initialize mortars between elements for exchanging fluxes.
Definition: InitializeMortars.hpp:75
dg::Actions::CollectDataForFluxes
Collect data that is needed to compute numerical fluxes and store it on mortars, projecting it if nec...
Definition: CollectDataForFluxes.hpp:68
dg::Events::Registrars::ObserveFields
Definition: ObserveFields.hpp:62
elliptic::dg::NumericalFluxes::FirstOrderInternalPenalty
The internal penalty flux for first-order elliptic equations.
Definition: FirstOrderInternalPenalty.hpp:136
Tags::EventsAndTriggers
Definition: Tags.hpp:57
dg::Initialization::exterior_compute_tags
tmpl::list< Tags... > exterior_compute_tags
Definition: InitializeInterfaces.hpp:47
observers::RegisterObservers
Passed to RegisterWithObservers action to register observer event.
Definition: RegisterObservers.hpp:15
dg::FirstOrderScheme::FirstOrderScheme
Boundary contributions for a first-order DG scheme.
Definition: FirstOrderScheme.hpp:67
domain::Tags::InternalDirections
Definition: Tags.hpp:231
db::wrap_tags_in
tmpl::transform< TagList, tmpl::bind< Wrapper, tmpl::_1, tmpl::pin< Args >... > > wrap_tags_in
Create a new list of Tags by wrapping each tag in TagList using the Wrapper.
Definition: PrefixHelpers.hpp:30
Parallel::register_derived_classes_with_charm
void register_derived_classes_with_charm() noexcept
Register derived classes of the Base class.
Definition: RegisterDerivedClassesWithCharm.hpp:31
elliptic::dg::Actions::ImposeHomogeneousDirichletBoundaryConditions
Set field data on external boundaries so that they represent homogeneous (zero) Dirichlet boundary co...
Definition: ImposeBoundaryConditions.hpp:91
SolveLinearEllipticProblem::OptionTags::LinearSolverGroup
Definition: SolveLinearEllipticProblem.hpp:58
elliptic::dg::Actions::InitializeFirstOrderOperator
Initialize DataBox tags for building the first-order elliptic DG operator.
Definition: InitializeFirstOrderOperator.hpp:45
dg::Actions::ReceiveDataForFluxes
Receive boundary data needed for fluxes from neighbors.
Definition: FluxCommunication.hpp:177
OptionString
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:30
TMPL.hpp
ConstGlobalCache.hpp
dg::Initialization::slice_tags_to_exterior
tmpl::list< Tags... > slice_tags_to_exterior
Definition: InitializeInterfaces.hpp:37
observers::Actions::RegisterWithObservers
Registers itself with the local observer parallel component so the observer knows to expect data from...
Definition: Actions.hpp:232