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 : if constexpr (std::is_same_v<InvDiagF, double>) {
308 : g = inverse_diagonal_m(g) * g;
309 : } else {
310 0 : g = inverse_diagonal_m(g);
311 : }
312 : }
313 0 : return g;
314 0 : }
315 :
316 0 : void init_fields(Field& b) override {
317 0 : layout_type& layout = b.getLayout();
318 0 : mesh_type& mesh = b.get_mesh();
319 :
320 0 : ULg_m = Field(mesh, layout);
321 0 : }
322 :
323 : protected:
324 : UpperAndLowerF upper_and_lower_m;
325 : InvDiagF inverse_diagonal_m;
326 : unsigned innerloops_m;
327 : Field ULg_m;
328 : };
329 :
330 : /*!
331 : * 2-step Gauss-Seidel preconditioner
332 : */
333 : template <typename Field, typename LowerF, typename UpperF, typename InvDiagF>
334 : struct gs_preconditioner : public preconditioner<Field> {
335 : constexpr static unsigned Dim = Field::dim;
336 : using mesh_type = typename Field::Mesh_t;
337 : using layout_type = typename Field::Layout_t;
338 :
339 0 : gs_preconditioner(LowerF&& lower, UpperF&& upper, InvDiagF&& inverse_diagonal,
340 : unsigned innerloops, unsigned outerloops)
341 : : preconditioner<Field>("Gauss-Seidel")
342 0 : , innerloops_m(innerloops)
343 0 : , outerloops_m(outerloops) {
344 0 : lower_m = std::move(lower);
345 0 : upper_m = std::move(upper);
346 0 : inverse_diagonal_m = std::move(inverse_diagonal);
347 0 : }
348 :
349 0 : Field operator()(Field& b) override {
350 0 : layout_type& layout = b.getLayout();
351 0 : mesh_type& mesh = b.get_mesh();
352 :
353 0 : Field x(mesh, layout);
354 :
355 0 : x = 0; // Initial guess
356 :
357 0 : for (unsigned int k = 0; k < outerloops_m; ++k) {
358 0 : UL_m = upper_m(x);
359 0 : r_m = b - UL_m;
360 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
361 0 : UL_m = lower_m(x);
362 0 : x = r_m - UL_m;
363 : if constexpr (std::is_same_v<InvDiagF, double>) {
364 : x = inverse_diagonal_m(x) * x;
365 : } else {
366 0 : x = inverse_diagonal_m(x);
367 : }
368 : }
369 0 : UL_m = lower_m(x);
370 0 : r_m = b - UL_m;
371 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
372 0 : UL_m = upper_m(x);
373 0 : x = r_m - UL_m;
374 : if constexpr (std::is_same_v<InvDiagF, double>) {
375 : x = inverse_diagonal_m(x) * x;
376 : } else {
377 0 : x = inverse_diagonal_m(x);
378 : }
379 : }
380 : }
381 0 : return x;
382 0 : }
383 :
384 0 : void init_fields(Field& b) override {
385 0 : layout_type& layout = b.getLayout();
386 0 : mesh_type& mesh = b.get_mesh();
387 :
388 0 : UL_m = Field(mesh, layout);
389 0 : r_m = Field(mesh, layout);
390 0 : }
391 :
392 : protected:
393 : LowerF lower_m;
394 : UpperF upper_m;
395 : InvDiagF inverse_diagonal_m;
396 : unsigned innerloops_m;
397 : unsigned outerloops_m;
398 : Field UL_m;
399 : Field r_m;
400 : };
401 :
402 : /*!
403 : * Symmetric successive over-relaxation
404 : */
405 : template <typename Field, typename LowerF, typename UpperF, typename InvDiagF, typename DiagF>
406 : struct ssor_preconditioner : public preconditioner<Field> {
407 : constexpr static unsigned Dim = Field::dim;
408 : using mesh_type = typename Field::Mesh_t;
409 : using layout_type = typename Field::Layout_t;
410 :
411 0 : ssor_preconditioner(LowerF&& lower, UpperF&& upper, InvDiagF&& inverse_diagonal,
412 : DiagF&& diagonal, unsigned innerloops, unsigned outerloops,
413 : double omega)
414 : : preconditioner<Field>("ssor")
415 0 : , innerloops_m(innerloops)
416 0 : , outerloops_m(outerloops)
417 0 : , omega_m(omega) {
418 0 : lower_m = std::move(lower);
419 0 : upper_m = std::move(upper);
420 0 : inverse_diagonal_m = std::move(inverse_diagonal);
421 0 : diagonal_m = std::move(diagonal);
422 0 : }
423 :
424 0 : Field operator()(Field& b) override {
425 0 : static IpplTimings::TimerRef initTimer = IpplTimings::getTimer("SSOR Init");
426 0 : IpplTimings::startTimer(initTimer);
427 :
428 : double D;
429 :
430 0 : layout_type& layout = b.getLayout();
431 0 : mesh_type& mesh = b.get_mesh();
432 :
433 0 : Field x(mesh, layout);
434 :
435 0 : x = 0; // Initial guess
436 :
437 0 : IpplTimings::stopTimer(initTimer);
438 :
439 0 : static IpplTimings::TimerRef loopTimer = IpplTimings::getTimer("SSOR loop");
440 0 : IpplTimings::startTimer(loopTimer);
441 :
442 0 : for (unsigned int k = 0; k < outerloops_m; ++k) {
443 : if constexpr (std::is_same_v<DiagF, double>) {
444 : UL_m = upper_m(x);
445 : D = diagonal_m(x);
446 : r_m = omega_m * (b - UL_m) + (1.0 - omega_m) * D * x;
447 :
448 : for (unsigned int j = 0; j < innerloops_m; ++j) {
449 : UL_m = lower_m(x);
450 : x = r_m - omega_m * UL_m;
451 : x = inverse_diagonal_m(x) * x;
452 : }
453 : UL_m = lower_m(x);
454 : D = diagonal_m(x);
455 : r_m = omega_m * (b - UL_m) + (1.0 - omega_m) * D * x;
456 : for (unsigned int j = 0; j < innerloops_m; ++j) {
457 : UL_m = upper_m(x);
458 : x = r_m - omega_m * UL_m;
459 : x = inverse_diagonal_m(x) * x;
460 : }
461 : } else {
462 0 : UL_m = upper_m(x);
463 0 : r_m = omega_m * (b - UL_m) + (1.0 - omega_m) * diagonal_m(x);
464 :
465 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
466 0 : UL_m = lower_m(x);
467 0 : x = r_m - omega_m * UL_m;
468 0 : x = inverse_diagonal_m(x);
469 : }
470 0 : UL_m = lower_m(x);
471 0 : r_m = omega_m * (b - UL_m) + (1.0 - omega_m) * diagonal_m(x);
472 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
473 0 : UL_m = upper_m(x);
474 0 : x = r_m - omega_m * UL_m;
475 0 : x = inverse_diagonal_m(x);
476 : }
477 : }
478 : }
479 0 : IpplTimings::stopTimer(loopTimer);
480 0 : return x;
481 0 : }
482 :
483 0 : void init_fields(Field& b) override {
484 0 : layout_type& layout = b.getLayout();
485 0 : mesh_type& mesh = b.get_mesh();
486 :
487 0 : UL_m = Field(mesh, layout);
488 0 : r_m = Field(mesh, layout);
489 0 : }
490 :
491 : protected:
492 : LowerF lower_m;
493 : UpperF upper_m;
494 : InvDiagF inverse_diagonal_m;
495 : DiagF diagonal_m;
496 : unsigned innerloops_m;
497 : unsigned outerloops_m;
498 : double omega_m;
499 : Field UL_m;
500 : Field r_m;
501 : };
502 :
503 : /*!
504 : * Computes the largest Eigenvalue of the Functor f
505 : * @param f Functor
506 : * @param x_0 initial guess
507 : * @param max_iter maximum number of iterations
508 : * @param tol tolerance
509 : */
510 : template <typename Field, typename Functor>
511 : double powermethod(Functor&& f, Field& x_0, unsigned int max_iter = 5000, double tol = 1e-3) {
512 : unsigned int i = 0;
513 : using mesh_type = typename Field::Mesh_t;
514 : using layout_type = typename Field::Layout_t;
515 : mesh_type& mesh = x_0.get_mesh();
516 : layout_type& layout = x_0.getLayout();
517 : Field x_new(mesh, layout);
518 : Field x_diff(mesh, layout);
519 : double error = 1.0;
520 : double lambda = 1.0;
521 : while (error > tol && i < max_iter) {
522 : x_new = f(x_0);
523 : lambda = norm(x_new);
524 : x_diff = x_new - lambda * x_0;
525 : error = norm(x_diff);
526 : x_new = x_new / lambda;
527 : x_0 = x_new.deepCopy();
528 : ++i;
529 : }
530 : if (i == max_iter) {
531 : std::cerr << "Powermethod did not converge, lambda_max : " << lambda
532 : << ", error : " << error << std::endl;
533 : }
534 : return lambda;
535 : }
536 :
537 : /*!
538 : * Computes the smallest Eigenvalue of the Functor f (f must be symmetric positive definite)
539 : * @param f Functor
540 : * @param x_0 initial guess
541 : * @param lambda_max largest Eigenvalue
542 : * @param max_iter maximum number of iterations
543 : * @param tol tolerance
544 : */
545 : template <typename Field, typename Functor>
546 : double adapted_powermethod(Functor&& f, Field& x_0, double lambda_max,
547 : unsigned int max_iter = 5000, double tol = 1e-3) {
548 : unsigned int i = 0;
549 : using mesh_type = typename Field::Mesh_t;
550 : using layout_type = typename Field::Layout_t;
551 : mesh_type& mesh = x_0.get_mesh();
552 : layout_type& layout = x_0.getLayout();
553 : Field x_new(mesh, layout);
554 : Field x_diff(mesh, layout);
555 : double error = 1.0;
556 : double lambda = 1.0;
557 : while (error > tol && i < max_iter) {
558 : x_new = f(x_0);
559 : x_new = x_new - lambda_max * x_0;
560 : lambda = -norm(x_new); // We know that lambda < 0;
561 : x_diff = x_new - lambda * x_0;
562 : error = norm(x_diff);
563 : x_new = x_new / -lambda;
564 : x_0 = x_new.deepCopy();
565 : ++i;
566 : }
567 : lambda = lambda + lambda_max;
568 : if (i == max_iter) {
569 : std::cerr << "Powermethod did not converge, lambda_min : " << lambda
570 : << ", error : " << error << std::endl;
571 : }
572 : return lambda;
573 : }
574 :
575 : /*
576 : // Use the powermethod to compute the eigenvalues if no analytical solution is known
577 : beta = powermethod(std::move(op_m), x_0);
578 : // Trick for computing the smallest Eigenvalue of SPD Matrix
579 : alpha = adapted_powermethod(std::move(op_m), x_0, beta);
580 : */
581 :
582 : } // namespace ippl
583 :
584 : #endif // IPPL_PRECONDITIONER_H
|
