SolveElasticity.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/DgElementArray.hpp"
16 #include "Elliptic/Systems/Elasticity/FirstOrderSystem.hpp"
17 #include "Elliptic/Systems/Elasticity/Tags.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 "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
25 #include "Options/Options.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/Actions/MutateApply.hpp"
34 #include "ParallelAlgorithms/DiscontinuousGalerkin/InitializeDomain.hpp"
35 #include "ParallelAlgorithms/Events/ObserveErrorNorms.hpp"
36 #include "ParallelAlgorithms/Events/ObserveFields.hpp"
37 #include "ParallelAlgorithms/Events/ObserveVolumeIntegrals.hpp"
38 #include "ParallelAlgorithms/EventsAndTriggers/Actions/RunEventsAndTriggers.hpp"
39 #include "ParallelAlgorithms/Initialization/Actions/AddComputeTags.hpp"
40 #include "ParallelAlgorithms/Initialization/Actions/RemoveOptionsAndTerminatePhase.hpp"
41 #include "ParallelAlgorithms/LinearSolver/Gmres/Gmres.hpp"
43 #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
44 #include "PointwiseFunctions/AnalyticData/Elasticity/AnalyticData.hpp"
45 #include "PointwiseFunctions/AnalyticSolutions/Elasticity/AnalyticSolution.hpp"
46 #include "PointwiseFunctions/AnalyticSolutions/Elasticity/BentBeam.hpp"
47 #include "PointwiseFunctions/AnalyticSolutions/Elasticity/HalfSpaceMirror.hpp"
48 #include "PointwiseFunctions/AnalyticSolutions/Elasticity/Zero.hpp"
49 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
50 #include "PointwiseFunctions/Elasticity/PotentialEnergy.hpp"
51 #include "PointwiseFunctions/Elasticity/Strain.hpp"
52 #include "Utilities/Blas.hpp"
54 #include "Utilities/Functional.hpp"
55 #include "Utilities/TMPL.hpp"
56 
57 /// \cond
58 namespace PUP {
59 class er;
60 } // namespace PUP
61 /// \endcond
62 
63 namespace SolveElasticity::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 SolveElasticity::OptionTags
75 
76 /// \cond
77 template <size_t Dim>
78 struct Metavariables {
79  static constexpr size_t volume_dim = Dim;
80  using system = Elasticity::FirstOrderSystem<Dim>;
81 
82  // List the possible backgrounds, i.e. the variable-independent part of the
83  // equations that define the problem to solve (along with the boundary
84  // conditions). We currently only have analytic solutions implemented, but
85  // will add non-solution backgrounds ASAP.
86  using analytic_solution_registrars = tmpl::flatten<tmpl::list<
87  tmpl::conditional_t<Dim == 2, Elasticity::Solutions::Registrars::BentBeam,
88  tmpl::list<>>,
89  tmpl::conditional_t<Dim == 3,
90  Elasticity::Solutions::Registrars::HalfSpaceMirror,
91  tmpl::list<>>>>;
92  using background_registrars = analytic_solution_registrars;
93  using background_tag = elliptic::Tags::Background<
95 
96  // List the possible initial guesses
97  using initial_guess_registrars =
98  tmpl::append<tmpl::list<Elasticity::Solutions::Registrars::Zero<Dim>>,
99  analytic_solution_registrars>;
100  using initial_guess_tag = elliptic::Tags::InitialGuess<
102 
103  static constexpr Options::String help{
104  "Find the solution to a linear elasticity problem."};
105 
106  // These are the fields we solve for
108  // These are the fixed sources, i.e. the RHS of the equations
109  using fixed_sources_tag = db::add_tag_prefix<::Tags::FixedSource, fields_tag>;
110  // This is the linear operator applied to the fields. We'll only use it to
111  // apply the operator to the initial guess, so an optimization would be to
112  // re-use the `operator_applied_to_vars_tag` below. This optimization needs a
113  // few minor changes to the parallel linear solver algorithm.
114  using operator_applied_to_fields_tag =
116 
117  // The linear solver algorithm. We must use GMRES since the operator is
118  // not guaranteed to be symmetric. It can be made symmetric by multiplying by
119  // the DG mass matrix.
120  using linear_solver =
121  LinearSolver::gmres::Gmres<Metavariables, fields_tag,
123  false>;
124  using linear_solver_iteration_id =
126  // For the GMRES linear solver we need to apply the DG operator to its
127  // internal "operand" in every iteration of the algorithm.
128  using vars_tag = typename linear_solver::operand_tag;
129  using operator_applied_to_vars_tag =
131  // We'll buffer the corresponding fluxes in this tag, but won't actually need
132  // to access them outside applying the operator
133  using fluxes_vars_tag =
135  typename system::primal_fluxes>>;
136 
137  // Collect events and triggers
138  // (public for use by the Charm++ registration code)
139  using analytic_solution_fields = typename system::primal_fields;
140  using observe_fields = tmpl::append<
141  analytic_solution_fields,
142  tmpl::list<Elasticity::Tags::Strain<volume_dim>,
144  using events = tmpl::list<
146  volume_dim, linear_solver_iteration_id, observe_fields,
147  analytic_solution_fields>,
148  dg::Events::Registrars::ObserveErrorNorms<linear_solver_iteration_id,
149  analytic_solution_fields>,
151  volume_dim, linear_solver_iteration_id,
152  tmpl::list<Elasticity::Tags::PotentialEnergyDensity<volume_dim>>>>;
153  using triggers = tmpl::list<elliptic::Triggers::Registrars::EveryNIterations<
154  linear_solver_iteration_id>>;
155 
156  // Collect all items to store in the cache.
157  using const_global_cache_tags =
158  tmpl::list<background_tag, initial_guess_tag,
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  typename Event<events>::creatable_classes, linear_solver>>>;
165 
166  // Specify all global synchronization points.
167  enum class Phase { Initialization, RegisterWithObserver, Solve, Exit };
168 
169  using initialization_actions = tmpl::list<
171  typename linear_solver::initialize_element,
176  background_tag>,
180  background_tag,
181  tmpl::append<typename system::primal_fields,
182  typename system::primal_fluxes>,
186  system, fixed_sources_tag>,
187  // Apply the DG operator to the initial guess
189  system, true, linear_solver_iteration_id, fields_tag, fluxes_vars_tag,
190  operator_applied_to_fields_tag, vars_tag, fluxes_vars_tag>,
192 
193  using build_linear_operator_actions = elliptic::dg::Actions::apply_operator<
194  system, true, linear_solver_iteration_id, vars_tag, fluxes_vars_tag,
195  operator_applied_to_vars_tag>;
196 
197  using register_actions =
200 
201  using solve_actions = tmpl::list<
202  typename linear_solver::template solve<tmpl::list<
203  Actions::RunEventsAndTriggers, build_linear_operator_actions>>,
205 
206  using dg_element_array = elliptic::DgElementArray<
207  Metavariables,
208  tmpl::list<Parallel::PhaseActions<Phase, Phase::Initialization,
209  initialization_actions>,
210  Parallel::PhaseActions<Phase, Phase::RegisterWithObserver,
211  register_actions>,
213 
214  // Specify all parallel components that will execute actions at some point.
215  using component_list = tmpl::flatten<
216  tmpl::list<dg_element_array, typename linear_solver::component_list,
219 
220  // Specify the transitions between phases.
221  template <typename... Tags>
222  static Phase determine_next_phase(
223  const gsl::not_null<
224  tuples::TaggedTuple<Tags...>*> /*phase_change_decision_data*/,
225  const Phase& current_phase,
226  const Parallel::CProxy_GlobalCache<
227  Metavariables>& /*cache_proxy*/) noexcept {
228  switch (current_phase) {
229  case Phase::Initialization:
230  return Phase::RegisterWithObserver;
231  case Phase::RegisterWithObserver:
232  return Phase::Solve;
233  case Phase::Solve:
234  return Phase::Exit;
235  case Phase::Exit:
236  ERROR(
237  "Should never call determine_next_phase with the current phase "
238  "being 'Exit'");
239  default:
240  ERROR(
241  "Unknown type of phase. Did you static_cast<Phase> an integral "
242  "value?");
243  }
244  }
245 
246  // NOLINTNEXTLINE(google-runtime-references)
247  void pup(PUP::er& /*p*/) noexcept {}
248 };
249 
250 static const std::vector<void (*)()> charm_init_node_funcs{
251  &setup_error_handling,
253  &domain::creators::register_derived_with_charm,
255  metavariables::background_tag::type::element_type>,
257  metavariables::initial_guess_tag::type::element_type>,
259  metavariables::system::boundary_conditions_base>,
264 static const std::vector<void (*)()> charm_init_proc_funcs{
266 /// \endcond
FloatingPointExceptions.hpp
Elasticity::Solutions::AnalyticSolution
Base class for analytic solutions of the linear Elasticity equations.
Definition: AnalyticSolution.hpp:34
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
observers::Actions::RegisterEventsWithObservers
Registers this element of a parallel component with the local Observer parallel component for each tr...
Definition: RegisterEvents.hpp:105
elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource
For linear systems, impose inhomogeneous boundary conditions as contributions to the fixed sources (i...
Definition: ApplyOperator.hpp:525
GlobalCache.hpp
Options.hpp
LinearSolver::Tags::Operand
The operand that the local linear operator is applied to.
Definition: Tags.hpp:41
Initialization::Actions::AddComputeTags
Initialize the list of compute tags in ComputeTagsList
Definition: AddComputeTags.hpp:42
Tags.hpp
std::vector
std::system
T system(T... args)
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:69
disable_openblas_multithreading
void disable_openblas_multithreading() noexcept
Disable OpenBLAS multithreading since it conflicts with Charm++ parallelism.
Tags.hpp
Elasticity::Tags::ConstitutiveRelationReference
Reference the constitutive relation provided by the ProviderTag
Definition: Tags.hpp:31
Registration::Registrar
A template for defining a registrar.
Definition: Registration.hpp:42
SolveElasticity::OptionTags::LinearSolverGroup
Definition: SolveElasticity.hpp:64
enable_floating_point_exceptions
void enable_floating_point_exceptions()
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
Event
Definition: Event.hpp:30
SolveElasticity::OptionTags::GmresGroup
Definition: SolveElasticity.hpp:69
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:67
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:513
elliptic::Actions::InitializeOptionalAnalyticSolution
Definition: InitializeAnalyticSolution.hpp:87
cstddef
Elasticity::AnalyticData::AnalyticData
Base class for the background of the Elasticity system, i.e. its variable-independent quantities....
Definition: AnalyticData.hpp:25
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:34
dg::Events::Registrars::ObserveFields
Definition: ObserveFields.hpp:68
AnalyticData
Provides analytic tensor data as a function of the spatial coordinates.
Definition: AnalyticData.hpp:30
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:57
Elasticity::Tags::PotentialEnergyDensity
The energy density stored in the deformation of the elastic material.
Definition: Tags.hpp:75
dg::Events::Registrars::ObserveVolumeIntegrals
Definition: ObserveVolumeIntegrals.hpp:51
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
Elasticity::FirstOrderSystem
The linear elasticity equation formulated as a set of coupled first-order PDEs.
Definition: FirstOrderSystem.hpp:55
elliptic::Tags::InitialGuess
The initial guess for the elliptic solve.
Definition: Tags.hpp:81
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
Elasticity::Tags::StrainCompute
The symmetric strain in the elastic material.
Definition: Strain.hpp:58
Elasticity::Tags::PotentialEnergyDensityCompute
Computes the energy density stored in the deformation of the elastic material.
Definition: PotentialEnergy.hpp:52
elliptic::Actions::InitializeFixedSources
Initialize the "fixed sources" of the elliptic equations, i.e. their variable-independent source term...
Definition: InitializeFixedSources.hpp:52
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
elliptic::dg::Actions::initialize_operator
tmpl::list< detail::InitializeFacesAndMortars< System > > initialize_operator
Initialize DataBox tags for applying the DG operator.
Definition: ApplyOperator.hpp:475