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