diff --git a/PCG_8h_source.html b/PCG_8h_source.html index 378a20d97..1691ea051 100644 --- a/PCG_8h_source.html +++ b/PCG_8h_source.html @@ -296,90 +296,94 @@
223 void operator()(lhs_type& lhs, rhs_type& rhs, const ParameterList& params) override {
224 constexpr unsigned Dim = lhs_type::dim;
225
-
226 typename lhs_type::Mesh_t mesh = lhs.get_mesh();
-
227 typename lhs_type::Layout_t layout = lhs.getLayout();
-
228
-
229 this->iterations_m = 0;
-
230 const int maxIterations = params.get<int>("max_iterations");
-
231
-
232 // Variable names mostly based on description in
-
233 // https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf
-
234 lhs_type r(mesh, layout);
-
235 lhs_type d(mesh, layout);
-
236 lhs_type s(mesh, layout);
-
237 lhs_type q(mesh, layout);
-
238
-
239 using bc_type = BConds<lhs_type, Dim>;
-
240 bc_type lhsBCs = lhs.getFieldBC();
-
241 bc_type bc;
+
226 if (preconditioner_m == nullptr) {
+
227 throw IpplException("PCG::operator()", "Preconditioner has not been set for PCG solver");
+
228 }
+
229
+
230 typename lhs_type::Mesh_t mesh = lhs.get_mesh();
+
231 typename lhs_type::Layout_t layout = lhs.getLayout();
+
232
+
233 this->iterations_m = 0;
+
234 const int maxIterations = params.get<int>("max_iterations");
+
235
+
236 // Variable names mostly based on description in
+
237 // https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf
+
238 lhs_type r(mesh, layout);
+
239 lhs_type d(mesh, layout);
+
240 lhs_type s(mesh, layout);
+
241 lhs_type q(mesh, layout);
242
-
243 bool allFacesPeriodic = true;
-
244 for (unsigned int i = 0; i < 2 * Dim; ++i) {
-
245 FieldBC bcType = lhsBCs[i]->getBCType();
-
246 if (bcType == PERIODIC_FACE) {
-
247 // If the LHS has periodic BCs, so does the residue
-
248 bc[i] = std::make_shared<PeriodicFace<lhs_type>>(i);
-
249 } else if (bcType & CONSTANT_FACE) {
-
250 // If the LHS has constant BCs, the residue is zero on the BCs
-
251 // Bitwise AND with CONSTANT_FACE will succeed for ZeroFace or ConstantFace
-
252 bc[i] = std::make_shared<ZeroFace<lhs_type>>(i);
-
253 allFacesPeriodic = false;
-
254 } else {
-
255 throw IpplException("PCG::operator()",
-
256 "Only periodic or constant BCs for LHS supported.");
-
257 return;
-
258 }
-
259 }
-
260
-
261 r = rhs - this->op_m(lhs);
-
262 d = preconditioner_m->operator()(r);
-
263 d.setFieldBC(bc);
+
243 using bc_type = BConds<lhs_type, Dim>;
+
244 bc_type lhsBCs = lhs.getFieldBC();
+
245 bc_type bc;
+
246
+
247 bool allFacesPeriodic = true;
+
248 for (unsigned int i = 0; i < 2 * Dim; ++i) {
+
249 FieldBC bcType = lhsBCs[i]->getBCType();
+
250 if (bcType == PERIODIC_FACE) {
+
251 // If the LHS has periodic BCs, so does the residue
+
252 bc[i] = std::make_shared<PeriodicFace<lhs_type>>(i);
+
253 } else if (bcType & CONSTANT_FACE) {
+
254 // If the LHS has constant BCs, the residue is zero on the BCs
+
255 // Bitwise AND with CONSTANT_FACE will succeed for ZeroFace or ConstantFace
+
256 bc[i] = std::make_shared<ZeroFace<lhs_type>>(i);
+
257 allFacesPeriodic = false;
+
258 } else {
+
259 throw IpplException("PCG::operator()",
+
260 "Only periodic or constant BCs for LHS supported.");
+
261 return;
+
262 }
+
263 }
264
-
265 T delta1 = innerProduct(r, d);
-
266 T delta0 = delta1;
-
267 this->residueNorm = std::sqrt(std::abs(delta1));
-
268 const T tolerance = params.get<T>("tolerance") * delta1;
-
269
-
270 while (this->iterations_m<maxIterations&& this->residueNorm> tolerance) {
-
271 q = this->op_m(d);
-
272 T alpha = delta1 / innerProduct(d, q);
-
273 lhs = lhs + alpha * d;
-
274
-
275 // The exact residue is given by
-
276 // r = rhs - BaseCG::op_m(lhs);
-
277 // This correction is generally not used in practice because
-
278 // applying the Laplacian is computationally expensive and
-
279 // the correction does not have a significant effect on accuracy;
-
280 // in some implementations, the correction may be applied every few
-
281 // iterations to offset accumulated floating point errors
-
282 r = r - alpha * q;
-
283 s = preconditioner_m->operator()(r);
-
284
-
285 delta0 = delta1;
-
286 delta1 = innerProduct(r, s);
-
287
-
288 T beta = delta1 / delta0;
-
289 this->residueNorm = std::sqrt(std::abs(delta1));
-
290
-
291 d = s + beta * d;
-
292 ++this->iterations_m;
-
293 }
+
265 r = rhs - this->op_m(lhs);
+
266 d = preconditioner_m->operator()(r);
+
267 d.setFieldBC(bc);
+
268
+
269 T delta1 = innerProduct(r, d);
+
270 T delta0 = delta1;
+
271 this->residueNorm = std::sqrt(std::abs(delta1));
+
272 const T tolerance = params.get<T>("tolerance") * delta1;
+
273
+
274 while (this->iterations_m<maxIterations&& this->residueNorm> tolerance) {
+
275 q = this->op_m(d);
+
276 T alpha = delta1 / innerProduct(d, q);
+
277 lhs = lhs + alpha * d;
+
278
+
279 // The exact residue is given by
+
280 // r = rhs - BaseCG::op_m(lhs);
+
281 // This correction is generally not used in practice because
+
282 // applying the Laplacian is computationally expensive and
+
283 // the correction does not have a significant effect on accuracy;
+
284 // in some implementations, the correction may be applied every few
+
285 // iterations to offset accumulated floating point errors
+
286 r = r - alpha * q;
+
287 s = preconditioner_m->operator()(r);
+
288
+
289 delta0 = delta1;
+
290 delta1 = innerProduct(r, s);
+
291
+
292 T beta = delta1 / delta0;
+
293 this->residueNorm = std::sqrt(std::abs(delta1));
294
-
295 if (allFacesPeriodic) {
-
296 T avg = lhs.getVolumeAverage();
-
297 lhs = lhs - avg;
-
298 }
-
299 }
-
300
-
301 protected:
-
302 std::unique_ptr<preconditioner<FieldLHS>> preconditioner_m;
-
303 };
+
295 d = s + beta * d;
+
296 ++this->iterations_m;
+
297 }
+
298
+
299 if (allFacesPeriodic) {
+
300 T avg = lhs.getVolumeAverage();
+
301 lhs = lhs - avg;
+
302 }
+
303 }
304
-
305}; // namespace ippl
-
306
-
307#endif
+
305 protected:
+
306 std::unique_ptr<preconditioner<FieldLHS>> preconditioner_m;
+
307 };
308
-
309
+
309}; // namespace ippl
+
310
+
311#endif
+
312
+
313
Dim
constexpr unsigned Dim
Definition: BumponTailInstability.cpp:22
Preconditioner.h
SolverAlgorithm.h
@@ -405,7 +409,7 @@
ippl::PCG
Definition: PCG.h:148
ippl::PCG::UpperLowerF
std::function< UpperLowerRet(lhs_type)> UpperLowerF
Definition: PCG.h:157
ippl::PCG::setPreconditioner
void setPreconditioner(OperatorF &&op, LowerF &&lower, UpperF &&upper, UpperLowerF &&upper_and_lower, InverseDiagF &&inverse_diagonal, double alpha, double beta, std::string preconditioner_type="", int level=5, int degree=31, int richardson_iterations=1, int inner=5, int outer=1) override
Definition: PCG.h:169
-
ippl::PCG::preconditioner_m
std::unique_ptr< preconditioner< FieldLHS > > preconditioner_m
Definition: PCG.h:302
+
ippl::PCG::preconditioner_m
std::unique_ptr< preconditioner< FieldLHS > > preconditioner_m
Definition: PCG.h:306
ippl::PCG::operator()
void operator()(lhs_type &lhs, rhs_type &rhs, const ParameterList &params) override
Definition: PCG.h:223
ippl::PCG::PCG
PCG()
Definition: PCG.h:160
ippl::PCG::rhs_type
FieldRHS rhs_type
Definition: SolverAlgorithm.h:19