Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <boost/functional/hash.hpp>
7 : #include <cstddef>
8 : #include <map>
9 : #include <memory>
10 : #include <optional>
11 : #include <tuple>
12 : #include <unordered_map>
13 : #include <utility>
14 : #include <vector>
15 :
16 : #include "DataStructures/DataBox/DataBox.hpp"
17 : #include "Domain/Creators/Tags/ExternalBoundaryConditions.hpp"
18 : #include "Domain/Structure/Direction.hpp"
19 : #include "Domain/Structure/Element.hpp"
20 : #include "Domain/Tags.hpp"
21 : #include "Elliptic/BoundaryConditions/BoundaryCondition.hpp"
22 : #include "Elliptic/BoundaryConditions/BoundaryConditionType.hpp"
23 : #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp"
24 : #include "Elliptic/Systems/Poisson/BoundaryConditions/Robin.hpp"
25 : #include "Elliptic/Systems/Poisson/FirstOrderSystem.hpp"
26 : #include "Elliptic/Systems/Poisson/Tags.hpp"
27 : #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
28 : #include "NumericalAlgorithms/LinearSolver/ExplicitInverse.hpp"
29 : #include "NumericalAlgorithms/LinearSolver/Gmres.hpp"
30 : #include "NumericalAlgorithms/LinearSolver/LinearSolver.hpp"
31 : #include "Options/Auto.hpp"
32 : #include "Options/String.hpp"
33 : #include "ParallelAlgorithms/LinearSolver/Schwarz/ElementCenteredSubdomainData.hpp"
34 : #include "ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp"
35 : #include "Utilities/ErrorHandling/Assert.hpp"
36 : #include "Utilities/MakeWithValue.hpp"
37 : #include "Utilities/Serialization/CharmPupable.hpp"
38 : #include "Utilities/TMPL.hpp"
39 :
40 : /// Linear solvers that approximately invert the
41 : /// `elliptic::dg::subdomain_operator::SubdomainOperator` to make the Schwarz
42 : /// subdomain solver converge faster.
43 : /// \see LinearSolver::Schwarz::Schwarz
44 1 : namespace elliptic::subdomain_preconditioners {
45 :
46 : /// \cond
47 : template <size_t Dim, typename OptionsGroup, typename Solver,
48 : typename LinearSolverRegistrars>
49 : struct MinusLaplacian;
50 : /// \endcond
51 :
52 0 : namespace Registrars {
53 : template <size_t Dim, typename OptionsGroup,
54 : typename Solver = LinearSolver::Serial::LinearSolver<tmpl::list<
55 : ::LinearSolver::Serial::Registrars::Gmres<
56 : ::LinearSolver::Schwarz::ElementCenteredSubdomainData<
57 : Dim, tmpl::list<Poisson::Tags::Field>>>,
58 : ::LinearSolver::Serial::Registrars::ExplicitInverse>>>
59 0 : struct MinusLaplacian {
60 : template <typename LinearSolverRegistrars>
61 0 : using f = subdomain_preconditioners::MinusLaplacian<Dim, OptionsGroup, Solver,
62 : LinearSolverRegistrars>;
63 : };
64 : } // namespace Registrars
65 :
66 : /*!
67 : * \brief Approximate the subdomain operator with a flat-space Laplacian
68 : * for every tensor component separately.
69 : *
70 : * This linear solver applies the `Solver` to every tensor component in
71 : * turn, approximating the subdomain operator with a flat-space Laplacian.
72 : * This can be a lot cheaper than solving the
73 : * full subdomain operator and can provide effective preconditioning for an
74 : * iterative subdomain solver. The approximation is better the closer the
75 : * original PDEs are to a set of decoupled flat-space Poisson equations.
76 : *
77 : * \par Boundary conditions
78 : * At external boundaries we impose homogeneous Dirichlet or Neumann boundary
79 : * conditions on the Poisson sub-problems. For tensor components and element
80 : * faces where the original boundary conditions are of Dirichlet-type we choose
81 : * homogeneous Dirichlet boundary conditions, and for those where the original
82 : * boundary conditions are of Neumann-type we choose homogeneous Neumann
83 : * boundary conditions. This means we may end up with more than one distinct
84 : * Poisson operator on subdomains with external boundaries, one per unique
85 : * combination of element face and boundary-condition type among the tensor
86 : * components.
87 : *
88 : * \tparam Dim Spatial dimension
89 : * \tparam OptionsGroup The options group identifying the
90 : * `LinearSolver::Schwarz::Schwarz` solver that defines the subdomain geometry.
91 : * \tparam Solver Any class that provides a `solve` and a `reset` function,
92 : * but typically a `LinearSolver::Serial::LinearSolver`. The solver will be
93 : * factory-created from input-file options.
94 : */
95 : template <size_t Dim, typename OptionsGroup,
96 : typename Solver = LinearSolver::Serial::LinearSolver<tmpl::list<
97 : ::LinearSolver::Serial::Registrars::Gmres<
98 : ::LinearSolver::Schwarz::ElementCenteredSubdomainData<
99 : Dim, tmpl::list<Poisson::Tags::Field>>>,
100 : ::LinearSolver::Serial::Registrars::ExplicitInverse>>,
101 : typename LinearSolverRegistrars =
102 : tmpl::list<Registrars::MinusLaplacian<Dim, OptionsGroup, Solver>>>
103 1 : class MinusLaplacian
104 : : public LinearSolver::Serial::LinearSolver<LinearSolverRegistrars> {
105 : private:
106 0 : using Base = LinearSolver::Serial::LinearSolver<LinearSolverRegistrars>;
107 0 : using StoredSolverType = tmpl::conditional_t<std::is_abstract_v<Solver>,
108 : std::unique_ptr<Solver>, Solver>;
109 : // Identifies an external block boundary, for boundary conditions
110 0 : using BoundaryId = std::pair<size_t, Direction<Dim>>;
111 :
112 : public:
113 0 : static constexpr size_t volume_dim = Dim;
114 0 : using options_group = OptionsGroup;
115 0 : using poisson_system =
116 : Poisson::FirstOrderSystem<Dim, Poisson::Geometry::FlatCartesian>;
117 0 : using BoundaryConditionsBase =
118 : typename poisson_system::boundary_conditions_base;
119 0 : using SubdomainOperator = elliptic::dg::subdomain_operator::SubdomainOperator<
120 : poisson_system, OptionsGroup,
121 : tmpl::list<Poisson::BoundaryConditions::Robin<Dim>>>;
122 0 : using SubdomainData = ::LinearSolver::Schwarz::ElementCenteredSubdomainData<
123 : Dim, tmpl::list<Poisson::Tags::Field>>;
124 : // Associates Dirichlet or Neumann conditions to every external block
125 : // boundary. For every configuration of this type we'll need a separate
126 : // solver, which may cache the Poisson operator matrix that imposes the
127 : // corresponding boundary conditions. This type is used as key in other maps,
128 : // so it must be hashable. `std::unordered_map` isn't supported by
129 : // `boost::hash`, but `std::map` is.
130 0 : using BoundaryConditionsSignature =
131 : std::map<BoundaryId, elliptic::BoundaryConditionType>;
132 0 : using solver_type = Solver;
133 :
134 0 : struct SolverOptionTag {
135 0 : static std::string name() { return "Solver"; }
136 0 : using type = StoredSolverType;
137 0 : static constexpr Options::String help =
138 : "The linear solver used to invert the Laplace operator. The solver is "
139 : "shared between tensor components with the same type of boundary "
140 : "conditions (Dirichlet-type or Neumann-type).";
141 : };
142 :
143 0 : struct BoundaryConditions {
144 0 : using type = Options::Auto<elliptic::BoundaryConditionType>;
145 0 : static constexpr Options::String help =
146 : "The boundary conditions imposed by the Laplace operator. Specify "
147 : "'Auto' to choose between homogeneous Dirichlet or Neumann boundary "
148 : "conditions automatically, based on the configuration of the the full "
149 : "operator.";
150 : };
151 :
152 0 : using options = tmpl::list<SolverOptionTag, BoundaryConditions>;
153 0 : static constexpr Options::String help =
154 : "Approximate the linear operator with a Laplace operator "
155 : "for every tensor component separately.";
156 :
157 0 : MinusLaplacian() = default;
158 0 : MinusLaplacian(MinusLaplacian&& /*rhs*/) = default;
159 0 : MinusLaplacian& operator=(MinusLaplacian&& /*rhs*/) = default;
160 0 : ~MinusLaplacian() = default;
161 0 : MinusLaplacian(const MinusLaplacian& rhs)
162 : : Base(rhs),
163 : solver_(rhs.clone_solver()),
164 : boundary_condition_type_(rhs.boundary_condition_type_) {}
165 0 : MinusLaplacian& operator=(const MinusLaplacian& rhs) {
166 : Base::operator=(rhs);
167 : solver_ = rhs.clone_solver();
168 : boundary_condition_type_ = rhs.boundary_condition_type_;
169 : return *this;
170 : }
171 :
172 : /// \cond
173 : explicit MinusLaplacian(CkMigrateMessage* m) : Base(m) {}
174 : using PUP::able::register_constructor;
175 : WRAPPED_PUPable_decl_template(MinusLaplacian); // NOLINT
176 : /// \endcond
177 :
178 0 : MinusLaplacian(
179 : StoredSolverType solver,
180 : std::optional<elliptic::BoundaryConditionType> boundary_condition_type)
181 : : solver_(std::move(solver)),
182 : boundary_condition_type_(boundary_condition_type) {}
183 :
184 0 : const Solver& solver() const {
185 : if constexpr (std::is_abstract_v<Solver>) {
186 : return *solver_;
187 : } else {
188 : return solver_;
189 : }
190 : }
191 :
192 : /// The cached solvers for each unique boundary-condition signature. These are
193 : /// only exposed to allow verifying their consistency.
194 1 : const auto& cached_solvers() const { return solvers_; }
195 :
196 : /// Solve the equation \f$Ax=b\f$ by approximating \f$A\f$ with a Laplace
197 : /// operator for every tensor component in \f$x\f$.
198 : template <typename LinearOperator, typename VarsType, typename SourceType,
199 : typename... OperatorArgs>
200 1 : Convergence::HasConverged solve(
201 : gsl::not_null<VarsType*> solution, LinearOperator&& linear_operator,
202 : const SourceType& source,
203 : const std::tuple<OperatorArgs...>& operator_args) const;
204 :
205 0 : void reset() override {
206 : mutable_solver().reset();
207 : // These buffers depend on the operator being solved, so they are reset when
208 : // the operator changes. The other buffers only hold workspace memory so
209 : // they don't need to be reset.
210 : bc_signatures_ = std::nullopt;
211 : solvers_.clear();
212 : boundary_conditions_.clear();
213 : }
214 :
215 : // NOLINTNEXTLINE(google-runtime-references)
216 0 : void pup(PUP::er& p) override {
217 : Base::pup(p);
218 : // Serialize object state
219 : p | solver_;
220 : p | boundary_condition_type_;
221 : // Serialize caches that are possibly expensive to re-create
222 : p | bc_signatures_;
223 : p | solvers_;
224 : // Other properties are memory buffers that don't need to be serialized.
225 : }
226 :
227 0 : std::unique_ptr<Base> get_clone() const override {
228 : return std::make_unique<MinusLaplacian>(*this);
229 : }
230 :
231 : private:
232 : // For each tensor component, get the set of boundary conditions that should
233 : // be imposed by the Laplacian operator on that component, based on the domain
234 : // configuration. These are cached for successive solves of the same operator.
235 : // An empty list means the subdomain has no external boundaries.
236 : template <typename DbTagsList>
237 : const std::vector<BoundaryConditionsSignature>&
238 0 : get_cached_boundary_conditions_signatures(
239 : const db::DataBox<DbTagsList>& box) const {
240 : if (bc_signatures_.has_value()) {
241 : // Use cached value
242 : return *bc_signatures_;
243 : }
244 : bc_signatures_ = std::vector<BoundaryConditionsSignature>{};
245 : const auto& all_boundary_conditions =
246 : db::get<domain::Tags::ExternalBoundaryConditions<Dim>>(box);
247 : const auto collect_bc_signatures = [&bc_signatures = *bc_signatures_,
248 : &override_boundary_condition_type =
249 : boundary_condition_type_,
250 : &all_boundary_conditions](
251 : const BoundaryId& boundary_id) {
252 : if (not bc_signatures.empty() and
253 : bc_signatures[0].find(boundary_id) != bc_signatures[0].end()) {
254 : // We have already handled this block boundary in an earlier
255 : // iteration of the loop. This can happend when the domain is
256 : // h-refined, or when multiple elements in a subdomain face the
257 : // block boundary.
258 : return;
259 : }
260 : // Get the type of boundary condition (Dirichlet or Neumann) for
261 : // each tensor component from the domain configuration
262 : const auto& [block_id, direction] = boundary_id;
263 : const auto original_boundary_condition = dynamic_cast<
264 : const elliptic::BoundaryConditions::BoundaryCondition<Dim>*>(
265 : all_boundary_conditions.at(block_id).at(direction).get());
266 : ASSERT(original_boundary_condition != nullptr,
267 : "The boundary condition in block "
268 : << block_id << ", direction " << direction
269 : << " is not of the expected type "
270 : "'elliptic::BoundaryConditions::BoundaryCondition<"
271 : << Dim << ">'.");
272 : const std::vector<elliptic::BoundaryConditionType> bc_types =
273 : original_boundary_condition->boundary_condition_types();
274 : const size_t num_components = bc_types.size();
275 : if (bc_signatures.empty()) {
276 : bc_signatures.resize(num_components);
277 : } else {
278 : ASSERT(bc_signatures.size() == num_components,
279 : "Unexpected number of boundary-condition types. The "
280 : "boundary condition in block "
281 : << block_id << ", direction " << direction << " returned "
282 : << num_components << " items, but another returned "
283 : << bc_signatures.size()
284 : << " (one for each tensor component).");
285 : }
286 : for (size_t component = 0; component < num_components; ++component) {
287 : bc_signatures[component].emplace(
288 : boundary_id,
289 : override_boundary_condition_type.value_or(bc_types.at(component)));
290 : }
291 : };
292 : // Collect boundary-condition signatures from all elements in the subdomain
293 : const auto& element = db::get<domain::Tags::Element<Dim>>(box);
294 : for (const auto& direction : element.external_boundaries()) {
295 : collect_bc_signatures({element.id().block_id(), direction});
296 : }
297 : const auto& overlap_elements =
298 : db::get<LinearSolver::Schwarz::Tags::Overlaps<
299 : domain::Tags::Element<Dim>, Dim, OptionsGroup>>(box);
300 : for (const auto& [overlap_id, overlap_element] : overlap_elements) {
301 : (void)overlap_id;
302 : for (const auto& direction : overlap_element.external_boundaries()) {
303 : collect_bc_signatures({overlap_element.id().block_id(), direction});
304 : }
305 : }
306 : return *bc_signatures_;
307 : }
308 :
309 : // For a tensor component with the given boundary-condition configuration, get
310 : // the solver and the set of boundary conditions that will be passed to the
311 : // Poisson subdomain operator. The solver is cached for successive solves of
312 : // the same operator. Caching is very important for performance, since the
313 : // solver may build up an explicit matrix representation that is expensive to
314 : // re-create.
315 : std::pair<const Solver&,
316 : const std::unordered_map<BoundaryId, const BoundaryConditionsBase&,
317 : boost::hash<BoundaryId>>&>
318 0 : get_cached_solver_and_boundary_conditions(
319 : const std::vector<BoundaryConditionsSignature>& bc_signatures,
320 : const size_t component) const {
321 : if (bc_signatures.empty()) {
322 : // The subdomain has no external boundaries. We use the original solver.
323 : return {solver(), boundary_conditions_[{}]};
324 : }
325 : // Get the cached solver corresponding to this component's
326 : // boundary-condition configuration
327 : const auto& bc_signature = bc_signatures.at(component);
328 : if (solvers_.find(bc_signature) == solvers_.end()) {
329 : solvers_.emplace(bc_signature, clone_solver());
330 : }
331 : const Solver& cached_solver = [this, &bc_signature]() -> const Solver& {
332 : if constexpr (std::is_abstract_v<Solver>) {
333 : return *solvers_.at(bc_signature);
334 : } else {
335 : return solvers_.at(bc_signature);
336 : }
337 : }();
338 : // Get the cached set of boundary conditions for this component
339 : if (boundary_conditions_.find(bc_signature) == boundary_conditions_.end()) {
340 : auto& bc = boundary_conditions_[bc_signature];
341 : for (const auto& [boundary_id, bc_type] : bc_signature) {
342 : bc.emplace(boundary_id,
343 : bc_type == elliptic::BoundaryConditionType::Dirichlet
344 : ? dirichlet_bc
345 : : neumann_bc);
346 : }
347 : }
348 : return {cached_solver, boundary_conditions_[bc_signature]};
349 : }
350 :
351 0 : Solver& mutable_solver() {
352 : if constexpr (std::is_abstract_v<Solver>) {
353 : return *solver_;
354 : } else {
355 : return solver_;
356 : }
357 : }
358 :
359 0 : StoredSolverType clone_solver() const {
360 : if constexpr (std::is_abstract_v<Solver>) {
361 : return solver_->get_clone();
362 : } else {
363 : return solver_;
364 : }
365 : }
366 :
367 : // The option-constructed solver for the separate Poisson problems. It is
368 : // cloned for each unique boundary-condition configuration of the tensor
369 : // components.
370 0 : StoredSolverType solver_{};
371 0 : std::optional<elliptic::BoundaryConditionType> boundary_condition_type_;
372 :
373 : // These are caches to keep track of the unique boundary-condition
374 : // configurations. Each configuration has its own solver, because the solver
375 : // may cache the operator to speed up repeated solves.
376 : // - The boundary-condition configuration for each tensor component, or an
377 : // empty list if the subdomain has no external boundaries, or `std::nullopt`
378 : // to signal a clean cache.
379 : // NOLINTNEXTLINE(spectre-mutable)
380 : mutable std::optional<std::vector<BoundaryConditionsSignature>>
381 0 : bc_signatures_{};
382 : // - A clone of the `solver_` for each unique boundary-condition configuration
383 : // NOLINTNEXTLINE(spectre-mutable)
384 : mutable std::unordered_map<BoundaryConditionsSignature, StoredSolverType,
385 : boost::hash<BoundaryConditionsSignature>>
386 0 : solvers_{};
387 : // NOLINTNEXTLINE(spectre-mutable)
388 : mutable std::unordered_map<
389 : BoundaryConditionsSignature,
390 : std::unordered_map<BoundaryId, const BoundaryConditionsBase&,
391 : boost::hash<BoundaryId>>,
392 : boost::hash<BoundaryConditionsSignature>>
393 0 : boundary_conditions_{};
394 :
395 : // These are memory buffers that are kept around for repeated solves. Possible
396 : // optimization: Free this memory at appropriate times, e.g. when the element
397 : // has completed a full subdomain solve and goes to the background. In some
398 : // cases the `subdomain_operator_` is never even used again in subsequent
399 : // subdomain solves because it is cached as a matrix (see
400 : // LinearSolver::Serial::ExplicitInverse), so we don't need the memory anymore
401 : // at all.
402 : // NOLINTNEXTLINE(spectre-mutable)
403 0 : mutable SubdomainOperator subdomain_operator_{};
404 : // NOLINTNEXTLINE(spectre-mutable)
405 0 : mutable SubdomainData source_{};
406 : // NOLINTNEXTLINE(spectre-mutable)
407 0 : mutable SubdomainData initial_guess_in_solution_out_{};
408 :
409 : // These boundary condition instances can be re-used for all tensor components
410 : // Note that the _linearized_ boundary conditions are applied by the subdomain
411 : // operator, so the constant (third argument) is irrelevant.
412 0 : const Poisson::BoundaryConditions::Robin<Dim> dirichlet_bc{1., 0., 0.};
413 0 : const Poisson::BoundaryConditions::Robin<Dim> neumann_bc{0., 1., 0.};
414 : };
415 :
416 : namespace detail {
417 : template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
418 : void assign_component(
419 : const gsl::not_null<::LinearSolver::Schwarz::ElementCenteredSubdomainData<
420 : Dim, LhsTagsList>*>
421 : lhs,
422 : const ::LinearSolver::Schwarz::ElementCenteredSubdomainData<
423 : Dim, RhsTagsList>& rhs,
424 : const size_t lhs_component, const size_t rhs_component) {
425 : // Possible optimization: Once we have non-owning Variables we can use a view
426 : // into the rhs here instead of copying.
427 : const size_t num_points_element = rhs.element_data.number_of_grid_points();
428 : for (size_t i = 0; i < num_points_element; ++i) {
429 : lhs->element_data.data()[lhs_component * num_points_element + i] =
430 : rhs.element_data.data()[rhs_component * num_points_element + i];
431 : }
432 : for (const auto& [overlap_id, rhs_data] : rhs.overlap_data) {
433 : const size_t num_points_overlap = rhs_data.number_of_grid_points();
434 : // The random-access operation is relatively slow because it computes a
435 : // hash, so it's important for performance to avoid repeating it in every
436 : // iteration of the loop below.
437 : auto& lhs_vars = lhs->overlap_data[overlap_id];
438 : for (size_t i = 0; i < num_points_overlap; ++i) {
439 : lhs_vars.data()[lhs_component * num_points_overlap + i] =
440 : rhs_data.data()[rhs_component * num_points_overlap + i];
441 : }
442 : }
443 : }
444 : } // namespace detail
445 :
446 : template <size_t Dim, typename OptionsGroup, typename Solver,
447 : typename LinearSolverRegistrars>
448 : template <typename LinearOperator, typename VarsType, typename SourceType,
449 : typename... OperatorArgs>
450 : Convergence::HasConverged
451 0 : MinusLaplacian<Dim, OptionsGroup, Solver, LinearSolverRegistrars>::solve(
452 : const gsl::not_null<VarsType*> initial_guess_in_solution_out,
453 : LinearOperator&& /*linear_operator*/, const SourceType& source,
454 : const std::tuple<OperatorArgs...>& operator_args) const {
455 : const auto& box = get<0>(operator_args);
456 : // Solve each component of the source variables in turn, assuming the operator
457 : // is a Laplacian. For each component we select either homogeneous Dirichlet
458 : // or Neumann boundary conditions, based on the type of boundary conditions
459 : // imposed by the full operator.
460 : static constexpr size_t num_components =
461 : VarsType::ElementData::number_of_independent_components;
462 : source_.destructive_resize(source);
463 : initial_guess_in_solution_out_.destructive_resize(source);
464 : const std::vector<BoundaryConditionsSignature>& bc_signatures =
465 : get_cached_boundary_conditions_signatures(box);
466 : ASSERT(bc_signatures.empty() or bc_signatures.size() == num_components,
467 : "Expected " << num_components
468 : << " boundary-condition signatures (one per tensor "
469 : "component), but got "
470 : << bc_signatures.size() << ".");
471 : // Possible optimization: elide when num_components is 1
472 : for (size_t component = 0; component < num_components; ++component) {
473 : detail::assign_component(make_not_null(&source_), source, 0, component);
474 : detail::assign_component(make_not_null(&initial_guess_in_solution_out_),
475 : *initial_guess_in_solution_out, 0, component);
476 : const auto& [solver, boundary_conditions] =
477 : get_cached_solver_and_boundary_conditions(bc_signatures, component);
478 : solver.solve(make_not_null(&initial_guess_in_solution_out_),
479 : subdomain_operator_, source_,
480 : std::forward_as_tuple(box, boundary_conditions));
481 : detail::assign_component(initial_guess_in_solution_out,
482 : initial_guess_in_solution_out_, component, 0);
483 : }
484 : return {0, 0};
485 : }
486 :
487 : /// \cond
488 : template <size_t Dim, typename OptionsGroup, typename Solver,
489 : typename LinearSolverRegistrars>
490 : // NOLINTNEXTLINE
491 : PUP::able::PUP_ID MinusLaplacian<Dim, OptionsGroup, Solver,
492 : LinearSolverRegistrars>::my_PUP_ID = 0;
493 : /// \endcond
494 :
495 : } // namespace elliptic::subdomain_preconditioners
|
