Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #include "gmg.hpp"
5 :
6 : #include <mfem.hpp>
7 : #include "linalg/chebyshev.hpp"
8 : #include "linalg/distrelaxation.hpp"
9 : #include "linalg/rap.hpp"
10 : #include "utils/timer.hpp"
11 :
12 : namespace palace
13 : {
14 :
15 : template <typename OperType>
16 0 : GeometricMultigridSolver<OperType>::GeometricMultigridSolver(
17 : MPI_Comm comm, std::unique_ptr<Solver<OperType>> &&coarse_solver,
18 : const std::vector<const Operator *> &P, const std::vector<const Operator *> *G,
19 : int cycle_it, int smooth_it, int cheby_order, double cheby_sf_max, double cheby_sf_min,
20 : bool cheby_4th_kind)
21 0 : : Solver<OperType>(), pc_it(cycle_it), P(P.begin(), P.end()), A(P.size() + 1),
22 0 : dbc_tdof_lists(P.size()), B(P.size() + 1), X(P.size() + 1), Y(P.size() + 1),
23 0 : R(P.size() + 1), use_timer(false)
24 : {
25 : // Configure levels of geometric coarsening. Multigrid vectors will be configured at first
26 : // call to Mult. The multigrid operator size is set based on the finest space dimension.
27 0 : const auto n_levels = P.size() + 1;
28 0 : MFEM_VERIFY(n_levels > 0,
29 : "Empty finite element space hierarchy during multigrid solver setup!");
30 0 : MFEM_VERIFY(!G || G->size() == n_levels,
31 : "Invalid input for distributive relaxation smoother auxiliary space transfer "
32 : "operators (mismatch in number of levels)!");
33 :
34 : // Use the supplied level 0 (coarse) solver.
35 : B[0] = std::move(coarse_solver);
36 :
37 : // Configure level smoothers. Use distributive relaxation smoothing if an auxiliary
38 : // finite element space was provided.
39 0 : for (std::size_t l = 1; l < n_levels; l++)
40 : {
41 0 : if (G)
42 : {
43 0 : const int cheby_smooth_it = 1;
44 0 : B[l] = std::make_unique<DistRelaxationSmoother<OperType>>(
45 : comm, *(*G)[l], smooth_it, cheby_smooth_it, cheby_order, cheby_sf_max,
46 : cheby_sf_min, cheby_4th_kind);
47 : }
48 : else
49 : {
50 0 : const int cheby_smooth_it = smooth_it;
51 0 : if (cheby_4th_kind)
52 : {
53 0 : B[l] = std::make_unique<ChebyshevSmoother<OperType>>(comm, cheby_smooth_it,
54 : cheby_order, cheby_sf_max);
55 : }
56 : else
57 : {
58 0 : B[l] = std::make_unique<ChebyshevSmoother1stKind<OperType>>(
59 : comm, cheby_smooth_it, cheby_order, cheby_sf_max, cheby_sf_min);
60 : }
61 : }
62 : }
63 0 : }
64 :
65 : template <typename OperType>
66 0 : void GeometricMultigridSolver<OperType>::SetOperator(const OperType &op)
67 : {
68 : using ParOperType =
69 : typename std::conditional<std::is_same<OperType, ComplexOperator>::value,
70 : ComplexParOperator, ParOperator>::type;
71 :
72 0 : const auto *mg_op = dynamic_cast<const BaseMultigridOperator<OperType> *>(&op);
73 0 : MFEM_VERIFY(mg_op, "GeometricMultigridSolver requires a MultigridOperator or "
74 : "ComplexMultigridOperator argument provided to SetOperator!");
75 :
76 : const auto n_levels = A.size();
77 0 : MFEM_VERIFY(
78 : mg_op->GetNumLevels() == n_levels &&
79 : (!mg_op->HasAuxiliaryOperators() || mg_op->GetNumAuxiliaryLevels() == n_levels),
80 : "Invalid number of levels for operators in multigrid solver setup!");
81 0 : for (std::size_t l = 0; l < n_levels; l++)
82 : {
83 0 : A[l] = &mg_op->GetOperatorAtLevel(l);
84 0 : MFEM_VERIFY(
85 : A[l]->Width() == A[l]->Height() &&
86 : (n_levels == 1 ||
87 : (A[l]->Height() == ((l < n_levels - 1) ? P[l]->Width() : P[l - 1]->Height()))),
88 : "Invalid operator sizes for GeometricMultigridSolver!");
89 :
90 0 : const auto *PtAP_l = dynamic_cast<const ParOperType *>(&mg_op->GetOperatorAtLevel(l));
91 0 : MFEM_VERIFY(
92 : PtAP_l,
93 : "GeometricMultigridSolver requires ParOperator or ComplexParOperator operators!");
94 0 : if (l < n_levels - 1)
95 : {
96 0 : dbc_tdof_lists[l] = PtAP_l->GetEssentialTrueDofs();
97 : }
98 :
99 0 : auto *dist_smoother = dynamic_cast<DistRelaxationSmoother<OperType> *>(B[l].get());
100 0 : if (dist_smoother)
101 : {
102 0 : MFEM_VERIFY(mg_op->HasAuxiliaryOperators(),
103 : "Distributive relaxation smoother relies on both primary space and "
104 : "auxiliary space operators for multigrid smoothing!");
105 0 : dist_smoother->SetOperators(mg_op->GetOperatorAtLevel(l),
106 : mg_op->GetAuxiliaryOperatorAtLevel(l));
107 : }
108 : else
109 : {
110 0 : B[l]->SetOperator(mg_op->GetOperatorAtLevel(l));
111 : }
112 :
113 0 : X[l].SetSize(A[l]->Height());
114 0 : Y[l].SetSize(A[l]->Height());
115 0 : R[l].SetSize(A[l]->Height());
116 0 : X[l].UseDevice(true);
117 0 : Y[l].UseDevice(true);
118 0 : R[l].UseDevice(true);
119 : }
120 :
121 0 : this->height = op.Height();
122 0 : this->width = op.Width();
123 0 : }
124 :
125 : template <typename OperType>
126 0 : void GeometricMultigridSolver<OperType>::Mult(const VecType &x, VecType &y) const
127 : {
128 : // Initialize.
129 : const auto n_levels = A.size();
130 : MFEM_ASSERT(!this->initial_guess,
131 : "Geometric multigrid solver does not use initial guess!");
132 : MFEM_ASSERT(n_levels > 1 || pc_it == 1,
133 : "Single-level geometric multigrid will not work with multiple iterations!");
134 :
135 : // Apply V-cycle. The initial guess for y is zero'd at the first pre-smooth iteration.
136 0 : X.back() = x;
137 0 : for (int it = 0; it < pc_it; it++)
138 : {
139 0 : VCycle(n_levels - 1, (it > 0));
140 : }
141 0 : y = Y.back();
142 0 : }
143 :
144 : namespace
145 : {
146 :
147 : inline void RealMult(const Operator &op, const Vector &x, Vector &y)
148 : {
149 0 : op.Mult(x, y);
150 : }
151 :
152 : inline void RealMult(const Operator &op, const ComplexVector &x, ComplexVector &y)
153 : {
154 0 : op.Mult(x.Real(), y.Real());
155 0 : op.Mult(x.Imag(), y.Imag());
156 : }
157 :
158 : inline void RealMultTranspose(const Operator &op, const Vector &x, Vector &y)
159 : {
160 0 : op.MultTranspose(x, y);
161 : }
162 :
163 : inline void RealMultTranspose(const Operator &op, const ComplexVector &x, ComplexVector &y)
164 : {
165 0 : op.MultTranspose(x.Real(), y.Real());
166 0 : op.MultTranspose(x.Imag(), y.Imag());
167 : }
168 :
169 : } // namespace
170 :
171 : template <typename OperType>
172 0 : void GeometricMultigridSolver<OperType>::VCycle(int l, bool initial_guess) const
173 : {
174 : // Pre-smooth, with zero initial guess (Y = 0 set inside). This is the coarse solve at
175 : // level 0. Important to note that the smoothers must respect the initial guess flag
176 : // correctly (given X, Y, compute Y <- Y + B (X - A Y)) .
177 0 : B[l]->SetInitialGuess(initial_guess);
178 0 : if (l == 0)
179 : {
180 0 : BlockTimer bt(Timer::KSP_COARSE_SOLVE, use_timer);
181 0 : B[l]->Mult(X[l], Y[l]);
182 : return;
183 0 : }
184 0 : B[l]->Mult2(X[l], Y[l], R[l]);
185 :
186 : // Compute residual.
187 0 : A[l]->Mult(Y[l], R[l]);
188 0 : linalg::AXPBY(1.0, X[l], -1.0, R[l]);
189 :
190 : // Coarse grid correction.
191 0 : RealMultTranspose(*P[l - 1], R[l], X[l - 1]);
192 0 : if (dbc_tdof_lists[l - 1])
193 : {
194 0 : linalg::SetSubVector(X[l - 1], *dbc_tdof_lists[l - 1], 0.0);
195 : }
196 0 : VCycle(l - 1, false);
197 :
198 : // Prolongate and add.
199 0 : RealMult(*P[l - 1], Y[l - 1], R[l]);
200 0 : Y[l] += R[l];
201 :
202 : // Post-smooth, with nonzero initial guess.
203 0 : B[l]->SetInitialGuess(true);
204 0 : B[l]->MultTranspose2(X[l], Y[l], R[l]);
205 : }
206 :
207 : template class GeometricMultigridSolver<Operator>;
208 : template class GeometricMultigridSolver<ComplexOperator>;
209 :
210 : } // namespace palace
|
