Line data Source code
1 : //
2 : // General pre-conditioners for the pre-conditioned Conjugate Gradient Solver
3 : // such as Jacobi, Polynomial Newton, Polynomial Chebyshev, Richardson, Gauss-Seidel
4 : // provides a function operator() that returns the preconditioned field
5 : //
6 :
7 : #ifndef IPPL_PRECONDITIONER_H
8 : #define IPPL_PRECONDITIONER_H
9 :
10 : #include "Expression/IpplOperations.h" // get the function apply()
11 :
12 : // Expands to a lambda that acts as a wrapper for a differential operator
13 : // fun: the function for which to create the wrapper, such as ippl::laplace
14 : // type: the argument type, which should match the LHS type for the solver
15 : #define IPPL_SOLVER_OPERATOR_WRAPPER(fun, type) \
16 : [](type arg) { \
17 : return fun(arg); \
18 : }
19 :
20 : namespace ippl {
21 : template <typename Field>
22 : struct preconditioner {
23 : constexpr static unsigned Dim = Field::dim;
24 : using mesh_type = typename Field::Mesh_t;
25 : using layout_type = typename Field::Layout_t;
26 :
27 0 : preconditioner()
28 0 : : type_m("Identity") {}
29 :
30 0 : preconditioner(std::string name)
31 0 : : type_m(name) {}
32 :
33 0 : virtual ~preconditioner() = default;
34 :
35 : // Placeholder for the function operator, actually implemented in the derived classes
36 0 : virtual Field operator()(Field& u) {
37 0 : Field res = u.deepCopy();
38 0 : return res;
39 : }
40 :
41 : // Placeholder for setting additional fields, actually implemented in the derived classes
42 0 : virtual void init_fields(Field& b) {
43 0 : Field res = b.deepCopy();
44 0 : return;
45 0 : }
46 :
47 : std::string get_type() { return type_m; };
48 :
49 : protected:
50 : std::string type_m;
51 : };
52 :
53 : /*!
54 : * Jacobi preconditioner
55 : * M = w*diag{A} // w is a damping factor
56 : */
57 : template <typename Field, typename InvDiagF>
58 : struct jacobi_preconditioner : public preconditioner<Field> {
59 : constexpr static unsigned Dim = Field::dim;
60 : using mesh_type = typename Field::Mesh_t;
61 : using layout_type = typename Field::Layout_t;
62 :
63 0 : jacobi_preconditioner(InvDiagF&& inverse_diagonal, double w = 1.0)
64 : : preconditioner<Field>("jacobi")
65 0 : , w_m(w) {
66 0 : inverse_diagonal_m = std::move(inverse_diagonal);
67 0 : }
68 :
69 0 : Field operator()(Field& u) override {
70 0 : mesh_type& mesh = u.get_mesh();
71 0 : layout_type& layout = u.getLayout();
72 0 : Field res(mesh, layout);
73 :
74 0 : res = inverse_diagonal_m(u);
75 0 : res = w_m * res;
76 0 : return res;
77 0 : }
78 :
79 : protected:
80 : InvDiagF inverse_diagonal_m;
81 : double w_m; // Damping factor
82 : };
83 :
84 : /*!
85 : * Polynomial Newton Preconditioner
86 : * Computes iteratively approximations for A^-1
87 : */
88 : template <typename Field, typename OperatorF>
89 : struct polynomial_newton_preconditioner : public preconditioner<Field> {
90 : constexpr static unsigned Dim = Field::dim;
91 : using mesh_type = typename Field::Mesh_t;
92 : using layout_type = typename Field::Layout_t;
93 :
94 0 : polynomial_newton_preconditioner(OperatorF&& op, double alpha, double beta,
95 : unsigned int max_level = 6, double zeta = 1e-3,
96 : double* eta = nullptr)
97 : : preconditioner<Field>("polynomial_newton")
98 0 : , alpha_m(alpha)
99 0 : , beta_m(beta)
100 0 : , level_m(max_level)
101 0 : , zeta_m(zeta)
102 0 : , eta_m(eta) {
103 0 : op_m = std::move(op);
104 0 : }
105 :
106 : // Memory management is needed because of runtime defined eta_m
107 0 : ~polynomial_newton_preconditioner() {
108 0 : if (eta_m != nullptr) {
109 0 : delete[] eta_m;
110 0 : eta_m = nullptr;
111 : }
112 0 : }
113 :
114 : polynomial_newton_preconditioner(const polynomial_newton_preconditioner& other)
115 : : preconditioner<Field>("polynomial_newton")
116 : , level_m(other.level_m)
117 : , alpha_m(other.alpha_m)
118 : , beta_m(other.beta_m)
119 : , zeta_m(other.zeta_m)
120 : , eta_m(other.eta_m) {
121 : op_m = std::move(other.op_m);
122 : }
123 :
124 : polynomial_newton_preconditioner& operator=(const polynomial_newton_preconditioner& other) {
125 : return *this = polynomial_newton_preconditioner(other);
126 : }
127 :
128 0 : Field recursive_preconditioner(Field& u, unsigned int level) {
129 0 : mesh_type& mesh = u.get_mesh();
130 0 : layout_type& layout = u.getLayout();
131 : // Define etas if not defined yet
132 0 : if (eta_m == nullptr) {
133 : // Precompute the etas for later use
134 0 : eta_m = new double[level_m + 1];
135 0 : eta_m[0] = 2.0 / ((alpha_m + beta_m) * (1.0 + zeta_m));
136 0 : if (level_m > 0) {
137 0 : eta_m[1] =
138 : 2.0
139 0 : / (1.0 + 2 * alpha_m * eta_m[0] - alpha_m * eta_m[0] * alpha_m * eta_m[0]);
140 : }
141 0 : for (unsigned int i = 2; i < level_m + 1; ++i) {
142 0 : eta_m[i] = 2.0 / (1.0 + 2 * eta_m[i - 1] - eta_m[i - 1] * eta_m[i - 1]);
143 : }
144 : }
145 :
146 0 : Field res(mesh, layout);
147 : // Base case
148 0 : if (level == 0) {
149 0 : res = eta_m[0] * u;
150 0 : return res;
151 : }
152 : // Recursive case
153 0 : Field PAPr(mesh, layout);
154 0 : Field Pr(mesh, layout);
155 :
156 0 : Pr = recursive_preconditioner(u, level - 1);
157 0 : PAPr = op_m(Pr);
158 0 : PAPr = recursive_preconditioner(PAPr, level - 1);
159 0 : res = eta_m[level] * (2.0 * Pr - PAPr);
160 0 : return res;
161 0 : }
162 :
163 0 : Field operator()(Field& u) override { return recursive_preconditioner(u, level_m); }
164 :
165 : protected:
166 : OperatorF op_m; // Operator to be preconditioned
167 : double alpha_m; // Smallest Eigenvalue
168 : double beta_m; // Largest Eigenvalue
169 : unsigned int level_m; // Number of recursive calls
170 : double zeta_m; // smallest (alpha + beta) is multiplied by (1+zeta) to avoid clustering of
171 : // Eigenvalues
172 : double* eta_m = nullptr; // Size is determined at runtime
173 : };
174 :
175 : /*!
176 : * Polynomial Chebyshev Preconditioner
177 : * Computes iteratively approximations for A^-1
178 : */
179 : template <typename Field, typename OperatorF>
180 : struct polynomial_chebyshev_preconditioner : public preconditioner<Field> {
181 : constexpr static unsigned Dim = Field::dim;
182 : using mesh_type = typename Field::Mesh_t;
183 : using layout_type = typename Field::Layout_t;
184 :
185 0 : polynomial_chebyshev_preconditioner(OperatorF&& op, double alpha, double beta,
186 : unsigned int degree = 63, double zeta = 1e-3)
187 : : preconditioner<Field>("polynomial_chebyshev")
188 0 : , alpha_m(alpha)
189 0 : , beta_m(beta)
190 0 : , degree_m(degree)
191 0 : , zeta_m(zeta)
192 0 : , rho_m(nullptr) {
193 0 : op_m = std::move(op);
194 0 : }
195 :
196 : // Memory management is needed because of runtime defined rho_m
197 0 : ~polynomial_chebyshev_preconditioner() {
198 0 : if (rho_m != nullptr) {
199 0 : delete[] rho_m;
200 0 : rho_m = nullptr;
201 : }
202 0 : }
203 :
204 : polynomial_chebyshev_preconditioner(const polynomial_chebyshev_preconditioner& other)
205 : : preconditioner<Field>("polynomial_chebyshev")
206 : , degree_m(other.degree_m)
207 : , theta_m(other.theta_m)
208 : , sigma_m(other.sigma_m)
209 : , delta_m(other.delta_m)
210 : , alpha_m(other.delta_m)
211 : , beta_m(other.delta_m)
212 : , zeta_m(other.zeta_m)
213 : , rho_m(other.rho_m) {
214 : op_m = std::move(other.op_m);
215 : }
216 :
217 : polynomial_chebyshev_preconditioner& operator=(
218 : const polynomial_chebyshev_preconditioner& other) {
219 : return *this = polynomial_chebyshev_preconditioner(other);
220 : }
221 :
222 0 : Field operator()(Field& r) override {
223 0 : mesh_type& mesh = r.get_mesh();
224 0 : layout_type& layout = r.getLayout();
225 :
226 0 : Field res(mesh, layout);
227 0 : Field x(mesh, layout);
228 0 : Field x_old(mesh, layout);
229 0 : Field A(mesh, layout);
230 0 : Field z(mesh, layout);
231 :
232 : // Precompute the coefficients if not done yet
233 0 : if (rho_m == nullptr) {
234 : // Start precomputing the coefficients
235 0 : theta_m = (beta_m + alpha_m) / 2.0 * (1.0 + zeta_m);
236 0 : delta_m = (beta_m - alpha_m) / 2.0;
237 0 : sigma_m = theta_m / delta_m;
238 :
239 0 : rho_m = new double[degree_m + 1];
240 0 : rho_m[0] = 1.0 / sigma_m;
241 0 : for (unsigned int i = 1; i < degree_m + 1; ++i) {
242 0 : rho_m[i] = 1.0 / (2.0 * sigma_m - rho_m[i - 1]);
243 : }
244 : } // End of precomputing the coefficients
245 :
246 0 : res = r.deepCopy();
247 :
248 0 : x_old = r / theta_m;
249 0 : A = op_m(r);
250 0 : x = 2.0 * rho_m[1] / delta_m * (2.0 * r - A / theta_m);
251 :
252 0 : if (degree_m == 0) {
253 0 : return x_old;
254 : }
255 :
256 0 : if (degree_m == 1) {
257 0 : return x;
258 : }
259 0 : for (unsigned int i = 2; i < degree_m + 1; ++i) {
260 0 : A = op_m(x);
261 0 : z = 2.0 / delta_m * (r - A);
262 0 : res = rho_m[i] * (2 * sigma_m * x - rho_m[i - 1] * x_old + z);
263 0 : x_old = x.deepCopy();
264 0 : x = res.deepCopy();
265 : }
266 0 : return res;
267 0 : }
268 :
269 : protected:
270 : OperatorF op_m;
271 : double alpha_m;
272 : double beta_m;
273 : double delta_m;
274 : double theta_m;
275 : double sigma_m;
276 : unsigned degree_m;
277 : double zeta_m;
278 : double* rho_m = nullptr; // Size is determined at runtime
279 : };
280 :
281 : /*!
282 : * Richardson preconditioner
283 : */
284 : template <typename Field, typename UpperAndLowerF, typename InvDiagF>
285 : struct richardson_preconditioner : public preconditioner<Field> {
286 : constexpr static unsigned Dim = Field::dim;
287 : using mesh_type = typename Field::Mesh_t;
288 : using layout_type = typename Field::Layout_t;
289 :
290 0 : richardson_preconditioner(UpperAndLowerF&& upper_and_lower, InvDiagF&& inverse_diagonal,
291 : unsigned innerloops = 5)
292 : : preconditioner<Field>("Richardson")
293 0 : , innerloops_m(innerloops) {
294 0 : upper_and_lower_m = std::move(upper_and_lower);
295 0 : inverse_diagonal_m = std::move(inverse_diagonal);
296 0 : }
297 :
298 0 : Field operator()(Field& r) override {
299 0 : mesh_type& mesh = r.get_mesh();
300 0 : layout_type& layout = r.getLayout();
301 0 : Field g(mesh, layout);
302 :
303 0 : g = 0;
304 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
305 0 : ULg_m = upper_and_lower_m(g);
306 0 : g = r - ULg_m;
307 :
308 : // The inverse diagonal is applied to the
309 : // vector itself to return the result usually.
310 : // However, the operator for FEM already
311 : // returns the result of inv_diag * itself
312 : // due to the matrix-free evaluation.
313 : // Therefore, we need this if to differentiate
314 : // the two cases.
315 : if constexpr (std::is_same_v<InvDiagF, double>) {
316 : g = inverse_diagonal_m(g) * g;
317 : } else {
318 0 : g = inverse_diagonal_m(g);
319 : }
320 : }
321 0 : return g;
322 0 : }
323 :
324 0 : void init_fields(Field& b) override {
325 0 : layout_type& layout = b.getLayout();
326 0 : mesh_type& mesh = b.get_mesh();
327 :
328 0 : ULg_m = Field(mesh, layout);
329 0 : }
330 :
331 : protected:
332 : UpperAndLowerF upper_and_lower_m;
333 : InvDiagF inverse_diagonal_m;
334 : unsigned innerloops_m;
335 : Field ULg_m;
336 : };
337 :
338 : /*!
339 : * 2-step Gauss-Seidel preconditioner
340 : */
341 : template <typename Field, typename LowerF, typename UpperF, typename InvDiagF>
342 : struct gs_preconditioner : public preconditioner<Field> {
343 : constexpr static unsigned Dim = Field::dim;
344 : using mesh_type = typename Field::Mesh_t;
345 : using layout_type = typename Field::Layout_t;
346 :
347 0 : gs_preconditioner(LowerF&& lower, UpperF&& upper, InvDiagF&& inverse_diagonal,
348 : unsigned innerloops, unsigned outerloops)
349 : : preconditioner<Field>("Gauss-Seidel")
350 0 : , innerloops_m(innerloops)
351 0 : , outerloops_m(outerloops) {
352 0 : lower_m = std::move(lower);
353 0 : upper_m = std::move(upper);
354 0 : inverse_diagonal_m = std::move(inverse_diagonal);
355 0 : }
356 :
357 0 : Field operator()(Field& b) override {
358 0 : layout_type& layout = b.getLayout();
359 0 : mesh_type& mesh = b.get_mesh();
360 :
361 0 : Field x(mesh, layout);
362 :
363 0 : x = 0; // Initial guess
364 :
365 0 : for (unsigned int k = 0; k < outerloops_m; ++k) {
366 0 : UL_m = upper_m(x);
367 0 : r_m = b - UL_m;
368 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
369 0 : UL_m = lower_m(x);
370 0 : x = r_m - UL_m;
371 : // The inverse diagonal is applied to the
372 : // vector itself to return the result usually.
373 : // However, the operator for FEM already
374 : // returns the result of inv_diag * itself
375 : // due to the matrix-free evaluation.
376 : // Therefore, we need this if to differentiate
377 : // the two cases.
378 : if constexpr (std::is_same_v<InvDiagF, double>) {
379 : x = inverse_diagonal_m(x) * x;
380 : } else {
381 0 : x = inverse_diagonal_m(x);
382 : }
383 : }
384 0 : UL_m = lower_m(x);
385 0 : r_m = b - UL_m;
386 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
387 0 : UL_m = upper_m(x);
388 0 : x = r_m - UL_m;
389 : // The inverse diagonal is applied to the
390 : // vector itself to return the result usually.
391 : // However, the operator for FEM already
392 : // returns the result of inv_diag * itself
393 : // due to the matrix-free evaluation.
394 : // Therefore, we need this if to differentiate
395 : // the two cases.
396 : if constexpr (std::is_same_v<InvDiagF, double>) {
397 : x = inverse_diagonal_m(x) * x;
398 : } else {
399 0 : x = inverse_diagonal_m(x);
400 : }
401 : }
402 : }
403 0 : return x;
404 0 : }
405 :
406 0 : void init_fields(Field& b) override {
407 0 : layout_type& layout = b.getLayout();
408 0 : mesh_type& mesh = b.get_mesh();
409 :
410 0 : UL_m = Field(mesh, layout);
411 0 : r_m = Field(mesh, layout);
412 0 : }
413 :
414 : protected:
415 : LowerF lower_m;
416 : UpperF upper_m;
417 : InvDiagF inverse_diagonal_m;
418 : unsigned innerloops_m;
419 : unsigned outerloops_m;
420 : Field UL_m;
421 : Field r_m;
422 : };
423 :
424 : /*!
425 : * Symmetric successive over-relaxation
426 : */
427 : template <typename Field, typename LowerF, typename UpperF, typename InvDiagF, typename DiagF>
428 : struct ssor_preconditioner : public preconditioner<Field> {
429 : constexpr static unsigned Dim = Field::dim;
430 : using mesh_type = typename Field::Mesh_t;
431 : using layout_type = typename Field::Layout_t;
432 :
433 0 : ssor_preconditioner(LowerF&& lower, UpperF&& upper, InvDiagF&& inverse_diagonal,
434 : DiagF&& diagonal, unsigned innerloops, unsigned outerloops,
435 : double omega)
436 : : preconditioner<Field>("ssor")
437 0 : , innerloops_m(innerloops)
438 0 : , outerloops_m(outerloops)
439 0 : , omega_m(omega) {
440 0 : lower_m = std::move(lower);
441 0 : upper_m = std::move(upper);
442 0 : inverse_diagonal_m = std::move(inverse_diagonal);
443 0 : diagonal_m = std::move(diagonal);
444 0 : }
445 :
446 0 : Field operator()(Field& b) override {
447 0 : static IpplTimings::TimerRef initTimer = IpplTimings::getTimer("SSOR Init");
448 0 : IpplTimings::startTimer(initTimer);
449 :
450 : double D;
451 :
452 0 : layout_type& layout = b.getLayout();
453 0 : mesh_type& mesh = b.get_mesh();
454 :
455 0 : Field x(mesh, layout);
456 :
457 0 : x = 0; // Initial guess
458 :
459 0 : IpplTimings::stopTimer(initTimer);
460 :
461 0 : static IpplTimings::TimerRef loopTimer = IpplTimings::getTimer("SSOR loop");
462 0 : IpplTimings::startTimer(loopTimer);
463 :
464 : // The inverse diagonal is applied to the
465 : // vector itself to return the result usually.
466 : // However, the operator for FEM already
467 : // returns the result of inv_diag * itself
468 : // due to the matrix-free evaluation.
469 : // Therefore, we need this if to differentiate
470 : // the two cases.
471 0 : for (unsigned int k = 0; k < outerloops_m; ++k) {
472 : if constexpr (std::is_same_v<DiagF, double>) {
473 : UL_m = upper_m(x);
474 : D = diagonal_m(x);
475 : r_m = omega_m * (b - UL_m) + (1.0 - omega_m) * D * x;
476 :
477 : for (unsigned int j = 0; j < innerloops_m; ++j) {
478 : UL_m = lower_m(x);
479 : x = r_m - omega_m * UL_m;
480 : x = inverse_diagonal_m(x) * x;
481 : }
482 : UL_m = lower_m(x);
483 : D = diagonal_m(x);
484 : r_m = omega_m * (b - UL_m) + (1.0 - omega_m) * D * x;
485 : for (unsigned int j = 0; j < innerloops_m; ++j) {
486 : UL_m = upper_m(x);
487 : x = r_m - omega_m * UL_m;
488 : x = inverse_diagonal_m(x) * x;
489 : }
490 : } else {
491 0 : UL_m = upper_m(x);
492 0 : r_m = omega_m * (b - UL_m) + (1.0 - omega_m) * diagonal_m(x);
493 :
494 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
495 0 : UL_m = lower_m(x);
496 0 : x = r_m - omega_m * UL_m;
497 0 : x = inverse_diagonal_m(x);
498 : }
499 0 : UL_m = lower_m(x);
500 0 : r_m = omega_m * (b - UL_m) + (1.0 - omega_m) * diagonal_m(x);
501 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
502 0 : UL_m = upper_m(x);
503 0 : x = r_m - omega_m * UL_m;
504 0 : x = inverse_diagonal_m(x);
505 : }
506 : }
507 : }
508 0 : IpplTimings::stopTimer(loopTimer);
509 0 : return x;
510 0 : }
511 :
512 0 : void init_fields(Field& b) override {
513 0 : layout_type& layout = b.getLayout();
514 0 : mesh_type& mesh = b.get_mesh();
515 :
516 0 : UL_m = Field(mesh, layout);
517 0 : r_m = Field(mesh, layout);
518 0 : }
519 :
520 : protected:
521 : LowerF lower_m;
522 : UpperF upper_m;
523 : InvDiagF inverse_diagonal_m;
524 : DiagF diagonal_m;
525 : unsigned innerloops_m;
526 : unsigned outerloops_m;
527 : double omega_m;
528 : Field UL_m;
529 : Field r_m;
530 : };
531 :
532 : /*!
533 : * Computes the largest Eigenvalue of the Functor f
534 : * @param f Functor
535 : * @param x_0 initial guess
536 : * @param max_iter maximum number of iterations
537 : * @param tol tolerance
538 : */
539 : template <typename Field, typename Functor>
540 : double powermethod(Functor&& f, Field& x_0, unsigned int max_iter = 5000, double tol = 1e-3) {
541 : unsigned int i = 0;
542 : using mesh_type = typename Field::Mesh_t;
543 : using layout_type = typename Field::Layout_t;
544 : mesh_type& mesh = x_0.get_mesh();
545 : layout_type& layout = x_0.getLayout();
546 : Field x_new(mesh, layout);
547 : Field x_diff(mesh, layout);
548 : double error = 1.0;
549 : double lambda = 1.0;
550 : while (error > tol && i < max_iter) {
551 : x_new = f(x_0);
552 : lambda = norm(x_new);
553 : x_diff = x_new - lambda * x_0;
554 : error = norm(x_diff);
555 : x_new = x_new / lambda;
556 : x_0 = x_new.deepCopy();
557 : ++i;
558 : }
559 : if (i == max_iter) {
560 : std::cerr << "Powermethod did not converge, lambda_max : " << lambda
561 : << ", error : " << error << std::endl;
562 : }
563 : return lambda;
564 : }
565 :
566 : /*!
567 : * Computes the smallest Eigenvalue of the Functor f (f must be symmetric positive definite)
568 : * @param f Functor
569 : * @param x_0 initial guess
570 : * @param lambda_max largest Eigenvalue
571 : * @param max_iter maximum number of iterations
572 : * @param tol tolerance
573 : */
574 : template <typename Field, typename Functor>
575 : double adapted_powermethod(Functor&& f, Field& x_0, double lambda_max,
576 : unsigned int max_iter = 5000, double tol = 1e-3) {
577 : unsigned int i = 0;
578 : using mesh_type = typename Field::Mesh_t;
579 : using layout_type = typename Field::Layout_t;
580 : mesh_type& mesh = x_0.get_mesh();
581 : layout_type& layout = x_0.getLayout();
582 : Field x_new(mesh, layout);
583 : Field x_diff(mesh, layout);
584 : double error = 1.0;
585 : double lambda = 1.0;
586 : while (error > tol && i < max_iter) {
587 : x_new = f(x_0);
588 : x_new = x_new - lambda_max * x_0;
589 : lambda = -norm(x_new); // We know that lambda < 0;
590 : x_diff = x_new - lambda * x_0;
591 : error = norm(x_diff);
592 : x_new = x_new / -lambda;
593 : x_0 = x_new.deepCopy();
594 : ++i;
595 : }
596 : lambda = lambda + lambda_max;
597 : if (i == max_iter) {
598 : std::cerr << "Powermethod did not converge, lambda_min : " << lambda
599 : << ", error : " << error << std::endl;
600 : }
601 : return lambda;
602 : }
603 :
604 : /*
605 : // Use the powermethod to compute the eigenvalues if no analytical solution is known
606 : beta = powermethod(std::move(op_m), x_0);
607 : // Trick for computing the smallest Eigenvalue of SPD Matrix
608 : alpha = adapted_powermethod(std::move(op_m), x_0, beta);
609 : */
610 :
611 : } // namespace ippl
612 :
613 : #endif // IPPL_PRECONDITIONER_H
|
