Branch data 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
|
