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/RandomizeVariables.hpp"
35 : #include "ParallelAlgorithms/Amr/Actions/Component.hpp"
36 : #include "ParallelAlgorithms/Amr/Actions/Initialize.hpp"
37 : #include "ParallelAlgorithms/Amr/Projectors/DefaultInitialize.hpp"
38 : #include "ParallelAlgorithms/Amr/Projectors/Variables.hpp"
39 : #include "ParallelAlgorithms/Amr/Tags.hpp"
40 : #include "ParallelAlgorithms/LinearSolver/Actions/BuildMatrix.hpp"
41 : #include "ParallelAlgorithms/LinearSolver/Actions/MakeIdentityIfSkipped.hpp"
42 : #include "ParallelAlgorithms/LinearSolver/Gmres/Gmres.hpp"
43 : #include "ParallelAlgorithms/LinearSolver/Multigrid/Actions/RestrictFields.hpp"
44 : #include "ParallelAlgorithms/LinearSolver/Multigrid/Multigrid.hpp"
45 : #include "ParallelAlgorithms/LinearSolver/Schwarz/Actions/CommunicateOverlapFields.hpp"
46 : #include "ParallelAlgorithms/LinearSolver/Schwarz/Actions/ResetSubdomainSolver.hpp"
47 : #include "ParallelAlgorithms/LinearSolver/Schwarz/Schwarz.hpp"
48 : #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
49 : #include "ParallelAlgorithms/NonlinearSolver/NewtonRaphson/NewtonRaphson.hpp"
50 : #include "ParallelAlgorithms/NonlinearSolver/Tags.hpp"
51 : #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
52 : #include "PointwiseFunctions/InitialDataUtilities/Background.hpp"
53 : #include "PointwiseFunctions/InitialDataUtilities/InitialGuess.hpp"
54 : #include "Utilities/ProtocolHelpers.hpp"
55 : #include "Utilities/TMPL.hpp"
56 :
57 : /// Items related to composing nonlinear elliptic solver executables
58 : namespace elliptic {
59 :
60 : /// Option tags for nonlinear elliptic solver executables
61 : namespace OptionTags {
62 :
63 0 : struct NonlinearSolverGroup {
64 0 : static std::string name() { return "NonlinearSolver"; }
65 0 : static constexpr Options::String help = "The iterative nonlinear solver";
66 : };
67 0 : struct NewtonRaphsonGroup {
68 0 : static std::string name() { return "NewtonRaphson"; }
69 0 : static constexpr Options::String help =
70 : "Options for the Newton-Raphson nonlinear solver";
71 0 : using group = NonlinearSolverGroup;
72 : };
73 0 : struct LinearSolverGroup {
74 0 : static std::string name() { return "LinearSolver"; }
75 0 : static constexpr Options::String help =
76 : "The iterative Krylov-subspace linear solver";
77 : };
78 0 : struct GmresGroup {
79 0 : static std::string name() { return "Gmres"; }
80 0 : static constexpr Options::String help = "Options for the GMRES linear solver";
81 0 : using group = LinearSolverGroup;
82 : };
83 0 : struct SchwarzSmootherGroup {
84 0 : static std::string name() { return "SchwarzSmoother"; }
85 0 : static constexpr Options::String help = "Options for the Schwarz smoother";
86 0 : using group = LinearSolverGroup;
87 : };
88 0 : struct MultigridGroup {
89 0 : static std::string name() { return "Multigrid"; }
90 0 : static constexpr Options::String help = "Options for the multigrid";
91 0 : using group = LinearSolverGroup;
92 : };
93 :
94 : } // namespace OptionTags
95 :
96 : /*!
97 : * \brief A complete nonlinear elliptic solver stack. Use to compose an
98 : * executable.
99 : *
100 : * The elliptic solver stack is described in detail in \cite Vu2021coj.
101 : *
102 : * Uses `Metavariables::volume_dim` and `Metavariables::system`. Also uses
103 : * `Metavariables` to instantiate parallel components.
104 : */
105 : template <typename Metavariables>
106 1 : struct Solver {
107 0 : static constexpr size_t volume_dim = Metavariables::volume_dim;
108 0 : using system = typename Metavariables::system;
109 : static_assert(
110 : tt::assert_conforms_to_v<system, elliptic::protocols::FirstOrderSystem>);
111 0 : static constexpr bool is_linear =
112 : std::is_same_v<elliptic::get_sources_computer<system, true>,
113 : typename system::sources_computer>;
114 :
115 0 : using background_tag =
116 : elliptic::Tags::Background<elliptic::analytic_data::Background>;
117 0 : using initial_guess_tag =
118 : elliptic::Tags::InitialGuess<elliptic::analytic_data::InitialGuess>;
119 :
120 : /// These are the fields we solve for
121 1 : using fields_tag = ::Tags::Variables<typename system::primal_fields>;
122 : /// These are the fluxes corresponding to the fields, i.e. essentially their
123 : /// first derivatives. These are background fields for the linearized sources.
124 1 : using fluxes_tag = ::Tags::Variables<typename system::primal_fluxes>;
125 : /// These are the fixed sources, i.e. the RHS of the equations
126 1 : using fixed_sources_tag = db::add_tag_prefix<::Tags::FixedSource, fields_tag>;
127 0 : using operator_applied_to_fields_tag = tmpl::conditional_t<
128 : is_linear,
129 : // This is the linear operator applied to the fields. We only use it to
130 : // apply the operator to the initial guess, so an optimization would be to
131 : // re-use the `operator_applied_to_vars_tag` below. This optimization
132 : // needs a few minor changes to the parallel linear solver algorithm.
133 : db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, fields_tag>,
134 : db::add_tag_prefix<NonlinearSolver::Tags::OperatorAppliedTo, fields_tag>>;
135 :
136 : // For AMR iterations
137 0 : using amr_iteration_id =
138 : Convergence::Tags::IterationId<::amr::OptionTags::AmrGroup>;
139 :
140 : /// The nonlinear solver algorithm
141 1 : using nonlinear_solver = NonlinearSolver::newton_raphson::NewtonRaphson<
142 : Metavariables, fields_tag, OptionTags::NewtonRaphsonGroup,
143 : fixed_sources_tag, LinearSolver::multigrid::Tags::IsFinestGrid>;
144 0 : using nonlinear_solver_iteration_id =
145 : Convergence::Tags::IterationId<typename nonlinear_solver::options_group>;
146 :
147 : /// The linear solver algorithm. We use GMRES since the operator is not
148 : /// necessarily symmetric. Using CG here would be an optimization for
149 : /// symmetric problems.
150 1 : using linear_solver = LinearSolver::gmres::Gmres<
151 : Metavariables,
152 : tmpl::conditional_t<is_linear, fields_tag,
153 : typename nonlinear_solver::linear_solver_fields_tag>,
154 : OptionTags::GmresGroup, true,
155 : tmpl::conditional_t<is_linear, fixed_sources_tag,
156 : typename nonlinear_solver::linear_solver_source_tag>,
157 : LinearSolver::multigrid::Tags::IsFinestGrid>;
158 0 : using linear_solver_iteration_id =
159 : Convergence::Tags::IterationId<typename linear_solver::options_group>;
160 :
161 : /// Precondition each linear solver iteration with a multigrid V-cycle
162 1 : using multigrid = LinearSolver::multigrid::Multigrid<
163 : Metavariables, typename linear_solver::operand_tag,
164 : OptionTags::MultigridGroup, elliptic::dg::Tags::Massive,
165 : typename linear_solver::preconditioner_source_tag>;
166 :
167 : /// Smooth each multigrid level with a number of Schwarz smoothing steps
168 1 : using subdomain_operator =
169 : elliptic::dg::subdomain_operator::SubdomainOperator<
170 : system, OptionTags::SchwarzSmootherGroup>;
171 0 : using subdomain_preconditioners = tmpl::list<
172 : elliptic::subdomain_preconditioners::Registrars::MinusLaplacian<
173 : volume_dim, OptionTags::SchwarzSmootherGroup>>;
174 0 : using schwarz_smoother = LinearSolver::Schwarz::Schwarz<
175 : typename multigrid::smooth_fields_tag, OptionTags::SchwarzSmootherGroup,
176 : subdomain_operator, subdomain_preconditioners,
177 : typename multigrid::smooth_source_tag,
178 : LinearSolver::multigrid::Tags::MultigridLevel>;
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 nonlinear_solver::linear_solver_fields_tag::tags_list,
202 : typename nonlinear_solver::linear_solver_source_tag::tags_list>>;
203 :
204 : /// Collect all reduction tags for observers
205 1 : using observed_reduction_data_tags =
206 : observers::collect_reduction_data_tags<tmpl::flatten<tmpl::list<
207 : tmpl::conditional_t<is_linear, tmpl::list<>, nonlinear_solver>,
208 : linear_solver, multigrid, schwarz_smoother>>>;
209 :
210 : template <bool Linearized>
211 0 : using dg_operator = elliptic::dg::Actions::DgOperator<
212 : system, Linearized,
213 : tmpl::conditional_t<Linearized, linear_solver_iteration_id,
214 : nonlinear_solver_iteration_id>,
215 : tmpl::conditional_t<Linearized, vars_tag, fields_tag>,
216 : tmpl::conditional_t<Linearized, fluxes_vars_tag, fluxes_tag>,
217 : tmpl::conditional_t<Linearized, operator_applied_to_vars_tag,
218 : operator_applied_to_fields_tag>>;
219 :
220 0 : using build_matrix = LinearSolver::Actions::BuildMatrix<
221 : linear_solver_iteration_id, vars_tag, operator_applied_to_vars_tag,
222 : domain::Tags::Coordinates<volume_dim, Frame::Inertial>,
223 : LinearSolver::multigrid::Tags::IsFinestGrid>;
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>,
247 : elliptic::amr::Actions::Initialize>,
248 : elliptic::Actions::InitializeFields<system, initial_guess_tag>,
249 : ::Actions::RandomizeVariables<fields_tag, RandomizeInitialGuess>,
250 : elliptic::Actions::InitializeFixedSources<system, background_tag>,
251 : init_analytic_solution_action,
252 : elliptic::dg::Actions::initialize_operator<system, background_tag>,
253 : tmpl::conditional_t<
254 : is_linear, tmpl::list<>,
255 : ::Initialization::Actions::AddComputeTags<tmpl::list<
256 : // For linearized boundary conditions
257 : elliptic::Tags::BoundaryFieldsCompute<volume_dim, fields_tag>,
258 : elliptic::Tags::BoundaryFluxesCompute<volume_dim, fields_tag,
259 : fluxes_tag>>>>>;
260 :
261 0 : using register_actions = tmpl::flatten<tmpl::list<
262 : tmpl::conditional_t<is_linear, typename build_matrix::register_actions,
263 : typename nonlinear_solver::register_element>,
264 : typename multigrid::register_element,
265 : typename schwarz_smoother::register_element>>;
266 :
267 : template <typename Label>
268 0 : using smooth_actions = typename schwarz_smoother::template solve<
269 : typename dg_operator<true>::apply_actions, Label>;
270 :
271 : // These tags are communicated on subdomain overlaps to initialize the
272 : // subdomain geometry. AMR updates these tags, so we have to communicate them
273 : // after each AMR step.
274 0 : using subdomain_init_tags =
275 : tmpl::list<domain::Tags::Mesh<volume_dim>,
276 : domain::Tags::Element<volume_dim>,
277 : domain::Tags::NeighborMesh<volume_dim>>;
278 :
279 : /// This data needs to be communicated on subdomain overlap regions
280 1 : using communicated_overlap_tags = tmpl::flatten<tmpl::list<
281 : // For linearized sources
282 : fields_tag, fluxes_tag,
283 : // For linearized boundary conditions
284 : domain::make_faces_tags<volume_dim, typename system::primal_fields>,
285 : domain::make_faces_tags<
286 : volume_dim, db::wrap_tags_in<::Tags::NormalDotFlux,
287 : typename system::primal_fields>>>>;
288 :
289 0 : using init_subdomain_action =
290 : elliptic::dg::subdomain_operator::Actions::InitializeSubdomain<
291 : system, background_tag, typename schwarz_smoother::options_group,
292 : false>;
293 :
294 : template <typename StepActions>
295 0 : using linear_solve_actions = typename linear_solver::template solve<
296 : tmpl::list<
297 : // Multigrid preconditioning
298 : typename multigrid::template solve<
299 : typename dg_operator<true>::apply_actions,
300 : // Schwarz smoothing on each multigrid level
301 : smooth_actions<LinearSolver::multigrid::VcycleDownLabel>,
302 : smooth_actions<LinearSolver::multigrid::VcycleUpLabel>>,
303 : // Support disabling the preconditioner
304 : ::LinearSolver::Actions::make_identity_if_skipped<
305 : multigrid, typename dg_operator<true>::apply_actions>>,
306 : StepActions>;
307 :
308 : template <typename StepActions>
309 0 : using nonlinear_solve_actions = typename nonlinear_solver::template solve<
310 : typename dg_operator<false>::apply_actions,
311 : tmpl::list<
312 : // Transfer data down the multigrid hierachy
313 : LinearSolver::multigrid::Actions::ReceiveFieldsFromFinerGrid<
314 : volume_dim, tmpl::list<fields_tag, fluxes_tag>,
315 : typename multigrid::options_group>,
316 : LinearSolver::multigrid::Actions::SendFieldsToCoarserGrid<
317 : tmpl::list<fields_tag, fluxes_tag>,
318 : typename multigrid::options_group, void>,
319 : // Communicate data on subdomain overlap regions
320 : LinearSolver::Schwarz::Actions::SendOverlapFields<
321 : communicated_overlap_tags,
322 : typename schwarz_smoother::options_group, false,
323 : nonlinear_solver_iteration_id>,
324 : LinearSolver::Schwarz::Actions::ReceiveOverlapFields<
325 : volume_dim, communicated_overlap_tags,
326 : typename schwarz_smoother::options_group, false,
327 : nonlinear_solver_iteration_id>,
328 : // Reset Schwarz subdomain solver
329 : LinearSolver::Schwarz::Actions::ResetSubdomainSolver<
330 : typename schwarz_smoother::options_group>,
331 : // Linear solve for correction
332 : linear_solve_actions<tmpl::list<>>>,
333 : StepActions>;
334 :
335 : template <typename StepActions>
336 0 : using solve_actions = tmpl::flatten<tmpl::list<
337 : // Observe initial state
338 : elliptic::Actions::RunEventsAndTriggers<amr_iteration_id>,
339 : // Stop AMR iterations if complete
340 : elliptic::amr::Actions::StopAmr,
341 : // Run AMR or other phase changes if requested
342 : PhaseControl::Actions::ExecutePhaseChange, StepActions,
343 : // Increment AMR iteration ID
344 : elliptic::amr::Actions::IncrementIterationId,
345 : // Communicate subdomain geometry and initialize subdomain to account for
346 : // domain changes
347 : LinearSolver::Schwarz::Actions::SendOverlapFields<
348 : subdomain_init_tags, typename schwarz_smoother::options_group, false,
349 : amr_iteration_id>,
350 : LinearSolver::Schwarz::Actions::ReceiveOverlapFields<
351 : volume_dim, subdomain_init_tags,
352 : typename schwarz_smoother::options_group, false, amr_iteration_id>,
353 : init_subdomain_action,
354 : // Linear or nonlinear solve
355 : tmpl::conditional_t<
356 : is_linear,
357 : // Linear solve
358 : tmpl::list<
359 : // Apply the DG operator to the initial guess
360 : typename elliptic::dg::Actions::DgOperator<
361 : system, true, linear_solver_iteration_id, fields_tag,
362 : fluxes_vars_tag, operator_applied_to_fields_tag, vars_tag,
363 : fluxes_vars_tag>::apply_actions,
364 : // Modify fixed sources with boundary conditions
365 : elliptic::dg::Actions::
366 : ImposeInhomogeneousBoundaryConditionsOnSource<
367 : system, fixed_sources_tag>,
368 : // Krylov solve
369 : linear_solve_actions<tmpl::list<>>>,
370 : // Nonlinear solve
371 : nonlinear_solve_actions<tmpl::list<>>>>>;
372 :
373 : template <typename Tag>
374 0 : using overlaps_tag =
375 : LinearSolver::Schwarz::Tags::Overlaps<Tag, volume_dim,
376 : OptionTags::SchwarzSmootherGroup>;
377 :
378 0 : using amr_projectors = tmpl::flatten<tmpl::list<
379 : elliptic::dg::ProjectGeometry<volume_dim>,
380 : tmpl::conditional_t<is_linear, tmpl::list<>,
381 : typename nonlinear_solver::amr_projectors>,
382 : typename linear_solver::amr_projectors,
383 : typename multigrid::amr_projectors,
384 : typename schwarz_smoother::amr_projectors,
385 : ::amr::projectors::DefaultInitialize<tmpl::append<
386 : tmpl::list<domain::Tags::InitialExtents<volume_dim>,
387 : domain::Tags::InitialRefinementLevels<volume_dim>>,
388 : // Tags communicated on subdomain overlaps. No need to project
389 : // these during AMR because they will be communicated.
390 : db::wrap_tags_in<
391 : overlaps_tag,
392 : tmpl::append<subdomain_init_tags,
393 : tmpl::conditional_t<is_linear, tmpl::list<>,
394 : communicated_overlap_tags>>>,
395 : // Tags initialized on subdomains. No need to project these during
396 : // AMR because they will get re-initialized after communication.
397 : typename init_subdomain_action::simple_tags>>,
398 : ::amr::projectors::ProjectVariables<volume_dim, fields_tag>,
399 : elliptic::Actions::InitializeFixedSources<system, background_tag>,
400 : init_analytic_solution_action,
401 : elliptic::dg::Actions::amr_projectors<system, background_tag>,
402 : typename dg_operator<true>::amr_projectors,
403 : tmpl::conditional_t<is_linear, typename build_matrix::amr_projectors,
404 : typename dg_operator<false>::amr_projectors>,
405 : elliptic::amr::Actions::Initialize>>;
406 :
407 0 : using component_list = tmpl::flatten<
408 : tmpl::list<tmpl::conditional_t<is_linear, tmpl::list<>,
409 : typename nonlinear_solver::component_list>,
410 : typename linear_solver::component_list,
411 : typename multigrid::component_list,
412 : typename schwarz_smoother::component_list,
413 : ::amr::Component<Metavariables>>>;
414 : };
415 :
416 : } // namespace elliptic
|
