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/Factory1D.hpp"
10 #include "Domain/Creators/Factory2D.hpp"
11 #include "Domain/Creators/Factory3D.hpp"
12 #include "Domain/Creators/RegisterDerivedWithCharm.hpp"
13 #include "Domain/Tags.hpp"
14 #include "Elliptic/Actions/InitializeAnalyticSolution.hpp"
15 #include "Elliptic/Actions/InitializeFields.hpp"
16 #include "Elliptic/Actions/InitializeFixedSources.hpp"
17 #include "Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp"
18 #include "Elliptic/DiscontinuousGalerkin/Actions/InitializeDomain.hpp"
19 #include "Elliptic/DiscontinuousGalerkin/DgElementArray.hpp"
20 #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/InitializeSubdomain.hpp"
21 #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp"
22 #include "Elliptic/Systems/Elasticity/FirstOrderSystem.hpp"
23 #include "Elliptic/Systems/Elasticity/Tags.hpp"
24 #include "Elliptic/Tags.hpp"
25 #include "Elliptic/Triggers/EveryNIterations.hpp"
26 #include "IO/Observer/Actions/RegisterEvents.hpp"
27 #include "IO/Observer/Helpers.hpp"
28 #include "IO/Observer/ObserverComponent.hpp"
29 #include "NumericalAlgorithms/Convergence/Tags.hpp"
30 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
31 #include "Options/Options.hpp"
32 #include "Options/Protocols/FactoryCreation.hpp"
33 #include "Parallel/Actions/SetupDataBox.hpp"
34 #include "Parallel/Actions/TerminatePhase.hpp"
35 #include "Parallel/GlobalCache.hpp"
36 #include "Parallel/InitializationFunctions.hpp"
37 #include "Parallel/PhaseDependentActionList.hpp"
38 #include "Parallel/Reduction.hpp"
40 #include "ParallelAlgorithms/Actions/MutateApply.hpp"
41 #include "ParallelAlgorithms/Events/Factory.hpp"
42 #include "ParallelAlgorithms/Events/ObserveVolumeIntegrals.hpp"
43 #include "ParallelAlgorithms/EventsAndTriggers/Actions/RunEventsAndTriggers.hpp"
44 #include "ParallelAlgorithms/EventsAndTriggers/Completion.hpp"
45 #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
46 #include "ParallelAlgorithms/EventsAndTriggers/LogicalTriggers.hpp"
48 #include "ParallelAlgorithms/EventsAndTriggers/Trigger.hpp"
49 #include "ParallelAlgorithms/Initialization/Actions/AddComputeTags.hpp"
50 #include "ParallelAlgorithms/Initialization/Actions/RemoveOptionsAndTerminatePhase.hpp"
51 #include "ParallelAlgorithms/LinearSolver/Gmres/Gmres.hpp"
52 #include "ParallelAlgorithms/LinearSolver/Schwarz/Schwarz.hpp"
54 #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
55 #include "PointwiseFunctions/AnalyticData/Elasticity/AnalyticData.hpp"
56 #include "PointwiseFunctions/AnalyticSolutions/Elasticity/AnalyticSolution.hpp"
57 #include "PointwiseFunctions/AnalyticSolutions/Elasticity/BentBeam.hpp"
58 #include "PointwiseFunctions/AnalyticSolutions/Elasticity/HalfSpaceMirror.hpp"
59 #include "PointwiseFunctions/AnalyticSolutions/Elasticity/Zero.hpp"
60 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
61 #include "PointwiseFunctions/Elasticity/PotentialEnergy.hpp"
62 #include "PointwiseFunctions/Elasticity/Strain.hpp"
63 #include "Utilities/Blas.hpp"
65 #include "Utilities/Functional.hpp"
66 #include "Utilities/MemoryHelpers.hpp"
67 #include "Utilities/ProtocolHelpers.hpp"
68 #include "Utilities/TMPL.hpp"
69 
70 /// \cond
71 namespace PUP {
72 class er;
73 } // namespace PUP
74 /// \endcond
75 
76 namespace SolveElasticity::OptionTags {
78  static std::string name() noexcept { return "LinearSolver"; }
79  static constexpr Options::String help =
80  "The iterative Krylov-subspace linear solver";
81 };
82 struct GmresGroup {
83  static std::string name() noexcept { return "GMRES"; }
84  static constexpr Options::String help = "Options for the GMRES linear solver";
85  using group = LinearSolverGroup;
86 };
88  static std::string name() noexcept { return "SchwarzSmoother"; }
89  static constexpr Options::String help = "Options for the Schwarz smoother";
90  using group = LinearSolverGroup;
91 };
92 } // namespace SolveElasticity::OptionTags
93 
94 /// \cond
95 template <size_t Dim>
96 struct Metavariables {
97  static constexpr size_t volume_dim = Dim;
98  using system = Elasticity::FirstOrderSystem<Dim>;
99 
100  // List the possible backgrounds, i.e. the variable-independent part of the
101  // equations that define the problem to solve (along with the boundary
102  // conditions). We currently only have analytic solutions implemented, but
103  // will add non-solution backgrounds ASAP.
104  using analytic_solution_registrars = tmpl::flatten<tmpl::list<
105  tmpl::conditional_t<Dim == 2, Elasticity::Solutions::Registrars::BentBeam,
106  tmpl::list<>>,
107  tmpl::conditional_t<Dim == 3,
108  Elasticity::Solutions::Registrars::HalfSpaceMirror,
109  tmpl::list<>>>>;
110  using background_registrars = analytic_solution_registrars;
111  using background_tag = elliptic::Tags::Background<
113 
114  // List the possible initial guesses
115  using initial_guess_registrars =
116  tmpl::append<tmpl::list<Elasticity::Solutions::Registrars::Zero<Dim>>,
117  analytic_solution_registrars>;
118  using initial_guess_tag = elliptic::Tags::InitialGuess<
120 
121  static constexpr Options::String help{
122  "Find the solution to a linear elasticity problem."};
123 
124  // These are the fields we solve for
126  // These are the fixed sources, i.e. the RHS of the equations
127  using fixed_sources_tag = db::add_tag_prefix<::Tags::FixedSource, fields_tag>;
128  // This is the linear operator applied to the fields. We'll only use it to
129  // apply the operator to the initial guess, so an optimization would be to
130  // re-use the `operator_applied_to_vars_tag` below. This optimization needs a
131  // few minor changes to the parallel linear solver algorithm.
132  using operator_applied_to_fields_tag =
134 
135  // The linear solver algorithm. We must use GMRES since the operator is
136  // not guaranteed to be symmetric. It can be made symmetric by multiplying by
137  // the DG mass matrix.
138  using linear_solver =
139  LinearSolver::gmres::Gmres<Metavariables, fields_tag,
141  using linear_solver_iteration_id =
143  // Precondition each linear solver iteration with a number of Schwarz
144  // smoothing steps
145  using subdomain_operator =
148  tmpl::list<Elasticity::Tags::ConstitutiveRelation<Dim>>>;
149  using schwarz_smoother = LinearSolver::Schwarz::Schwarz<
150  typename linear_solver::operand_tag,
152  typename linear_solver::preconditioner_source_tag>;
153  // For the GMRES linear solver we need to apply the DG operator to its
154  // internal "operand" in every iteration of the algorithm.
155  using vars_tag = typename linear_solver::operand_tag;
156  using operator_applied_to_vars_tag =
158  // We'll buffer the corresponding fluxes in this tag, but won't actually need
159  // to access them outside applying the operator
160  using fluxes_vars_tag =
162  typename system::primal_fluxes>>;
163 
164  using analytic_solution_fields = typename system::primal_fields;
165  using observe_fields = tmpl::append<
166  analytic_solution_fields,
167  tmpl::list<Elasticity::Tags::Strain<volume_dim>,
169 
170  // Collect all items to store in the cache.
171  using const_global_cache_tags =
172  tmpl::list<background_tag, initial_guess_tag, Tags::EventsAndTriggers>;
173 
174  struct factory_creation
175  : tt::ConformsTo<Options::protocols::FactoryCreation> {
176  using factory_classes = tmpl::map<
177  tmpl::pair<DomainCreator<volume_dim>, domain_creators<volume_dim>>,
178  tmpl::pair<Event,
179  tmpl::flatten<tmpl::list<
181  dg::Events::field_observations<
182  volume_dim, linear_solver_iteration_id,
183  observe_fields, analytic_solution_fields>,
185  volume_dim, linear_solver_iteration_id,
187  volume_dim>>>>>>,
188  tmpl::pair<Trigger,
189  tmpl::push_back<Triggers::logical_triggers,
191  linear_solver_iteration_id>>>>;
192  };
193 
194  // Collect all reduction tags for observers
195  using observed_reduction_data_tags =
196  observers::collect_reduction_data_tags<tmpl::flatten<tmpl::list<
197  tmpl::at<typename factory_creation::factory_classes, Event>,
198  linear_solver, schwarz_smoother>>>;
199 
200  // Specify all global synchronization points.
201  enum class Phase { Initialization, RegisterWithObserver, Solve, Exit };
202 
203  using initialization_actions = tmpl::list<
206  typename linear_solver::initialize_element,
207  typename schwarz_smoother::initialize_element,
212  background_tag>,
216  background_tag,
217  tmpl::append<typename system::primal_fields,
218  typename system::primal_fluxes>,
222  system, background_tag, typename schwarz_smoother::options_group>,
224  system, fixed_sources_tag>,
225  // Apply the DG operator to the initial guess
227  system, true, linear_solver_iteration_id, fields_tag, fluxes_vars_tag,
228  operator_applied_to_fields_tag, vars_tag, fluxes_vars_tag>,
230 
231  using build_linear_operator_actions = elliptic::dg::Actions::apply_operator<
232  system, true, linear_solver_iteration_id, vars_tag, fluxes_vars_tag,
233  operator_applied_to_vars_tag>;
234 
235  using register_actions =
237  typename schwarz_smoother::register_element,
239 
240  using solve_actions =
241  tmpl::list<typename linear_solver::template solve<
243  typename schwarz_smoother::template solve<
244  build_linear_operator_actions>>>,
247 
248  using dg_element_array = elliptic::DgElementArray<
249  Metavariables,
250  tmpl::list<Parallel::PhaseActions<Phase, Phase::Initialization,
251  initialization_actions>,
252  Parallel::PhaseActions<Phase, Phase::RegisterWithObserver,
253  register_actions>,
255 
256  // Specify all parallel components that will execute actions at some point.
257  using component_list = tmpl::flatten<
258  tmpl::list<dg_element_array, typename linear_solver::component_list,
259  typename schwarz_smoother::component_list,
262 
263  // Specify the transitions between phases.
264  template <typename... Tags>
265  static Phase determine_next_phase(
266  const gsl::not_null<
267  tuples::TaggedTuple<Tags...>*> /*phase_change_decision_data*/,
268  const Phase& current_phase,
269  const Parallel::CProxy_GlobalCache<
270  Metavariables>& /*cache_proxy*/) noexcept {
271  switch (current_phase) {
272  case Phase::Initialization:
273  return Phase::RegisterWithObserver;
274  case Phase::RegisterWithObserver:
275  return Phase::Solve;
276  case Phase::Solve:
277  return Phase::Exit;
278  case Phase::Exit:
279  ERROR(
280  "Should never call determine_next_phase with the current phase "
281  "being 'Exit'");
282  default:
283  ERROR(
284  "Unknown type of phase. Did you static_cast<Phase> an integral "
285  "value?");
286  }
287  }
288 
289  // NOLINTNEXTLINE(google-runtime-references)
290  void pup(PUP::er& /*p*/) noexcept {}
291 };
292 
293 static const std::vector<void (*)()> charm_init_node_funcs{
294  &setup_error_handling, &setup_memory_allocation_failure_reporting,
296  &domain::creators::register_derived_with_charm,
298  metavariables::background_tag::type::element_type>,
300  metavariables::initial_guess_tag::type::element_type>,
302  metavariables::system::boundary_conditions_base>,
304  metavariables::schwarz_smoother::subdomain_solver>,
305  &Parallel::register_factory_classes_with_charm<metavariables>};
306 static const std::vector<void (*)()> charm_init_proc_funcs{
308 /// \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:89
std::string
observers::Actions::RegisterEventsWithObservers
Registers this element of a parallel component with the local Observer parallel component for each tr...
Definition: RegisterEvents.hpp:114
Tags.hpp
elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource
For linear systems, impose inhomogeneous boundary conditions as contributions to the fixed sources (i...
Definition: ApplyOperator.hpp:500
GlobalCache.hpp
Options.hpp
LinearSolver::Tags::Operand
The operand that the local linear operator is applied to.
Definition: Tags.hpp:40
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)
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
elliptic::dg::subdomain_operator::Actions::InitializeSubdomain
Initialize the geometry for the DG subdomain operator.
Definition: InitializeSubdomain.hpp:117
Initialization::Actions::RemoveOptionsAndTerminatePhase
Definition: RemoveOptionsAndTerminatePhase.hpp:27
Parallel::Actions::TerminatePhase
Terminate the algorithm to proceed to the next phase.
Definition: TerminatePhase.hpp:26
dg::Events::ObserveVolumeIntegrals
Definition: ObserveVolumeIntegrals.hpp:43
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
Elasticity::Tags::ConstitutiveRelationReference
Reference the constitutive relation provided by the ProviderTag
Definition: Tags.hpp:31
SolveElasticity::OptionTags::LinearSolverGroup
Definition: SolveElasticity.hpp:77
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()
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
Event
Definition: Event.hpp:19
LinearSolver::Schwarz::Schwarz
An additive Schwarz subdomain solver for linear systems of equations .
Definition: Schwarz.hpp:147
SolveElasticity::OptionTags::GmresGroup
Definition: SolveElasticity.hpp:82
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:488
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: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:450
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
Elasticity::Tags::PotentialEnergyDensity
The energy density stored in the deformation of the elastic material.
Definition: Tags.hpp:75
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
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: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
elliptic::dg::subdomain_operator::SubdomainOperator
The elliptic DG operator on an element-centered subdomain.
Definition: SubdomainOperator.hpp:99
SolveElasticity::OptionTags::SchwarzSmootherGroup
Definition: SolveElasticity.hpp:87
tt::ConformsTo
Indicate a class conforms to the Protocol.
Definition: ProtocolHelpers.hpp:22
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:54
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13