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/DgElementArray.hpp"
15 #include "Elliptic/DiscontinuousGalerkin/ImposeBoundaryConditions.hpp"
16 #include "Elliptic/DiscontinuousGalerkin/ImposeInhomogeneousBoundaryConditionsOnSource.hpp"
17 #include "Elliptic/DiscontinuousGalerkin/InitializeFirstOrderOperator.hpp"
18 #include "Elliptic/DiscontinuousGalerkin/NumericalFluxes/FirstOrderInternalPenalty.hpp"
19 #include "Elliptic/FirstOrderOperator.hpp"
20 #include "Elliptic/Systems/Elasticity/FirstOrderSystem.hpp"
21 #include "Elliptic/Systems/Elasticity/Tags.hpp"
22 #include "Elliptic/Tags.hpp"
23 #include "Elliptic/Triggers/EveryNIterations.hpp"
24 #include "IO/Observer/Actions/RegisterEvents.hpp"
25 #include "IO/Observer/Helpers.hpp"
26 #include "IO/Observer/ObserverComponent.hpp"
27 #include "NumericalAlgorithms/Convergence/Tags.hpp"
28 #include "NumericalAlgorithms/DiscontinuousGalerkin/BoundarySchemes/FirstOrder/FirstOrderScheme.hpp"
29 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
30 #include "Options/Options.hpp"
31 #include "Parallel/Actions/SetupDataBox.hpp"
32 #include "Parallel/Actions/TerminatePhase.hpp"
33 #include "Parallel/GlobalCache.hpp"
34 #include "Parallel/InitializationFunctions.hpp"
35 #include "Parallel/PhaseDependentActionList.hpp"
36 #include "Parallel/Reduction.hpp"
37 #include "Parallel/RegisterDerivedClassesWithCharm.hpp"
38 #include "ParallelAlgorithms/Actions/MutateApply.hpp"
39 #include "ParallelAlgorithms/DiscontinuousGalerkin/CollectDataForFluxes.hpp"
40 #include "ParallelAlgorithms/DiscontinuousGalerkin/FluxCommunication.hpp"
41 #include "ParallelAlgorithms/DiscontinuousGalerkin/InitializeDomain.hpp"
42 #include "ParallelAlgorithms/DiscontinuousGalerkin/InitializeInterfaces.hpp"
43 #include "ParallelAlgorithms/DiscontinuousGalerkin/InitializeMortars.hpp"
44 #include "ParallelAlgorithms/Events/ObserveErrorNorms.hpp"
45 #include "ParallelAlgorithms/Events/ObserveFields.hpp"
46 #include "ParallelAlgorithms/Events/ObserveVolumeIntegrals.hpp"
47 #include "ParallelAlgorithms/EventsAndTriggers/Actions/RunEventsAndTriggers.hpp"
48 #include "ParallelAlgorithms/Initialization/Actions/AddComputeTags.hpp"
49 #include "ParallelAlgorithms/Initialization/Actions/RemoveOptionsAndTerminatePhase.hpp"
50 #include "ParallelAlgorithms/LinearSolver/Gmres/Gmres.hpp"
51 #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
52 #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
53 #include "PointwiseFunctions/AnalyticData/Elasticity/AnalyticData.hpp"
54 #include "PointwiseFunctions/AnalyticSolutions/Elasticity/AnalyticSolution.hpp"
55 #include "PointwiseFunctions/AnalyticSolutions/Elasticity/BentBeam.hpp"
56 #include "PointwiseFunctions/AnalyticSolutions/Elasticity/HalfSpaceMirror.hpp"
57 #include "PointwiseFunctions/AnalyticSolutions/Elasticity/Zero.hpp"
58 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
59 #include "PointwiseFunctions/Elasticity/PotentialEnergy.hpp"
60 #include "Utilities/Blas.hpp"
61 #include "Utilities/ErrorHandling/FloatingPointExceptions.hpp"
62 #include "Utilities/Functional.hpp"
63 #include "Utilities/TMPL.hpp"
64 
65 namespace SolveElasticity::OptionTags {
66 struct LinearSolverGroup {
67  static std::string name() noexcept { return "LinearSolver"; }
68  static constexpr Options::String help =
69  "The iterative Krylov-subspace linear solver";
70 };
71 struct GmresGroup {
72  static std::string name() noexcept { return "GMRES"; }
73  static constexpr Options::String help = "Options for the GMRES linear solver";
74  using group = LinearSolverGroup;
75 };
76 } // namespace SolveElasticity::OptionTags
77 
78 /// \cond
79 template <size_t Dim>
80 struct Metavariables {
81  static constexpr size_t volume_dim = Dim;
82  using system = Elasticity::FirstOrderSystem<Dim>;
83 
84  // List the possible backgrounds, i.e. the variable-independent part of the
85  // equations that define the problem to solve (along with the boundary
86  // conditions). We currently only have analytic solutions implemented, but
87  // will add non-solution backgrounds ASAP.
88  using analytic_solution_registrars = tmpl::flatten<tmpl::list<
89  tmpl::conditional_t<Dim == 2, Elasticity::Solutions::Registrars::BentBeam,
90  tmpl::list<>>,
91  tmpl::conditional_t<Dim == 3,
92  Elasticity::Solutions::Registrars::HalfSpaceMirror,
93  tmpl::list<>>>>;
94  using background_registrars = analytic_solution_registrars;
95  using background_tag = elliptic::Tags::Background<
96  Elasticity::AnalyticData::AnalyticData<Dim, background_registrars>>;
97 
98  // List the possible initial guesses. We currently only support the trivial
99  // "zero" initial guess. This will be generalized ASAP.
100  using initial_guess_registrars =
101  tmpl::list<Elasticity::Solutions::Registrars::Zero<Dim>>;
102  using initial_guess_tag = elliptic::Tags::InitialGuess<
103  ::AnalyticData<Dim, initial_guess_registrars>>;
104 
105  static constexpr Options::String help{
106  "Find the solution to a linear elasticity problem."};
107 
108  using fluxes_computer_tag =
109  elliptic::Tags::FluxesComputer<typename system::fluxes_computer>;
110 
111  // Only Dirichlet boundary conditions are currently supported, and they are
112  // are all imposed by analytic solutions right now. This will be generalized
113  // ASAP and this alias deleted.
114  using analytic_solution_tag = background_tag;
115 
116  // The linear solver algorithm. We must use GMRES since the operator is
117  // not positive-definite for the first-order system.
118  using linear_solver =
119  LinearSolver::gmres::Gmres<Metavariables, typename system::fields_tag,
120  SolveElasticity::OptionTags::GmresGroup,
121  false>;
122  using linear_solver_iteration_id =
123  Convergence::Tags::IterationId<typename linear_solver::options_group>;
124  // For the GMRES linear solver we need to apply the DG operator to its
125  // internal "operand" in every iteration of the algorithm.
126  using linear_operand_tag = db::add_tag_prefix<LinearSolver::Tags::Operand,
127  typename system::fields_tag>;
128  using primal_variables = db::wrap_tags_in<LinearSolver::Tags::Operand,
129  typename system::primal_fields>;
130  using auxiliary_variables =
131  db::wrap_tags_in<LinearSolver::Tags::Operand,
132  typename system::auxiliary_fields>;
133 
134  // Parse numerical flux parameters from the input file to store in the cache.
135  using normal_dot_numerical_flux = Tags::NumericalFlux<
136  elliptic::dg::NumericalFluxes::FirstOrderInternalPenalty<
137  volume_dim, fluxes_computer_tag, primal_variables,
138  auxiliary_variables>>;
139  // Specify the DG boundary scheme. We use the strong first-order scheme here
140  // that only requires us to compute normals dotted into the first-order
141  // fluxes.
142  using boundary_scheme = dg::FirstOrderScheme::FirstOrderScheme<
143  volume_dim, linear_operand_tag,
144  db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo,
145  linear_operand_tag>,
146  normal_dot_numerical_flux, linear_solver_iteration_id>;
147 
148  // Collect events and triggers
149  // (public for use by the Charm++ registration code)
150  using observe_fields = typename system::fields_tag::tags_list;
151  using analytic_solution_fields = observe_fields;
152  using events = tmpl::list<
153  dg::Events::Registrars::ObserveFields<
154  volume_dim, linear_solver_iteration_id, observe_fields,
155  analytic_solution_fields>,
156  dg::Events::Registrars::ObserveErrorNorms<linear_solver_iteration_id,
157  analytic_solution_fields>,
158  dg::Events::Registrars::ObserveVolumeIntegrals<
159  volume_dim, linear_solver_iteration_id,
160  tmpl::list<Elasticity::Tags::PotentialEnergyDensity<volume_dim>>>>;
161  using triggers = tmpl::list<elliptic::Triggers::Registrars::EveryNIterations<
162  linear_solver_iteration_id>>;
163 
164  // Collect all items to store in the cache.
165  using const_global_cache_tags =
166  tmpl::list<background_tag, initial_guess_tag, fluxes_computer_tag,
167  normal_dot_numerical_flux,
168  Tags::EventsAndTriggers<events, triggers>>;
169 
170  // Collect all reduction tags for observers
171  using observed_reduction_data_tags =
172  observers::collect_reduction_data_tags<tmpl::flatten<tmpl::list<
173  typename Event<events>::creatable_classes, linear_solver>>>;
174 
175  // Specify all global synchronization points.
176  enum class Phase { Initialization, RegisterWithObserver, Solve, Exit };
177 
178  using initialization_actions = tmpl::list<
179  Actions::SetupDataBox, dg::Actions::InitializeDomain<volume_dim>,
180  dg::Actions::InitializeInterfaces<
181  system, dg::Initialization::slice_tags_to_face<>,
182  dg::Initialization::slice_tags_to_exterior<>,
183  dg::Initialization::face_compute_tags<
184  domain::Tags::BoundaryCoordinates<volume_dim>>,
185  dg::Initialization::exterior_compute_tags<>, false, false>,
186  typename linear_solver::initialize_element,
187  elliptic::Actions::InitializeFields<system, initial_guess_tag>,
188  elliptic::Actions::InitializeFixedSources<system, background_tag>,
189  Initialization::Actions::AddComputeTags<tmpl::list<
190  Elasticity::Tags::ConstitutiveRelationReference<volume_dim,
191  background_tag>,
192  Elasticity::Tags::PotentialEnergyDensityCompute<volume_dim>>>,
193  elliptic::Actions::InitializeOptionalAnalyticSolution<
194  background_tag, analytic_solution_fields,
195  Elasticity::Solutions::AnalyticSolution<Dim, background_registrars>>,
196  elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource<
197  Metavariables>,
198  dg::Actions::InitializeMortars<boundary_scheme>,
199  elliptic::dg::Actions::InitializeFirstOrderOperator<
200  volume_dim, typename system::fluxes_computer,
201  typename system::sources_computer, linear_operand_tag,
202  primal_variables, auxiliary_variables>,
203  Initialization::Actions::RemoveOptionsAndTerminatePhase>;
204 
205  using build_linear_operator_actions = tmpl::list<
206  dg::Actions::CollectDataForFluxes<
207  boundary_scheme, domain::Tags::InternalDirections<volume_dim>>,
208  dg::Actions::SendDataForFluxes<boundary_scheme>,
209  Actions::MutateApply<elliptic::FirstOrderOperator<
210  volume_dim, LinearSolver::Tags::OperatorAppliedTo,
211  linear_operand_tag>>,
212  elliptic::dg::Actions::ImposeHomogeneousDirichletBoundaryConditions<
213  linear_operand_tag, primal_variables>,
214  dg::Actions::CollectDataForFluxes<
215  boundary_scheme,
216  domain::Tags::BoundaryDirectionsInterior<volume_dim>>,
217  dg::Actions::ReceiveDataForFluxes<boundary_scheme>,
218  Actions::MutateApply<boundary_scheme>>;
219 
220  using register_actions =
221  tmpl::list<observers::Actions::RegisterEventsWithObservers,
222  Parallel::Actions::TerminatePhase>;
223 
224  using solve_actions = tmpl::list<
225  typename linear_solver::template solve<tmpl::list<
226  Actions::RunEventsAndTriggers, build_linear_operator_actions>>,
227  Actions::RunEventsAndTriggers, Parallel::Actions::TerminatePhase>;
228 
229  using dg_element_array = elliptic::DgElementArray<
230  Metavariables,
231  tmpl::list<Parallel::PhaseActions<Phase, Phase::Initialization,
232  initialization_actions>,
233  Parallel::PhaseActions<Phase, Phase::RegisterWithObserver,
234  register_actions>,
235  Parallel::PhaseActions<Phase, Phase::Solve, solve_actions>>>;
236 
237  // Specify all parallel components that will execute actions at some point.
238  using component_list = tmpl::flatten<
239  tmpl::list<dg_element_array, typename linear_solver::component_list,
240  observers::Observer<Metavariables>,
241  observers::ObserverWriter<Metavariables>>>;
242 
243  // Specify the transitions between phases.
244  static Phase determine_next_phase(
245  const Phase& current_phase,
246  const Parallel::CProxy_GlobalCache<
247  Metavariables>& /*cache_proxy*/) noexcept {
248  switch (current_phase) {
249  case Phase::Initialization:
250  return Phase::RegisterWithObserver;
251  case Phase::RegisterWithObserver:
252  return Phase::Solve;
253  case Phase::Solve:
254  return Phase::Exit;
255  case Phase::Exit:
256  ERROR(
257  "Should never call determine_next_phase with the current phase "
258  "being 'Exit'");
259  default:
260  ERROR(
261  "Unknown type of phase. Did you static_cast<Phase> an integral "
262  "value?");
263  }
264  }
265 };
266 
267 static const std::vector<void (*)()> charm_init_node_funcs{
268  &setup_error_handling,
269  &disable_openblas_multithreading,
270  &domain::creators::register_derived_with_charm,
271  &Parallel::register_derived_classes_with_charm<
272  metavariables::background_tag::type::element_type>,
273  &Parallel::register_derived_classes_with_charm<
274  metavariables::initial_guess_tag::type::element_type>,
275  &Parallel::register_derived_classes_with_charm<
276  metavariables::system::boundary_conditions_base>,
277  &Parallel::register_derived_classes_with_charm<
278  Event<metavariables::events>>,
279  &Parallel::register_derived_classes_with_charm<
280  Trigger<metavariables::triggers>>};
281 static const std::vector<void (*)()> charm_init_proc_funcs{
282  &enable_floating_point_exceptions};
283 /// \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:97
elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource
Adds boundary contributions to the sources.
Definition: ImposeInhomogeneousBoundaryConditionsOnSource.hpp:71
LinearSolver::Tags::OperatorAppliedTo
The linear operator applied to the data in Tag
Definition: Tags.hpp:54
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
dg::Actions::SendDataForFluxes
Send local boundary data needed for fluxes to neighbors.
Definition: FluxCommunication.hpp:78
Initialization::Actions::RemoveOptionsAndTerminatePhase
Definition: RemoveOptionsAndTerminatePhase.hpp:27
dg::Initialization::face_compute_tags
tmpl::list< Tags... > face_compute_tags
Definition: InitializeInterfaces.hpp:41
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:84
elliptic::Tags::Background
The variable-independent part of the elliptic equations, e.g. the fixed-sources in a Poisson equatio...
Definition: Tags.hpp:69
dg::Initialization::slice_tags_to_face
tmpl::list< Tags... > slice_tags_to_face
Definition: InitializeInterfaces.hpp:31
disable_openblas_multithreading
void disable_openblas_multithreading() noexcept
Disable OpenBLAS multithreading since it conflicts with Charm++ parallelism.
Definition: Blas.cpp:15
domain::Tags::BoundaryDirectionsInterior
Definition: Tags.hpp:296
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:66
elliptic::FirstOrderOperator
Mutating DataBox invokable to compute the bulk contribution to the operator represented by the Operat...
Definition: FirstOrderOperator.hpp:224
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:51
dg::Actions::InitializeInterfaces
Initialize items related to the interfaces between Elements and on external boundaries.
Definition: InitializeInterfaces.hpp:120
Event
Definition: Event.hpp:30
SolveElasticity::OptionTags::GmresGroup
Definition: SolveElasticity.hpp:71
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::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
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:52
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::Actions::InitializeMortars
Initialize mortars between elements for exchanging fluxes.
Definition: InitializeMortars.hpp:78
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:68
elliptic::dg::NumericalFluxes::FirstOrderInternalPenalty
The internal penalty flux for first-order elliptic equations.
Definition: FirstOrderInternalPenalty.hpp:147
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
dg::Initialization::exterior_compute_tags
tmpl::list< Tags... > exterior_compute_tags
Definition: InitializeInterfaces.hpp:46
dg::Events::Registrars::ObserveVolumeIntegrals
Definition: ObserveVolumeIntegrals.hpp:51
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
dg::FirstOrderScheme::FirstOrderScheme
Boundary contributions for a first-order DG scheme.
Definition: FirstOrderScheme.hpp:66
domain::Tags::InternalDirections
Definition: Tags.hpp:269
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:54
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:31
Blas.hpp
elliptic::dg::Actions::ImposeHomogeneousDirichletBoundaryConditions
Set field data on external boundaries so that they represent homogeneous (zero) Dirichlet boundary co...
Definition: ImposeBoundaryConditions.hpp:91
elliptic::dg::Actions::InitializeFirstOrderOperator
Initialize DataBox tags for building the first-order elliptic DG operator.
Definition: InitializeFirstOrderOperator.hpp:49
dg::Actions::ReceiveDataForFluxes
Receive boundary data needed for fluxes from neighbors.
Definition: FluxCommunication.hpp:175
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:53
TMPL.hpp
domain::Tags::BoundaryCoordinates
Definition: InterfaceComputeTags.hpp:186
dg::Initialization::slice_tags_to_exterior
tmpl::list< Tags... > slice_tags_to_exterior
Definition: InitializeInterfaces.hpp:36
© Copyright 2017 - 2021 SXS Collaboration, MIT License