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/InitializeFluxes.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/Actions/ApplyFluxes.hpp"
28 #include "NumericalAlgorithms/DiscontinuousGalerkin/Actions/ComputeNonconservativeBoundaryFluxes.hpp"
30 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
31 #include "Options/Options.hpp"
32 #include "Parallel/Actions/TerminatePhase.hpp"
34 #include "Parallel/InitializationFunctions.hpp"
35 #include "Parallel/PhaseDependentActionList.hpp"
36 #include "Parallel/Reduction.hpp"
38 #include "ParallelAlgorithms/Actions/MutateApply.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 /// \cond
57 template <typename System, typename InitialGuess, typename BoundaryConditions>
58 struct Metavariables {
59  using system = System;
60  static constexpr size_t volume_dim = system::volume_dim;
61  using initial_guess = InitialGuess;
62  using boundary_conditions = BoundaryConditions;
63 
64  static constexpr OptionString help{
65  "Find the solution to a linear elliptic problem.\n"
66  "Linear solver: GMRES\n"
67  "Numerical flux: FirstOrderInternalPenaltyFlux"};
68 
69  using fluxes_computer_tag =
71 
72  // Only Dirichlet boundary conditions are currently supported, and they are
73  // are all imposed by analytic solutions right now.
74  // This will be generalized ASAP. We will also support numeric initial guesses
75  // and analytic initial guesses that aren't solutions ("analytic data").
76  using analytic_solution_tag = Tags::AnalyticSolution<boundary_conditions>;
77 
78  // The linear solver algorithm. We must use GMRES since the operator is
79  // not positive-definite for the first-order system.
80  using linear_solver =
82  using temporal_id = LinearSolver::Tags::IterationId;
83 
84  // This is needed for InitializeMortars and will be removed ASAP.
85  static constexpr bool local_time_stepping = false;
86 
87  // Parse numerical flux parameters from the input file to store in the cache.
88  using normal_dot_numerical_flux = Tags::NumericalFlux<
90  volume_dim, fluxes_computer_tag, typename system::primal_variables,
91  typename system::auxiliary_variables>>;
92 
93  // Collect events and triggers
94  // (public for use by the Charm++ registration code)
95  using observe_fields =
96  db::get_variables_tags_list<typename system::fields_tag>;
97  using analytic_solution_fields = observe_fields;
98  using events = tmpl::list<
100  volume_dim, LinearSolver::Tags::IterationId, observe_fields,
101  analytic_solution_fields>,
102  dg::Events::Registrars::ObserveErrorNorms<LinearSolver::Tags::IterationId,
103  analytic_solution_fields>>;
104  using triggers = tmpl::list<elliptic::Triggers::Registrars::EveryNIterations<
105  LinearSolver::Tags::IterationId>>;
106 
107  // Collect all items to store in the cache.
108  using const_global_cache_tags =
109  tmpl::list<analytic_solution_tag, fluxes_computer_tag,
111 
112  // Collect all reduction tags for observers
113  struct element_observation_type {};
114  using observed_reduction_data_tags =
115  observers::collect_reduction_data_tags<tmpl::flatten<tmpl::list<
116  typename Event<events>::creatable_classes, linear_solver>>>;
117 
118  // Specify all global synchronization points.
119  enum class Phase { Initialization, RegisterWithObserver, Solve, Exit };
120 
121  using initialization_actions = tmpl::list<
124  elliptic::Actions::InitializeAnalyticSolution<analytic_solution_tag,
125  analytic_solution_fields>,
127  system,
129  typename system::variables_tag>,
132  Metavariables>,
133  typename linear_solver::initialize_element,
136  // Initialization is done. Avoid introducing an extra phase by
137  // advancing the linear solver to the first step here.
138  typename linear_solver::prepare_step,
140 
141  // Specify all parallel components that will execute actions at some point.
142  using component_list = tmpl::append<
143  tmpl::list<elliptic::DgElementArray<
144  Metavariables,
145  tmpl::list<
146  Parallel::PhaseActions<Phase, Phase::Initialization,
147  initialization_actions>,
148 
150  Phase, Phase::RegisterWithObserver,
153  LinearSolver::Tags::IterationId,
154  element_observation_type>>,
156 
158  Phase, Phase::Solve,
159  tmpl::list<
165  typename system::variables_tag>>,
166  elliptic::dg::Actions::
167  ImposeHomogeneousDirichletBoundaryConditions<
168  Metavariables>,
171  typename linear_solver::perform_step,
172  typename linear_solver::prepare_step>>>>>,
173  typename linear_solver::component_list,
174  tmpl::list<observers::Observer<Metavariables>,
176 
177  // Specify the transitions between phases.
178  static Phase determine_next_phase(
179  const Phase& current_phase,
180  const Parallel::CProxy_ConstGlobalCache<
181  Metavariables>& /*cache_proxy*/) noexcept {
182  switch (current_phase) {
183  case Phase::Initialization:
184  return Phase::RegisterWithObserver;
185  case Phase::RegisterWithObserver:
186  return Phase::Solve;
187  case Phase::Solve:
188  return Phase::Exit;
189  case Phase::Exit:
190  ERROR(
191  "Should never call determine_next_phase with the current phase "
192  "being 'Exit'");
193  default:
194  ERROR(
195  "Unknown type of phase. Did you static_cast<Phase> an integral "
196  "value?");
197  }
198  }
199 };
200 
201 static const std::vector<void (*)()> charm_init_node_funcs{
202  &setup_error_handling, &domain::creators::register_derived_with_charm,
207 static const std::vector<void (*)()> charm_init_proc_funcs{
209 /// \endcond
The global cache tag for the numerical flux.
Definition: Tags.hpp:88
A GMRES solver for nonsymmetric linear systems of equations .
Definition: Gmres.hpp:65
Functions for serializing factory-created classes.
Adds boundary contributions to the sources.
Definition: ImposeInhomogeneousBoundaryConditionsOnSource.hpp:72
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
Definition: ConservativeSystem.hpp:50
Initializes the DataBox tags related to the elliptic system.
Definition: InitializeSystem.hpp:73
T system(T... args)
Defines class Poisson::FirstOrderSystem.
Contains the events and triggers.
Definition: Tags.hpp:57
Registers itself with the local observer parallel component so the observer knows to expect data from...
Definition: Actions.hpp:232
Initialize mortars between elements for exchanging fluxes.
Definition: InitializeMortars.hpp:79
Apply the function Mutator::apply to the DataBox.
Definition: MutateApply.hpp:40
Receive boundary data needed for fluxes from neighbors.
Definition: FluxCommunication.hpp:60
Defines classes and functions for making classes creatable from input files.
Defines actions SendDataForFluxes and ReceiveDataForFluxes.
The nodegroup parallel component that is responsible for writing data to disk.
Definition: ObserverComponent.hpp:48
Initialize items related to the basic structure of the element.
Definition: InitializeDomain.hpp:63
Place the analytic solution of the system fields in the DataBox.
Definition: InitializeAnalyticSolution.hpp:45
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:30
Holds an IterationId that identifies a step in the linear solver algorithm.
Definition: Tags.hpp:79
void register_derived_classes_with_charm() noexcept
Register derived classes of the Base class.
Definition: RegisterDerivedClassesWithCharm.hpp:31
void enable_floating_point_exceptions()
After a call to this function, the code will terminate with a floating point exception on overflow...
Definition: FloatingPointExceptions.cpp:27
tmpl::list< Tags... > slice_tags_to_exterior
Tags that are to be sliced to the exterior side of the faces of the element
Definition: InitializeInterfaces.hpp:37
Terminate the algorithm if the linear solver has converged.
Definition: TerminateIfConverged.hpp:32
Definition: RemoveOptionsAndTerminatePhase.hpp:27
A template for defining a registrar.
Definition: Registration.hpp:42
The parallel component responsible for managing the DG elements that compose the computational domain...
Definition: DgElementArray.hpp:37
Send local boundary data needed for fluxes to neighbors.
Definition: FluxCommunication.hpp:206
The internal penalty flux for first-order elliptic equations.
Definition: FirstOrderInternalPenalty.hpp:134
Functions to enable/disable termination on floating point exceptions.
Run the events and triggers.
Definition: RunEventsAndTriggers.hpp:27
Terminate the algorithm to proceed to the next phase.
Definition: TerminatePhase.hpp:26
Mutating DataBox invokable to compute the bulk contribution to the operator represented by the Operat...
Definition: FirstOrderOperator.hpp:201
Wraps the template metaprogramming library used (brigand)
Base class for something that can happen during a simulation (such as an observation).
Definition: Event.hpp:30
List of all the actions to be executed in the specified phase.
Definition: PhaseDependentActionList.hpp:16
Passed to RegisterWithObservers action to register observer event.
Definition: RegisterObservers.hpp:15
Defines tags related to domain quantities.
The analytic solution, with the type of the analytic solution set as the template parameter...
Definition: Tags.hpp:51
Holds an object that computes the principal part of the elliptic PDEs.
Definition: Tags.hpp:27
tmpl::list< Tags... > slice_tags_to_face
Tags that are to be sliced to the faces of the element
Definition: InitializeInterfaces.hpp:32
Compute element boundary contributions to the temporal step of the variables.
Definition: ApplyFluxes.hpp:70
Definition: ObserveFields.hpp:62
The linear operator applied to the data in Tag
Definition: Tags.hpp:66
Defines class template ConstGlobalCache.
Defines DataBox tags for the linear solver.
Base class for checking whether to run an Event.
Definition: Trigger.hpp:34
Initialize items related to the interfaces between Elements and on external boundaries.
Definition: InitializeInterfaces.hpp:141
Initialize DataBox tags related to discontinuous Galerkin fluxes.
Definition: InitializeFluxes.hpp:45