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<DataVector>>>>,
58 : ::LinearSolver::Serial::Registrars::ExplicitInverse<double>>>>
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 : * \par Complex variables
89 : * This preconditioner can also be applied to complex variables, in which case
90 : * the Laplacian just applies to the real and imaginary part separately.
91 : *
92 : * \tparam Dim Spatial dimension
93 : * \tparam OptionsGroup The options group identifying the
94 : * `LinearSolver::Schwarz::Schwarz` solver that defines the subdomain geometry.
95 : * \tparam Solver Any class that provides a `solve` and a `reset` function,
96 : * but typically a `LinearSolver::Serial::LinearSolver`. The solver will be
97 : * factory-created from input-file options.
98 : */
99 : template <size_t Dim, typename OptionsGroup,
100 : typename Solver = LinearSolver::Serial::LinearSolver<tmpl::list<
101 : ::LinearSolver::Serial::Registrars::Gmres<
102 : ::LinearSolver::Schwarz::ElementCenteredSubdomainData<
103 : Dim, tmpl::list<Poisson::Tags::Field<DataVector>>>>,
104 : ::LinearSolver::Serial::Registrars::ExplicitInverse<double>>>,
105 : typename LinearSolverRegistrars =
106 : tmpl::list<Registrars::MinusLaplacian<Dim, OptionsGroup, Solver>>>
107 1 : class MinusLaplacian
108 : : public LinearSolver::Serial::LinearSolver<LinearSolverRegistrars> {
109 : private:
110 0 : using Base = LinearSolver::Serial::LinearSolver<LinearSolverRegistrars>;
111 0 : using StoredSolverType = tmpl::conditional_t<std::is_abstract_v<Solver>,
112 : std::unique_ptr<Solver>, Solver>;
113 : // Identifies an external block boundary, for boundary conditions
114 0 : using BoundaryId = std::pair<size_t, Direction<Dim>>;
115 :
116 : public:
117 0 : static constexpr size_t volume_dim = Dim;
118 0 : using options_group = OptionsGroup;
119 0 : using poisson_system =
120 : Poisson::FirstOrderSystem<Dim, Poisson::Geometry::FlatCartesian>;
121 0 : using BoundaryConditionsBase =
122 : typename poisson_system::boundary_conditions_base;
123 0 : using SubdomainOperator = elliptic::dg::subdomain_operator::SubdomainOperator<
124 : poisson_system, OptionsGroup,
125 : tmpl::list<Poisson::BoundaryConditions::Robin<Dim>>>;
126 0 : using SubdomainData = ::LinearSolver::Schwarz::ElementCenteredSubdomainData<
127 : Dim, tmpl::list<Poisson::Tags::Field<DataVector>>>;
128 : // Associates Dirichlet or Neumann conditions to every external block
129 : // boundary. For every configuration of this type we'll need a separate
130 : // solver, which may cache the Poisson operator matrix that imposes the
131 : // corresponding boundary conditions. This type is used as key in other maps,
132 : // so it must be hashable. `std::unordered_map` isn't supported by
133 : // `boost::hash`, but `std::map` is.
134 0 : using BoundaryConditionsSignature =
135 : std::map<BoundaryId, elliptic::BoundaryConditionType>;
136 0 : using solver_type = Solver;
137 :
138 0 : struct SolverOptionTag {
139 0 : static std::string name() { return "Solver"; }
140 0 : using type = StoredSolverType;
141 0 : static constexpr Options::String help =
142 : "The linear solver used to invert the Laplace operator. The solver is "
143 : "shared between tensor components with the same type of boundary "
144 : "conditions (Dirichlet-type or Neumann-type).";
145 : };
146 :
147 0 : struct BoundaryConditions {
148 0 : using type = Options::Auto<elliptic::BoundaryConditionType>;
149 0 : static constexpr Options::String help =
150 : "The boundary conditions imposed by the Laplace operator. Specify "
151 : "'Auto' to choose between homogeneous Dirichlet or Neumann boundary "
152 : "conditions automatically, based on the configuration of the the full "
153 : "operator.";
154 : };
155 :
156 0 : using options = tmpl::list<SolverOptionTag, BoundaryConditions>;
157 0 : static constexpr Options::String help =
158 : "Approximate the linear operator with a Laplace operator "
159 : "for every tensor component separately.";
160 :
161 0 : MinusLaplacian() = default;
162 0 : MinusLaplacian(MinusLaplacian&& /*rhs*/) = default;
163 0 : MinusLaplacian& operator=(MinusLaplacian&& /*rhs*/) = default;
164 0 : ~MinusLaplacian() = default;
165 0 : MinusLaplacian(const MinusLaplacian& rhs)
166 : : Base(rhs),
167 : solver_(rhs.clone_solver()),
168 : boundary_condition_type_(rhs.boundary_condition_type_) {}
169 0 : MinusLaplacian& operator=(const MinusLaplacian& rhs) {
170 : Base::operator=(rhs);
171 : solver_ = rhs.clone_solver();
172 : boundary_condition_type_ = rhs.boundary_condition_type_;
173 : return *this;
174 : }
175 :
176 : /// \cond
177 : explicit MinusLaplacian(CkMigrateMessage* m) : Base(m) {}
178 : using PUP::able::register_constructor;
179 : WRAPPED_PUPable_decl_template(MinusLaplacian); // NOLINT
180 : /// \endcond
181 :
182 0 : MinusLaplacian(
183 : StoredSolverType solver,
184 : std::optional<elliptic::BoundaryConditionType> boundary_condition_type)
185 : : solver_(std::move(solver)),
186 : boundary_condition_type_(boundary_condition_type) {}
187 :
188 0 : const Solver& solver() const {
189 : if constexpr (std::is_abstract_v<Solver>) {
190 : return *solver_;
191 : } else {
192 : return solver_;
193 : }
194 : }
195 :
196 : /// The cached solvers for each unique boundary-condition signature. These are
197 : /// only exposed to allow verifying their consistency.
198 1 : const auto& cached_solvers() const { return solvers_; }
199 :
200 : /// Solve the equation \f$Ax=b\f$ by approximating \f$A\f$ with a Laplace
201 : /// operator for every tensor component in \f$x\f$.
202 : template <typename LinearOperator, typename VarsType, typename SourceType,
203 : typename... OperatorArgs>
204 1 : Convergence::HasConverged solve(
205 : gsl::not_null<VarsType*> solution, LinearOperator&& linear_operator,
206 : const SourceType& source,
207 : const std::tuple<OperatorArgs...>& operator_args) const;
208 :
209 0 : void reset() override {
210 : mutable_solver().reset();
211 : // These buffers depend on the operator being solved, so they are reset when
212 : // the operator changes. The other buffers only hold workspace memory so
213 : // they don't need to be reset.
214 : bc_signatures_ = std::nullopt;
215 : solvers_.clear();
216 : boundary_conditions_.clear();
217 : }
218 :
219 : // NOLINTNEXTLINE(google-runtime-references)
220 0 : void pup(PUP::er& p) override {
221 : Base::pup(p);
222 : // Serialize object state
223 : p | solver_;
224 : p | boundary_condition_type_;
225 : // Serialize caches that are possibly expensive to re-create
226 : p | bc_signatures_;
227 : p | solvers_;
228 : // Other properties are memory buffers that don't need to be serialized.
229 : }
230 :
231 0 : std::unique_ptr<Base> get_clone() const override {
232 : return std::make_unique<MinusLaplacian>(*this);
233 : }
234 :
235 : private:
236 : // For each tensor component, get the set of boundary conditions that should
237 : // be imposed by the Laplacian operator on that component, based on the domain
238 : // configuration. These are cached for successive solves of the same operator.
239 : // An empty list means the subdomain has no external boundaries.
240 : template <typename DbTagsList>
241 : const std::vector<BoundaryConditionsSignature>&
242 0 : get_cached_boundary_conditions_signatures(
243 : const db::DataBox<DbTagsList>& box) const {
244 : if (bc_signatures_.has_value()) {
245 : // Use cached value
246 : return *bc_signatures_;
247 : }
248 : bc_signatures_ = std::vector<BoundaryConditionsSignature>{};
249 : const auto& all_boundary_conditions =
250 : db::get<domain::Tags::ExternalBoundaryConditions<Dim>>(box);
251 : const auto collect_bc_signatures = [&bc_signatures = *bc_signatures_,
252 : &override_boundary_condition_type =
253 : boundary_condition_type_,
254 : &all_boundary_conditions](
255 : const BoundaryId& boundary_id) {
256 : if (not bc_signatures.empty() and
257 : bc_signatures[0].find(boundary_id) != bc_signatures[0].end()) {
258 : // We have already handled this block boundary in an earlier
259 : // iteration of the loop. This can happend when the domain is
260 : // h-refined, or when multiple elements in a subdomain face the
261 : // block boundary.
262 : return;
263 : }
264 : // Get the type of boundary condition (Dirichlet or Neumann) for
265 : // each tensor component from the domain configuration
266 : const auto& [block_id, direction] = boundary_id;
267 : const auto original_boundary_condition = dynamic_cast<
268 : const elliptic::BoundaryConditions::BoundaryCondition<Dim>*>(
269 : all_boundary_conditions.at(block_id).at(direction).get());
270 : ASSERT(original_boundary_condition != nullptr,
271 : "The boundary condition in block "
272 : << block_id << ", direction " << direction
273 : << " is not of the expected type "
274 : "'elliptic::BoundaryConditions::BoundaryCondition<"
275 : << Dim << ">'.");
276 : const std::vector<elliptic::BoundaryConditionType> bc_types =
277 : original_boundary_condition->boundary_condition_types();
278 : const size_t num_components = bc_types.size();
279 : if (bc_signatures.empty()) {
280 : bc_signatures.resize(num_components);
281 : } else {
282 : ASSERT(bc_signatures.size() == num_components,
283 : "Unexpected number of boundary-condition types. The "
284 : "boundary condition in block "
285 : << block_id << ", direction " << direction << " returned "
286 : << num_components << " items, but another returned "
287 : << bc_signatures.size()
288 : << " (one for each tensor component).");
289 : }
290 : for (size_t component = 0; component < num_components; ++component) {
291 : bc_signatures[component].emplace(
292 : boundary_id,
293 : override_boundary_condition_type.value_or(bc_types.at(component)));
294 : }
295 : };
296 : // Collect boundary-condition signatures from all elements in the subdomain
297 : const auto& element = db::get<domain::Tags::Element<Dim>>(box);
298 : for (const auto& direction : element.external_boundaries()) {
299 : collect_bc_signatures({element.id().block_id(), direction});
300 : }
301 : const auto& overlap_elements =
302 : db::get<LinearSolver::Schwarz::Tags::Overlaps<
303 : domain::Tags::Element<Dim>, Dim, OptionsGroup>>(box);
304 : for (const auto& [overlap_id, overlap_element] : overlap_elements) {
305 : (void)overlap_id;
306 : for (const auto& direction : overlap_element.external_boundaries()) {
307 : collect_bc_signatures({overlap_element.id().block_id(), direction});
308 : }
309 : }
310 : return *bc_signatures_;
311 : }
312 :
313 : // For a tensor component with the given boundary-condition configuration, get
314 : // the solver and the set of boundary conditions that will be passed to the
315 : // Poisson subdomain operator. The solver is cached for successive solves of
316 : // the same operator. Caching is very important for performance, since the
317 : // solver may build up an explicit matrix representation that is expensive to
318 : // re-create.
319 : std::pair<const Solver&,
320 : const std::unordered_map<BoundaryId, const BoundaryConditionsBase&,
321 : boost::hash<BoundaryId>>&>
322 0 : get_cached_solver_and_boundary_conditions(
323 : const std::vector<BoundaryConditionsSignature>& bc_signatures,
324 : const size_t component) const {
325 : if (bc_signatures.empty()) {
326 : // The subdomain has no external boundaries. We use the original solver.
327 : return {solver(), boundary_conditions_[{}]};
328 : }
329 : // Get the cached solver corresponding to this component's
330 : // boundary-condition configuration
331 : const auto& bc_signature = bc_signatures.at(component);
332 : if (solvers_.find(bc_signature) == solvers_.end()) {
333 : solvers_.emplace(bc_signature, clone_solver());
334 : }
335 : const Solver& cached_solver = [this, &bc_signature]() -> const Solver& {
336 : if constexpr (std::is_abstract_v<Solver>) {
337 : return *solvers_.at(bc_signature);
338 : } else {
339 : return solvers_.at(bc_signature);
340 : }
341 : }();
342 : // Get the cached set of boundary conditions for this component
343 : if (boundary_conditions_.find(bc_signature) == boundary_conditions_.end()) {
344 : auto& bc = boundary_conditions_[bc_signature];
345 : for (const auto& [boundary_id, bc_type] : bc_signature) {
346 : bc.emplace(boundary_id,
347 : bc_type == elliptic::BoundaryConditionType::Dirichlet
348 : ? dirichlet_bc
349 : : neumann_bc);
350 : }
351 : }
352 : return {cached_solver, boundary_conditions_[bc_signature]};
353 : }
354 :
355 0 : Solver& mutable_solver() {
356 : if constexpr (std::is_abstract_v<Solver>) {
357 : return *solver_;
358 : } else {
359 : return solver_;
360 : }
361 : }
362 :
363 0 : StoredSolverType clone_solver() const {
364 : if constexpr (std::is_abstract_v<Solver>) {
365 : return solver_->get_clone();
366 : } else {
367 : return solver_;
368 : }
369 : }
370 :
371 : // The option-constructed solver for the separate Poisson problems. It is
372 : // cloned for each unique boundary-condition configuration of the tensor
373 : // components.
374 0 : StoredSolverType solver_{};
375 0 : std::optional<elliptic::BoundaryConditionType> boundary_condition_type_;
376 :
377 : // These are caches to keep track of the unique boundary-condition
378 : // configurations. Each configuration has its own solver, because the solver
379 : // may cache the operator to speed up repeated solves.
380 : // - The boundary-condition configuration for each tensor component, or an
381 : // empty list if the subdomain has no external boundaries, or `std::nullopt`
382 : // to signal a clean cache.
383 : // NOLINTNEXTLINE(spectre-mutable)
384 : mutable std::optional<std::vector<BoundaryConditionsSignature>>
385 0 : bc_signatures_{};
386 : // - A clone of the `solver_` for each unique boundary-condition configuration
387 : // NOLINTNEXTLINE(spectre-mutable)
388 : mutable std::unordered_map<BoundaryConditionsSignature, StoredSolverType,
389 : boost::hash<BoundaryConditionsSignature>>
390 0 : solvers_{};
391 : // NOLINTNEXTLINE(spectre-mutable)
392 : mutable std::unordered_map<
393 : BoundaryConditionsSignature,
394 : std::unordered_map<BoundaryId, const BoundaryConditionsBase&,
395 : boost::hash<BoundaryId>>,
396 : boost::hash<BoundaryConditionsSignature>>
397 0 : boundary_conditions_{};
398 :
399 : // These are memory buffers that are kept around for repeated solves. Possible
400 : // optimization: Free this memory at appropriate times, e.g. when the element
401 : // has completed a full subdomain solve and goes to the background. In some
402 : // cases the `subdomain_operator_` is never even used again in subsequent
403 : // subdomain solves because it is cached as a matrix (see
404 : // LinearSolver::Serial::ExplicitInverse), so we don't need the memory anymore
405 : // at all.
406 : // NOLINTNEXTLINE(spectre-mutable)
407 0 : mutable SubdomainOperator subdomain_operator_{};
408 : // NOLINTNEXTLINE(spectre-mutable)
409 0 : mutable SubdomainData source_{};
410 : // NOLINTNEXTLINE(spectre-mutable)
411 0 : mutable SubdomainData initial_guess_in_solution_out_{};
412 :
413 : // These boundary condition instances can be re-used for all tensor components
414 : // Note that the _linearized_ boundary conditions are applied by the subdomain
415 : // operator, so the constant (third argument) is irrelevant.
416 0 : const Poisson::BoundaryConditions::Robin<Dim> dirichlet_bc{1., 0., 0.};
417 0 : const Poisson::BoundaryConditions::Robin<Dim> neumann_bc{0., 1., 0.};
418 : };
419 :
420 : namespace detail {
421 : template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
422 : void assign_component(
423 : const gsl::not_null<::LinearSolver::Schwarz::ElementCenteredSubdomainData<
424 : Dim, LhsTagsList>*>
425 : lhs,
426 : const ::LinearSolver::Schwarz::ElementCenteredSubdomainData<
427 : Dim, RhsTagsList>& rhs,
428 : const size_t lhs_component, const size_t rhs_component,
429 : // Use real or imaginary part of the lhs or rhs complex value. This is used
430 : // to apply the Laplacian to the real or imaginary part separately.
431 : const bool lhs_imag = false, const bool rhs_imag = false) {
432 : using LhsValueType =
433 : typename ::LinearSolver::Schwarz::ElementCenteredSubdomainData<
434 : Dim, LhsTagsList>::value_type;
435 : // Possible optimization: Once we have non-owning Variables we can use a view
436 : // into the rhs here instead of copying.
437 : const size_t num_points_element = rhs.element_data.number_of_grid_points();
438 : for (size_t i = 0; i < num_points_element; ++i) {
439 : auto& lhs_i =
440 : lhs->element_data.data()[lhs_component * num_points_element + i];
441 : const auto& rhs_i =
442 : rhs.element_data.data()[rhs_component * num_points_element + i];
443 : const double rhs_i_fundamental = rhs_imag ? imag(rhs_i) : real(rhs_i);
444 : if constexpr (std::is_same_v<LhsValueType, std::complex<double>>) {
445 : if (lhs_imag) {
446 : lhs_i.imag(rhs_i_fundamental);
447 : } else {
448 : lhs_i.real(rhs_i_fundamental);
449 : }
450 : } else {
451 : lhs_i = rhs_i_fundamental;
452 : }
453 : }
454 : for (const auto& [overlap_id, rhs_data] : rhs.overlap_data) {
455 : const size_t num_points_overlap = rhs_data.number_of_grid_points();
456 : // The random-access operation is relatively slow because it computes a
457 : // hash, so it's important for performance to avoid repeating it in every
458 : // iteration of the loop below.
459 : auto& lhs_vars = lhs->overlap_data[overlap_id];
460 : for (size_t i = 0; i < num_points_overlap; ++i) {
461 : auto& lhs_i = lhs_vars.data()[lhs_component * num_points_overlap + i];
462 : const auto& rhs_i =
463 : rhs_data.data()[rhs_component * num_points_overlap + i];
464 : const double rhs_i_fundamental = rhs_imag ? imag(rhs_i) : real(rhs_i);
465 : if constexpr (std::is_same_v<LhsValueType, std::complex<double>>) {
466 : if (lhs_imag) {
467 : lhs_i.imag(rhs_i_fundamental);
468 : } else {
469 : lhs_i.real(rhs_i_fundamental);
470 : }
471 : } else {
472 : lhs_i = rhs_i_fundamental;
473 : }
474 : }
475 : }
476 : }
477 : } // namespace detail
478 :
479 : template <size_t Dim, typename OptionsGroup, typename Solver,
480 : typename LinearSolverRegistrars>
481 : template <typename LinearOperator, typename VarsType, typename SourceType,
482 : typename... OperatorArgs>
483 : Convergence::HasConverged
484 0 : MinusLaplacian<Dim, OptionsGroup, Solver, LinearSolverRegistrars>::solve(
485 : const gsl::not_null<VarsType*> initial_guess_in_solution_out,
486 : LinearOperator&& /*linear_operator*/, const SourceType& source,
487 : const std::tuple<OperatorArgs...>& operator_args) const {
488 : const auto& box = get<0>(operator_args);
489 : // Solve each component of the source variables in turn, assuming the operator
490 : // is a Laplacian. For each component we select either homogeneous Dirichlet
491 : // or Neumann boundary conditions, based on the type of boundary conditions
492 : // imposed by the full operator.
493 : static constexpr size_t num_components =
494 : VarsType::ElementData::number_of_independent_components;
495 : source_.destructive_resize(source);
496 : initial_guess_in_solution_out_.destructive_resize(source);
497 : const std::vector<BoundaryConditionsSignature>& bc_signatures =
498 : get_cached_boundary_conditions_signatures(box);
499 : ASSERT(bc_signatures.empty() or bc_signatures.size() == num_components,
500 : "Expected " << num_components
501 : << " boundary-condition signatures (one per tensor "
502 : "component), but got "
503 : << bc_signatures.size() << ".");
504 : // Possible optimization: elide when num_components is 1
505 : using ValueType = typename VarsType::value_type;
506 : constexpr auto complex_components = []() {
507 : if constexpr (std::is_same_v<ValueType, std::complex<double>>) {
508 : return std::array<bool, 2>{{false, true}};
509 : } else {
510 : return std::array<bool, 1>{{false}};
511 : }
512 : }();
513 : for (size_t component = 0; component < num_components; ++component) {
514 : const auto& [solver, boundary_conditions] =
515 : get_cached_solver_and_boundary_conditions(bc_signatures, component);
516 : // Solve real and imaginary parts separately, if applicable
517 : for (const bool imag_component : complex_components) {
518 : detail::assign_component(make_not_null(&source_), source, 0, component,
519 : false, imag_component);
520 : detail::assign_component(make_not_null(&initial_guess_in_solution_out_),
521 : *initial_guess_in_solution_out, 0, component,
522 : false, imag_component);
523 : solver.solve(make_not_null(&initial_guess_in_solution_out_),
524 : subdomain_operator_, source_,
525 : std::forward_as_tuple(box, boundary_conditions));
526 : detail::assign_component(initial_guess_in_solution_out,
527 : initial_guess_in_solution_out_, component, 0,
528 : imag_component, false);
529 : }
530 : }
531 : return {0, 0};
532 : }
533 :
534 : /// \cond
535 : template <size_t Dim, typename OptionsGroup, typename Solver,
536 : typename LinearSolverRegistrars>
537 : // NOLINTNEXTLINE
538 : PUP::able::PUP_ID MinusLaplacian<Dim, OptionsGroup, Solver,
539 : LinearSolverRegistrars>::my_PUP_ID = 0;
540 : /// \endcond
541 :
542 : } // namespace elliptic::subdomain_preconditioners
|
