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