Line data Source code
1 : //
2 : // Class PCG
3 : // Preconditioned Conjugate Gradient solver algorithm
4 : //
5 :
6 : #ifndef IPPL_PCG_H
7 : #define IPPL_PCG_H
8 :
9 : #include "Preconditioner.h"
10 : #include "SolverAlgorithm.h"
11 :
12 : namespace ippl {
13 : template <typename OperatorRet, typename LowerRet, typename UpperRet, typename UpperLowerRet,
14 : typename InverseDiagRet, typename DiagRet, typename FieldLHS,
15 : typename FieldRHS = FieldLHS>
16 : class CG : public SolverAlgorithm<FieldLHS, FieldRHS> {
17 : using Base = SolverAlgorithm<FieldLHS, FieldRHS>;
18 : typedef typename Base::lhs_type::value_type T;
19 :
20 : public:
21 : using typename Base::lhs_type, typename Base::rhs_type;
22 : using OperatorF = std::function<OperatorRet(lhs_type)>;
23 : using LowerF = std::function<LowerRet(lhs_type)>;
24 : using UpperF = std::function<UpperRet(lhs_type)>;
25 : using UpperLowerF = std::function<UpperLowerRet(lhs_type)>;
26 : using InverseDiagF = std::function<InverseDiagRet(lhs_type)>;
27 : using DiagF = std::function<DiagRet(lhs_type)>;
28 :
29 0 : virtual ~CG() = default;
30 :
31 : /*!
32 : * Sets the differential operator for the conjugate gradient algorithm
33 : * @param op A function that returns OpRet and takes a field of the LHS type
34 : */
35 0 : virtual void setOperator(OperatorF op) { op_m = std::move(op); }
36 0 : virtual void setPreconditioner(
37 : [[maybe_unused]] OperatorF&& op, // Operator passed to chebyshev and newton
38 : [[maybe_unused]] LowerF&& lower, // Operator passed to 2-step gauss-seidel and ssor
39 : [[maybe_unused]] UpperF&& upper, // Operator passed to 2-step gauss-seidel and ssor
40 : [[maybe_unused]] UpperLowerF&&
41 : upper_and_lower, // Operator passed to 2-step gauss-seidel
42 : [[maybe_unused]] InverseDiagF&&
43 : inverse_diagonal, // Operator passed to jacobi, 2-step gauss-seidel and ssor
44 : [[maybe_unused]] DiagF&& diagonal, // Operator passed to SSOR
45 : [[maybe_unused]] double alpha, // smallest eigenvalue of the operator
46 : [[maybe_unused]] double beta, // largest eigenvalue of the operator
47 : [[maybe_unused]] std::string preconditioner_type =
48 : "", // Name of the preconditioner that should be used
49 : [[maybe_unused]] int level =
50 : 5, // This is a dummy default parameter, actual default parameter should be
51 : // set in main
52 : [[maybe_unused]] int degree =
53 : 31, // This is a dummy default parameter, actual default parameter should
54 : // be set in main
55 : [[maybe_unused]] int richardson_iterations =
56 : 1, // This is a dummy default parameter, actual default
57 : // parameter should be set in main
58 : [[maybe_unused]] int inner =
59 : 5, // This is a dummy default parameter, actual default parameter should be
60 : // set in main
61 : [[maybe_unused]] int outer =
62 : 1, // This is a dummy default parameter, actual default parameter should be
63 : [[maybe_unused]] double omega =
64 : 1 // This is a dummy default parameter, actual default parameter should be
65 : // set in main
66 0 : ) {}
67 : /*!
68 : * Query how many iterations were required to obtain the solution
69 : * the last time this solver was used
70 : * @return Iteration count of last solve
71 : */
72 0 : virtual int getIterationCount() { return iterations_m; }
73 :
74 0 : virtual void operator()(lhs_type& lhs, rhs_type& rhs,
75 : const ParameterList& params) override {
76 0 : constexpr unsigned Dim = lhs_type::dim;
77 0 : typename lhs_type::Mesh_t mesh = lhs.get_mesh();
78 0 : typename lhs_type::Layout_t layout = lhs.getLayout();
79 :
80 0 : iterations_m = 0;
81 0 : const int maxIterations = params.get<int>("max_iterations");
82 :
83 : // Variable names mostly based on description in
84 : // https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf
85 0 : lhs_type r(mesh, layout);
86 0 : lhs_type d(mesh, layout);
87 :
88 : using bc_type = BConds<lhs_type, Dim>;
89 0 : bc_type lhsBCs = lhs.getFieldBC();
90 0 : bc_type bc;
91 :
92 0 : bool allFacesPeriodic = true;
93 0 : for (unsigned int i = 0; i < 2 * Dim; ++i) {
94 0 : FieldBC bcType = lhsBCs[i]->getBCType();
95 0 : if (bcType == PERIODIC_FACE) {
96 : // If the LHS has periodic BCs, so does the residue
97 0 : bc[i] = std::make_shared<PeriodicFace<lhs_type>>(i);
98 0 : } else if (bcType & CONSTANT_FACE) {
99 : // If the LHS has constant BCs, the residue is zero on the BCs
100 : // Bitwise AND with CONSTANT_FACE will succeed for ZeroFace or ConstantFace
101 0 : bc[i] = std::make_shared<ZeroFace<lhs_type>>(i);
102 0 : allFacesPeriodic = false;
103 : } else {
104 0 : throw IpplException("PCG::operator()",
105 : "Only periodic or constant BCs for LHS supported.");
106 : return;
107 : }
108 : }
109 :
110 0 : r = rhs - op_m(lhs);
111 0 : d = r.deepCopy();
112 0 : d.setFieldBC(bc);
113 :
114 0 : T delta1 = innerProduct(r, d);
115 0 : T delta0 = delta1;
116 0 : residueNorm = std::sqrt(delta1);
117 0 : const T tolerance = params.get<T>("tolerance") * norm(rhs);
118 :
119 0 : lhs_type q(mesh, layout);
120 :
121 0 : while (iterations_m < maxIterations && residueNorm > tolerance) {
122 0 : q = op_m(d);
123 :
124 0 : T alpha = delta1 / innerProduct(d, q);
125 0 : lhs = lhs + alpha * d;
126 :
127 : // The exact residue is given by
128 : // r = rhs - op_m(lhs);
129 : // This correction is generally not used in practice because
130 : // applying the Laplacian is computationally expensive and
131 : // the correction does not have a significant effect on accuracy;
132 : // in some implementations, the correction may be applied every few
133 : // iterations to offset accumulated floating point errors
134 0 : r = r - alpha * q;
135 0 : delta0 = delta1;
136 0 : delta1 = innerProduct(r, r);
137 0 : T beta = delta1 / delta0;
138 :
139 0 : residueNorm = std::sqrt(delta1);
140 0 : d = r + beta * d;
141 0 : ++iterations_m;
142 : }
143 :
144 0 : if (allFacesPeriodic) {
145 0 : T avg = lhs.getVolumeAverage();
146 0 : lhs = lhs - avg;
147 : }
148 0 : }
149 :
150 0 : virtual T getResidue() const { return residueNorm; }
151 :
152 : protected:
153 : OperatorF op_m;
154 : T residueNorm = 0;
155 : int iterations_m = 0;
156 : };
157 :
158 : template <typename OperatorRet, typename LowerRet, typename UpperRet, typename UpperLowerRet,
159 : typename InverseDiagRet, typename DiagRet, typename FieldLHS,
160 : typename FieldRHS = FieldLHS>
161 : class PCG : public CG<OperatorRet, LowerRet, UpperRet, UpperLowerRet, InverseDiagRet, DiagRet,
162 : FieldLHS, FieldRHS> {
163 : using Base = SolverAlgorithm<FieldLHS, FieldRHS>;
164 : typedef typename Base::lhs_type::value_type T;
165 :
166 : public:
167 : using typename Base::lhs_type, typename Base::rhs_type;
168 : using OperatorF = std::function<OperatorRet(lhs_type)>;
169 : using LowerF = std::function<LowerRet(lhs_type)>;
170 : using UpperF = std::function<UpperRet(lhs_type)>;
171 : using UpperLowerF = std::function<UpperLowerRet(lhs_type)>;
172 : using InverseDiagF = std::function<InverseDiagRet(lhs_type)>;
173 : using DiagF = std::function<DiagRet(lhs_type)>;
174 :
175 0 : PCG()
176 : : CG<OperatorRet, LowerRet, UpperRet, UpperLowerRet, InverseDiagRet, DiagRet, FieldLHS,
177 : FieldRHS>()
178 0 : , preconditioner_m(nullptr){};
179 :
180 : /*!
181 : * Sets the differential operator for the conjugate gradient algorithm
182 : * @param op A function that returns OpRet and takes a field of the LHS type
183 : */
184 0 : void setPreconditioner(
185 : OperatorF&& op, // Operator passed to chebyshev and newton
186 : LowerF&& lower, // Operator passed to 2-step gauss-seidel
187 : UpperF&& upper, // Operator passed to 2-step gauss-seidel
188 : UpperLowerF&& upper_and_lower, // Operator passed to 2-step gauss-seidel
189 : InverseDiagF&& inverse_diagonal, // Operator passed to jacobi and 2-step gauss-seidel
190 : DiagF&& diagonal, // Operator passed to ssor
191 : double alpha, // smallest eigenvalue of the operator
192 : double beta, // largest eigenvalue of the operator
193 : std::string preconditioner_type = "", // Name of the preconditioner that should be used
194 : int level = 5, // This is a dummy default parameter, actual default parameter should be
195 : // set in main
196 : int degree = 31, // This is a dummy default parameter, actual default parameter should
197 : // be set in main
198 : int richardson_iterations = 4, // This is a dummy default parameter, actual default
199 : // parameter should be set in main
200 : int inner = 2, // This is a dummy default parameter, actual default parameter should be
201 : // set in main
202 : int outer = 2, // This is a dummy default parameter, actual default parameter should be
203 : // set in main
204 : double omega = 1.57079632679 // This is a dummy default parameter, actual default
205 : // parameter should be set in main
206 : // default = pi/2 as this was found optimal during hyperparameter scan for test case
207 : // (see https://amas.web.psi.ch/people/aadelmann/ETH-Accel-Lecture-1/projectscompleted/cse/BSc-mbolliger.pdf)
208 : ) override {
209 0 : if (preconditioner_type == "jacobi") {
210 : // Turn on damping parameter
211 : /*
212 : double w = 2.0 / ((alpha + beta));
213 : preconditioner_m = std::move(std::make_unique<jacobi_preconditioner<FieldLHS ,
214 : InverseDiagF>>(std::move(inverse_diagonal), w));
215 : */
216 0 : preconditioner_m =
217 0 : std::move(std::make_unique<jacobi_preconditioner<FieldLHS, InverseDiagF>>(
218 0 : std::move(inverse_diagonal)));
219 0 : } else if (preconditioner_type == "newton") {
220 0 : preconditioner_m = std::move(
221 0 : std::make_unique<polynomial_newton_preconditioner<FieldLHS, OperatorF>>(
222 0 : std::move(op), alpha, beta, level, 1e-3));
223 0 : } else if (preconditioner_type == "chebyshev") {
224 0 : preconditioner_m = std::move(
225 0 : std::make_unique<polynomial_chebyshev_preconditioner<FieldLHS, OperatorF>>(
226 0 : std::move(op), alpha, beta, degree, 1e-3));
227 0 : } else if (preconditioner_type == "richardson") {
228 0 : preconditioner_m =
229 0 : std::move(std::make_unique<
230 0 : richardson_preconditioner<FieldLHS, UpperLowerF, InverseDiagF>>(
231 0 : std::move(upper_and_lower), std::move(inverse_diagonal),
232 : richardson_iterations));
233 0 : } else if (preconditioner_type == "gauss-seidel") {
234 0 : preconditioner_m = std::move(
235 0 : std::make_unique<gs_preconditioner<FieldLHS, LowerF, UpperF, InverseDiagF>>(
236 0 : std::move(lower), std::move(upper), std::move(inverse_diagonal), inner,
237 : outer));
238 0 : } else if (preconditioner_type == "ssor") {
239 0 : preconditioner_m =
240 0 : std::move(std::make_unique<
241 0 : ssor_preconditioner<FieldLHS, LowerF, UpperF, InverseDiagF, DiagF>>(
242 0 : std::move(lower), std::move(upper), std::move(inverse_diagonal),
243 0 : std::move(diagonal), inner, outer, omega));
244 : } else {
245 0 : preconditioner_m = std::move(std::make_unique<preconditioner<FieldLHS>>());
246 : }
247 0 : }
248 :
249 0 : void operator()(lhs_type& lhs, rhs_type& rhs, const ParameterList& params) override {
250 0 : constexpr unsigned Dim = lhs_type::dim;
251 :
252 0 : if (preconditioner_m == nullptr) {
253 0 : throw IpplException("PCG::operator()",
254 : "Preconditioner has not been set for PCG solver");
255 : }
256 :
257 0 : typename lhs_type::Mesh_t mesh = lhs.get_mesh();
258 0 : typename lhs_type::Layout_t layout = lhs.getLayout();
259 :
260 0 : this->iterations_m = 0;
261 0 : const int maxIterations = params.get<int>("max_iterations");
262 :
263 : // Variable names mostly based on description in
264 : // https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf
265 0 : lhs_type r(mesh, layout);
266 0 : lhs_type d(mesh, layout);
267 0 : lhs_type s(mesh, layout);
268 0 : lhs_type q(mesh, layout);
269 :
270 0 : preconditioner_m->init_fields(lhs);
271 :
272 : using bc_type = BConds<lhs_type, Dim>;
273 0 : bc_type lhsBCs = lhs.getFieldBC();
274 0 : bc_type bc;
275 :
276 0 : bool allFacesPeriodic = true;
277 0 : for (unsigned int i = 0; i < 2 * Dim; ++i) {
278 0 : FieldBC bcType = lhsBCs[i]->getBCType();
279 0 : if (bcType == PERIODIC_FACE) {
280 : // If the LHS has periodic BCs, so does the residue
281 0 : bc[i] = std::make_shared<PeriodicFace<lhs_type>>(i);
282 0 : } else if (bcType & CONSTANT_FACE) {
283 : // If the LHS has constant BCs, the residue is zero on the BCs
284 : // Bitwise AND with CONSTANT_FACE will succeed for ZeroFace or ConstantFace
285 0 : bc[i] = std::make_shared<ZeroFace<lhs_type>>(i);
286 0 : allFacesPeriodic = false;
287 : } else {
288 0 : throw IpplException("PCG::operator()",
289 : "Only periodic or constant BCs for LHS supported.");
290 : return;
291 : }
292 : }
293 :
294 0 : r = rhs - this->op_m(lhs);
295 0 : d = preconditioner_m->operator()(r);
296 0 : d.setFieldBC(bc);
297 :
298 0 : T delta1 = innerProduct(r, d);
299 0 : T delta0 = delta1;
300 0 : this->residueNorm = Kokkos::sqrt(Kokkos::abs(delta1));
301 0 : const T tolerance = params.get<T>("tolerance") * delta1;
302 :
303 0 : while (this->iterations_m<maxIterations&& this->residueNorm> tolerance) {
304 0 : q = this->op_m(d);
305 0 : T alpha = delta1 / innerProduct(d, q);
306 0 : lhs = lhs + alpha * d;
307 :
308 : // The exact residue is given by
309 : // r = rhs - BaseCG::op_m(lhs);
310 : // This correction is generally not used in practice because
311 : // applying the Laplacian is computationally expensive and
312 : // the correction does not have a significant effect on accuracy;
313 : // in some implementations, the correction may be applied every few
314 : // iterations to offset accumulated floating point errors
315 0 : r = r - alpha * q;
316 0 : s = preconditioner_m->operator()(r);
317 :
318 0 : delta0 = delta1;
319 0 : delta1 = innerProduct(r, s);
320 :
321 0 : T beta = delta1 / delta0;
322 0 : this->residueNorm = Kokkos::sqrt(Kokkos::abs(delta1));
323 :
324 0 : d = s + beta * d;
325 0 : ++this->iterations_m;
326 : }
327 :
328 0 : if (allFacesPeriodic) {
329 0 : T avg = lhs.getVolumeAverage();
330 0 : lhs = lhs - avg;
331 : }
332 0 : }
333 :
334 : protected:
335 : std::unique_ptr<preconditioner<FieldLHS>> preconditioner_m;
336 : };
337 :
338 : }; // namespace ippl
339 :
340 : #endif
|
