Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <string>
8 :
9 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
10 : #include "DataStructures/DataBox/Prefixes.hpp"
11 : #include "DataStructures/Variables.hpp"
12 : #include "Domain/Tags/Faces.hpp"
13 : #include "Elliptic/Actions/InitializeAnalyticSolution.hpp"
14 : #include "Elliptic/Actions/InitializeFields.hpp"
15 : #include "Elliptic/Actions/InitializeFixedSources.hpp"
16 : #include "Elliptic/Actions/RunEventsAndTriggers.hpp"
17 : #include "Elliptic/Amr/Actions.hpp"
18 : #include "Elliptic/BoundaryConditions/Tags/BoundaryFields.hpp"
19 : #include "Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp"
20 : #include "Elliptic/DiscontinuousGalerkin/Actions/InitializeDomain.hpp"
21 : #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/InitializeSubdomain.hpp"
22 : #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp"
23 : #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
24 : #include "Elliptic/Protocols/FirstOrderSystem.hpp"
25 : #include "Elliptic/SubdomainPreconditioners/MinusLaplacian.hpp"
26 : #include "Elliptic/Systems/GetSourcesComputer.hpp"
27 : #include "Elliptic/Tags.hpp"
28 : #include "IO/Observer/Helpers.hpp"
29 : #include "NumericalAlgorithms/Convergence/Tags.hpp"
30 : #include "Options/String.hpp"
31 : #include "Parallel/PhaseControl/ExecutePhaseChange.hpp"
32 : #include "ParallelAlgorithms/Actions/AddComputeTags.hpp"
33 : #include "ParallelAlgorithms/Actions/InitializeItems.hpp"
34 : #include "ParallelAlgorithms/Actions/MutateApply.hpp"
35 : #include "ParallelAlgorithms/Actions/RandomizeVariables.hpp"
36 : #include "ParallelAlgorithms/Amr/Actions/Component.hpp"
37 : #include "ParallelAlgorithms/Amr/Actions/ElementsRegistration.hpp"
38 : #include "ParallelAlgorithms/Amr/Actions/Initialize.hpp"
39 : #include "ParallelAlgorithms/Amr/Projectors/DefaultInitialize.hpp"
40 : #include "ParallelAlgorithms/Amr/Projectors/Variables.hpp"
41 : #include "ParallelAlgorithms/Amr/Tags.hpp"
42 : #include "ParallelAlgorithms/LinearSolver/Actions/BuildMatrix.hpp"
43 : #include "ParallelAlgorithms/LinearSolver/Actions/MakeIdentityIfSkipped.hpp"
44 : #include "ParallelAlgorithms/LinearSolver/Gmres/Gmres.hpp"
45 : #include "ParallelAlgorithms/LinearSolver/Multigrid/Actions/RestrictFields.hpp"
46 : #include "ParallelAlgorithms/LinearSolver/Multigrid/Multigrid.hpp"
47 : #include "ParallelAlgorithms/LinearSolver/Schwarz/Actions/CommunicateOverlapFields.hpp"
48 : #include "ParallelAlgorithms/LinearSolver/Schwarz/Actions/ResetSubdomainSolver.hpp"
49 : #include "ParallelAlgorithms/LinearSolver/Schwarz/Schwarz.hpp"
50 : #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
51 : #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/NewtonRaphson.hpp"
52 : #include "ParallelAlgorithms/NonlinearSolver/Tags.hpp"
53 : #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
54 : #include "PointwiseFunctions/InitialDataUtilities/Background.hpp"
55 : #include "PointwiseFunctions/InitialDataUtilities/InitialGuess.hpp"
56 : #include "Utilities/ProtocolHelpers.hpp"
57 : #include "Utilities/TMPL.hpp"
58 :
59 : /// Items related to composing nonlinear elliptic solver executables
60 : namespace elliptic {
61 :
62 : /// Option tags for nonlinear elliptic solver executables
63 : namespace OptionTags {
64 :
65 0 : struct NonlinearSolverGroup {
66 0 : static std::string name() { return "NonlinearSolver"; }
67 0 : static constexpr Options::String help = "The iterative nonlinear solver";
68 : };
69 0 : struct NewtonRaphsonGroup {
70 0 : static std::string name() { return "NewtonRaphson"; }
71 0 : static constexpr Options::String help =
72 : "Options for the Newton-Raphson nonlinear solver";
73 0 : using group = NonlinearSolverGroup;
74 : };
75 0 : struct LinearSolverGroup {
76 0 : static std::string name() { return "LinearSolver"; }
77 0 : static constexpr Options::String help =
78 : "The iterative Krylov-subspace linear solver";
79 : };
80 0 : struct GmresGroup {
81 0 : static std::string name() { return "Gmres"; }
82 0 : static constexpr Options::String help = "Options for the GMRES linear solver";
83 0 : using group = LinearSolverGroup;
84 : };
85 0 : struct SchwarzSmootherGroup {
86 0 : static std::string name() { return "SchwarzSmoother"; }
87 0 : static constexpr Options::String help = "Options for the Schwarz smoother";
88 0 : using group = LinearSolverGroup;
89 : };
90 0 : struct MultigridGroup {
91 0 : static std::string name() { return "Multigrid"; }
92 0 : static constexpr Options::String help = "Options for the multigrid";
93 0 : using group = LinearSolverGroup;
94 : };
95 :
96 : } // namespace OptionTags
97 :
98 : /*!
99 : * \brief A complete nonlinear elliptic solver stack. Use to compose an
100 : * executable.
101 : *
102 : * The elliptic solver stack is described in detail in \cite Vu2021coj.
103 : *
104 : * Uses `Metavariables` to instantiate parallel components.
105 : */
106 : template <typename Metavariables, size_t Dim, typename System>
107 1 : struct Solver {
108 0 : static constexpr size_t volume_dim = Dim;
109 0 : using system = System;
110 : static_assert(
111 : tt::assert_conforms_to_v<system, elliptic::protocols::FirstOrderSystem>);
112 0 : static constexpr bool is_linear =
113 : std::is_same_v<elliptic::get_sources_computer<system, true>,
114 : typename system::sources_computer>;
115 :
116 0 : using background_tag =
117 : elliptic::Tags::Background<elliptic::analytic_data::Background>;
118 0 : using initial_guess_tag =
119 : elliptic::Tags::InitialGuess<elliptic::analytic_data::InitialGuess>;
120 :
121 : /// These are the fields we solve for
122 1 : using fields_tag = ::Tags::Variables<typename system::primal_fields>;
123 : /// These are the fluxes corresponding to the fields, i.e. essentially their
124 : /// first derivatives. These are background fields for the linearized sources.
125 1 : using fluxes_tag = ::Tags::Variables<typename system::primal_fluxes>;
126 : /// These are the fixed sources, i.e. the RHS of the equations
127 1 : using fixed_sources_tag = db::add_tag_prefix<::Tags::FixedSource, fields_tag>;
128 0 : using operator_applied_to_fields_tag = tmpl::conditional_t<
129 : is_linear,
130 : // This is the linear operator applied to the fields. We only use it to
131 : // apply the operator to the initial guess, so an optimization would be to
132 : // re-use the `operator_applied_to_vars_tag` below. This optimization
133 : // needs a few minor changes to the parallel linear solver algorithm.
134 : db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, fields_tag>,
135 : db::add_tag_prefix<NonlinearSolver::Tags::OperatorAppliedTo, fields_tag>>;
136 :
137 : // For AMR iterations
138 0 : using amr_iteration_id =
139 : Convergence::Tags::IterationId<::amr::OptionTags::AmrGroup>;
140 :
141 : /// The nonlinear solver algorithm
142 1 : using nonlinear_solver = NonlinearSolver::newton_raphson::NewtonRaphson<
143 : Metavariables, fields_tag, OptionTags::NewtonRaphsonGroup,
144 : fixed_sources_tag, ::amr::Tags::IsFinestGrid>;
145 0 : using nonlinear_solver_iteration_id =
146 : Convergence::Tags::IterationId<typename nonlinear_solver::options_group>;
147 :
148 : /// The linear solver algorithm. We use GMRES since the operator is not
149 : /// necessarily symmetric. Using CG here would be an optimization for
150 : /// symmetric problems.
151 1 : using linear_solver = LinearSolver::gmres::Gmres<
152 : Metavariables,
153 : tmpl::conditional_t<is_linear, fields_tag,
154 : typename nonlinear_solver::linear_solver_fields_tag>,
155 : OptionTags::GmresGroup, true,
156 : tmpl::conditional_t<is_linear, fixed_sources_tag,
157 : typename nonlinear_solver::linear_solver_source_tag>,
158 : ::amr::Tags::IsFinestGrid>;
159 0 : using linear_solver_iteration_id =
160 : Convergence::Tags::IterationId<typename linear_solver::options_group>;
161 :
162 : /// Precondition each linear solver iteration with a multigrid V-cycle
163 1 : using multigrid = LinearSolver::multigrid::Multigrid<
164 : Metavariables, volume_dim, typename linear_solver::operand_tag,
165 : OptionTags::MultigridGroup, elliptic::dg::Tags::Massive,
166 : typename linear_solver::preconditioner_source_tag>;
167 :
168 : /// Smooth each multigrid level with a number of Schwarz smoothing steps
169 1 : using subdomain_operator =
170 : elliptic::dg::subdomain_operator::SubdomainOperator<
171 : system, OptionTags::SchwarzSmootherGroup>;
172 0 : using subdomain_preconditioners = tmpl::list<
173 : elliptic::subdomain_preconditioners::Registrars::MinusLaplacian<
174 : volume_dim, OptionTags::SchwarzSmootherGroup>>;
175 0 : using schwarz_smoother = LinearSolver::Schwarz::Schwarz<
176 : typename multigrid::smooth_fields_tag, OptionTags::SchwarzSmootherGroup,
177 : subdomain_operator, subdomain_preconditioners,
178 : typename multigrid::smooth_source_tag, ::amr::Tags::GridIndex>;
179 :
180 : /// For the GMRES linear solver we need to apply the DG operator to its
181 : /// internal "operand" in every iteration of the algorithm.
182 1 : using vars_tag = typename linear_solver::operand_tag;
183 0 : using operator_applied_to_vars_tag =
184 : db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, vars_tag>;
185 : /// The correction fluxes can be stored in an arbitrary tag. We don't need to
186 : /// access them anywhere, they're just a memory buffer for the linearized
187 : /// operator.
188 1 : using fluxes_vars_tag =
189 : db::add_tag_prefix<NonlinearSolver::Tags::Correction, fluxes_tag>;
190 :
191 : /// Fields that may be observed to monitor the state of the solver
192 1 : using observe_fields = tmpl::conditional_t<
193 : is_linear,
194 : tmpl::append<
195 : typename fixed_sources_tag::tags_list,
196 : typename db::add_tag_prefix<LinearSolver::Tags::Operand,
197 : fields_tag>::tags_list,
198 : typename db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo,
199 : fields_tag>::tags_list>,
200 : tmpl::append<
201 : typename fixed_sources_tag::tags_list,
202 : typename nonlinear_solver::linear_solver_fields_tag::tags_list,
203 : typename nonlinear_solver::linear_solver_source_tag::tags_list>>;
204 :
205 : /// Collect all reduction tags for observers
206 1 : using observed_reduction_data_tags =
207 : observers::collect_reduction_data_tags<tmpl::flatten<tmpl::list<
208 : tmpl::conditional_t<is_linear, tmpl::list<>, nonlinear_solver>,
209 : linear_solver, multigrid, schwarz_smoother>>>;
210 :
211 : template <bool Linearized>
212 0 : using dg_operator = elliptic::dg::Actions::DgOperator<
213 : system, Linearized, tmpl::conditional_t<Linearized, vars_tag, fields_tag>,
214 : tmpl::conditional_t<Linearized, fluxes_vars_tag, fluxes_tag>,
215 : tmpl::conditional_t<Linearized, operator_applied_to_vars_tag,
216 : operator_applied_to_fields_tag>>;
217 :
218 0 : using build_matrix = LinearSolver::Actions::BuildMatrix<
219 : typename multigrid::smooth_fields_tag,
220 : typename multigrid::smooth_source_tag, vars_tag,
221 : operator_applied_to_vars_tag,
222 : domain::Tags::Coordinates<volume_dim, Frame::Inertial>,
223 : ::amr::Tags::GridIndex>;
224 :
225 0 : using build_matrix_actions = typename build_matrix::template actions<
226 : typename dg_operator<true>::apply_actions>;
227 :
228 : // For labeling the yaml option for RandomizeVariables
229 0 : struct RandomizeInitialGuess {};
230 :
231 0 : using init_analytic_solution_action =
232 : elliptic::Actions::InitializeOptionalAnalyticSolution<
233 : volume_dim, background_tag,
234 : tmpl::append<typename system::primal_fields,
235 : typename system::primal_fluxes>,
236 : elliptic::analytic_data::AnalyticSolution>;
237 :
238 0 : using initialization_actions = tmpl::list<
239 : elliptic::dg::Actions::InitializeDomain<volume_dim>,
240 : tmpl::conditional_t<is_linear, tmpl::list<>,
241 : typename nonlinear_solver::initialize_element>,
242 : typename linear_solver::initialize_element,
243 : typename multigrid::initialize_element,
244 : typename schwarz_smoother::initialize_element,
245 : Initialization::Actions::InitializeItems<
246 : ::amr::Initialization::Initialize<volume_dim, Metavariables>,
247 : elliptic::amr::Actions::Initialize>,
248 : elliptic::Actions::InitializeFields<system, initial_guess_tag>,
249 : ::Actions::RandomizeVariables<fields_tag, RandomizeInitialGuess>,
250 : init_analytic_solution_action,
251 : elliptic::dg::Actions::initialize_operator<system, background_tag>,
252 : // Actions below may use faces and normals
253 : elliptic::Actions::InitializeFixedSources<system, background_tag>,
254 : tmpl::conditional_t<is_linear,
255 : tmpl::list<::Actions::MutateApply<
256 : elliptic::dg::Actions::
257 : ImposeInhomogeneousBoundaryConditionsOnSource<
258 : system, fixed_sources_tag>>>,
259 : tmpl::list<>>,
260 : tmpl::conditional_t<
261 : is_linear, tmpl::list<>,
262 : ::Initialization::Actions::AddComputeTags<tmpl::list<
263 : // For linearized boundary conditions
264 : elliptic::Tags::BoundaryFieldsCompute<volume_dim, fields_tag>,
265 : elliptic::Tags::BoundaryFluxesCompute<volume_dim, fields_tag,
266 : fluxes_tag>>>>>;
267 :
268 0 : using register_actions = tmpl::flatten<tmpl::list<
269 : ::amr::Actions::RegisterElement,
270 : tmpl::conditional_t<is_linear, tmpl::list<>,
271 : typename nonlinear_solver::register_element>,
272 : typename multigrid::register_element,
273 : typename schwarz_smoother::register_element,
274 : typename build_matrix::register_actions>>;
275 :
276 : template <typename Label>
277 0 : using smooth_actions = typename schwarz_smoother::template solve<
278 : typename dg_operator<true>::apply_actions, Label>;
279 :
280 : // These tags are communicated on subdomain overlaps to initialize the
281 : // subdomain geometry. AMR updates these tags, so we have to communicate them
282 : // after each AMR step.
283 0 : using subdomain_init_tags =
284 : tmpl::list<domain::Tags::Mesh<volume_dim>,
285 : domain::Tags::Element<volume_dim>,
286 : domain::Tags::NeighborMesh<volume_dim>>;
287 :
288 : /// This data needs to be communicated on subdomain overlap regions
289 1 : using communicated_overlap_tags = tmpl::flatten<tmpl::list<
290 : // For linearized sources
291 : fields_tag, fluxes_tag,
292 : // For linearized boundary conditions
293 : domain::make_faces_tags<volume_dim, typename system::primal_fields>,
294 : domain::make_faces_tags<
295 : volume_dim, db::wrap_tags_in<::Tags::NormalDotFlux,
296 : typename system::primal_fields>>>>;
297 :
298 0 : using init_subdomain_action =
299 : elliptic::dg::subdomain_operator::Actions::InitializeSubdomain<
300 : system, background_tag, typename schwarz_smoother::options_group,
301 : false>;
302 :
303 : template <typename StepActions>
304 0 : using linear_solve_actions = typename linear_solver::template solve<
305 : tmpl::list<
306 : // Multigrid preconditioning
307 : typename multigrid::template solve<
308 : typename dg_operator<true>::apply_actions,
309 : // Schwarz smoothing on each multigrid level
310 : smooth_actions<LinearSolver::multigrid::VcycleDownLabel>,
311 : smooth_actions<LinearSolver::multigrid::VcycleUpLabel>,
312 : build_matrix_actions>,
313 : // Support disabling the preconditioner
314 : ::LinearSolver::Actions::make_identity_if_skipped<
315 : multigrid, typename dg_operator<true>::apply_actions>>,
316 : StepActions>;
317 :
318 : template <typename StepActions>
319 0 : using nonlinear_solve_actions = typename nonlinear_solver::template solve<
320 : typename dg_operator<false>::apply_actions,
321 : tmpl::list<
322 : // Transfer data down the multigrid hierachy
323 : LinearSolver::multigrid::Actions::ReceiveFieldsFromFinerGrid<
324 : volume_dim, tmpl::list<fields_tag, fluxes_tag>,
325 : typename multigrid::options_group>,
326 : LinearSolver::multigrid::Actions::SendFieldsToCoarserGrid<
327 : tmpl::list<fields_tag, fluxes_tag>,
328 : typename multigrid::options_group, void>,
329 : // Communicate data on subdomain overlap regions
330 : LinearSolver::Schwarz::Actions::SendOverlapFields<
331 : communicated_overlap_tags,
332 : typename schwarz_smoother::options_group, false,
333 : nonlinear_solver_iteration_id>,
334 : LinearSolver::Schwarz::Actions::ReceiveOverlapFields<
335 : volume_dim, communicated_overlap_tags,
336 : typename schwarz_smoother::options_group, false,
337 : nonlinear_solver_iteration_id>,
338 : // Reset Schwarz subdomain solver
339 : LinearSolver::Schwarz::Actions::ResetSubdomainSolver<
340 : typename schwarz_smoother::subdomain_solver,
341 : typename schwarz_smoother::options_group>,
342 : // Reset explicitly built matrix
343 : typename build_matrix::reset_actions,
344 : // Linear solve for correction
345 : linear_solve_actions<tmpl::list<>>>,
346 : StepActions>;
347 :
348 : template <typename StepActions>
349 0 : using solve_actions = tmpl::flatten<tmpl::list<
350 : // Observe initial state
351 : elliptic::Actions::RunEventsAndTriggers<amr_iteration_id>,
352 : // Stop AMR iterations if complete
353 : elliptic::amr::Actions::StopAmr,
354 : // Run AMR or other phase changes if requested
355 : PhaseControl::Actions::ExecutePhaseChange, StepActions,
356 : // Increment AMR iteration ID
357 : elliptic::amr::Actions::IncrementIterationId,
358 : // Communicate subdomain geometry and initialize subdomain to account for
359 : // domain changes
360 : LinearSolver::Schwarz::Actions::SendOverlapFields<
361 : subdomain_init_tags, typename schwarz_smoother::options_group, false,
362 : amr_iteration_id>,
363 : LinearSolver::Schwarz::Actions::ReceiveOverlapFields<
364 : volume_dim, subdomain_init_tags,
365 : typename schwarz_smoother::options_group, false, amr_iteration_id>,
366 : init_subdomain_action,
367 : // Linear or nonlinear solve
368 : tmpl::conditional_t<is_linear,
369 : // Linear solve
370 : tmpl::list<
371 : // Apply the DG operator to the initial guess
372 : typename elliptic::dg::Actions::DgOperator<
373 : system, true, fields_tag, fluxes_vars_tag,
374 : operator_applied_to_fields_tag, vars_tag,
375 : fluxes_vars_tag>::apply_actions,
376 : // Krylov solve
377 : linear_solve_actions<tmpl::list<>>>,
378 : // Nonlinear solve
379 : nonlinear_solve_actions<tmpl::list<>>>>>;
380 :
381 : template <typename Tag>
382 0 : using overlaps_tag =
383 : LinearSolver::Schwarz::Tags::Overlaps<Tag, volume_dim,
384 : OptionTags::SchwarzSmootherGroup>;
385 :
386 0 : using amr_projectors = tmpl::flatten<tmpl::list<
387 : elliptic::dg::ProjectGeometry<volume_dim>,
388 : tmpl::conditional_t<is_linear, tmpl::list<>,
389 : typename nonlinear_solver::amr_projectors>,
390 : typename linear_solver::amr_projectors,
391 : typename multigrid::amr_projectors,
392 : typename schwarz_smoother::amr_projectors,
393 : typename build_matrix::amr_projectors,
394 : ::amr::projectors::DefaultInitialize<tmpl::append<
395 : tmpl::list<
396 : domain::Tags::InitialExtents<volume_dim>,
397 : domain::Tags::InitialRefinementLevels<volume_dim>,
398 : db::add_tag_prefix<::Tags::deriv, fields_tag,
399 : tmpl::size_t<volume_dim>, Frame::Inertial>>,
400 : // Tags communicated on subdomain overlaps. No need to project
401 : // these during AMR because they will be communicated.
402 : db::wrap_tags_in<
403 : overlaps_tag,
404 : tmpl::append<subdomain_init_tags,
405 : tmpl::conditional_t<is_linear, tmpl::list<>,
406 : communicated_overlap_tags>>>,
407 : // Tags initialized on subdomains. No need to project these during
408 : // AMR because they will get re-initialized after communication.
409 : typename init_subdomain_action::simple_tags>>,
410 : ::amr::projectors::ProjectVariables<volume_dim, fields_tag>,
411 : init_analytic_solution_action,
412 : elliptic::dg::Actions::amr_projectors<system, background_tag>,
413 : typename dg_operator<true>::amr_projectors,
414 : tmpl::conditional_t<is_linear, tmpl::list<>,
415 : typename dg_operator<false>::amr_projectors>,
416 : // Actions below may use faces and normals
417 : elliptic::Actions::InitializeFixedSources<system, background_tag>,
418 : tmpl::conditional_t<
419 : is_linear,
420 : tmpl::list<elliptic::dg::Actions::
421 : ImposeInhomogeneousBoundaryConditionsOnSource<
422 : system, fixed_sources_tag>>,
423 : tmpl::list<>>,
424 : elliptic::amr::Actions::Initialize>>;
425 :
426 0 : using component_list = tmpl::flatten<
427 : tmpl::list<tmpl::conditional_t<is_linear, tmpl::list<>,
428 : typename nonlinear_solver::component_list>,
429 : typename linear_solver::component_list,
430 : typename multigrid::component_list,
431 : typename schwarz_smoother::component_list,
432 : typename build_matrix::template component_list<Metavariables>,
433 : ::amr::Component<Metavariables>>>;
434 : };
435 :
436 : } // namespace elliptic
|
