Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
petsc_precondition.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 
18 #ifdef DEAL_II_WITH_PETSC
19 
20 # include <deal.II/base/utilities.h>
21 
22 # include <deal.II/lac/exceptions.h>
27 
28 # include <petscconf.h>
29 
30 # include <cmath>
31 
33 
34 namespace PETScWrappers
35 {
37  : pc(nullptr)
38  , matrix(nullptr)
39  {}
40 
42  {
43  try
44  {
45  clear();
46  }
47  catch (...)
48  {}
49  }
50 
51  void
53  {
54  matrix = nullptr;
55 
56  if (pc != nullptr)
57  {
58  PetscErrorCode ierr = PCDestroy(&pc);
59  pc = nullptr;
60  AssertThrow(ierr == 0, ExcPETScError(ierr));
61  }
62  }
63 
64 
65  void
67  {
69 
70  const PetscErrorCode ierr = PCApply(pc, src, dst);
71  AssertThrow(ierr == 0, ExcPETScError(ierr));
72  }
73 
74 
75  void
77  {
79 
80  const PetscErrorCode ierr = PCApplyTranspose(pc, src, dst);
81  AssertThrow(ierr == 0, ExcPETScError(ierr));
82  }
83 
84 
85  void
87  {
88  // only allow the creation of the
89  // preconditioner once
91 
92  MPI_Comm comm;
93  // this ugly cast is necessary because the
94  // type Mat and PETScObject are
95  // unrelated.
96  PetscErrorCode ierr =
97  PetscObjectGetComm(reinterpret_cast<PetscObject>(matrix), &comm);
98  AssertThrow(ierr == 0, ExcPETScError(ierr));
99 
100  ierr = PCCreate(comm, &pc);
101  AssertThrow(ierr == 0, ExcPETScError(ierr));
102 
103 # if DEAL_II_PETSC_VERSION_LT(3, 5, 0)
104  ierr = PCSetOperators(pc, matrix, matrix, SAME_PRECONDITIONER);
105 # else
106  ierr = PCSetOperators(pc, matrix, matrix);
107 # endif
108  AssertThrow(ierr == 0, ExcPETScError(ierr));
109  }
110 
111 
112  const PC &
114  {
115  return pc;
116  }
117 
118 
119  PreconditionerBase::operator Mat() const
120  {
121  return matrix;
122  }
123 
124 
125  /* ----------------- PreconditionJacobi -------------------- */
127  const AdditionalData &additional_data_)
128  {
129  additional_data = additional_data_;
130 
131  PetscErrorCode ierr = PCCreate(comm, &pc);
132  AssertThrow(ierr == 0, ExcPETScError(ierr));
133 
134  initialize();
135  }
136 
137 
138 
140  const AdditionalData &additional_data)
141  {
142  initialize(matrix, additional_data);
143  }
144 
145  void
147  {
149 
150  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCJACOBI));
151  AssertThrow(ierr == 0, ExcPETScError(ierr));
152 
153  ierr = PCSetFromOptions(pc);
154  AssertThrow(ierr == 0, ExcPETScError(ierr));
155  }
156 
157  void
159  const AdditionalData &additional_data_)
160  {
161  clear();
162 
163  matrix = static_cast<Mat>(matrix_);
164  additional_data = additional_data_;
165 
166  create_pc();
167  initialize();
168 
169  PetscErrorCode ierr = PCSetUp(pc);
170  AssertThrow(ierr == 0, ExcPETScError(ierr));
171  }
172 
173 
174  /* ----------------- PreconditionBlockJacobi -------------------- */
176  const MPI_Comm comm,
177  const AdditionalData &additional_data_)
178  {
179  additional_data = additional_data_;
180 
181  PetscErrorCode ierr = PCCreate(comm, &pc);
182  AssertThrow(ierr == 0, ExcPETScError(ierr));
183 
184  initialize();
185  }
186 
187 
188 
190  const MatrixBase & matrix,
191  const AdditionalData &additional_data)
192  {
193  initialize(matrix, additional_data);
194  }
195 
196  void
198  {
199  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCBJACOBI));
200  AssertThrow(ierr == 0, ExcPETScError(ierr));
201 
202  ierr = PCSetFromOptions(pc);
203  AssertThrow(ierr == 0, ExcPETScError(ierr));
204  }
205 
206 
207  void
209  const AdditionalData &additional_data_)
210  {
211  clear();
212 
213  matrix = static_cast<Mat>(matrix_);
214  additional_data = additional_data_;
215 
216  create_pc();
217  initialize();
218 
219  PetscErrorCode ierr = PCSetUp(pc);
220  AssertThrow(ierr == 0, ExcPETScError(ierr));
221  }
222 
223 
224  /* ----------------- PreconditionSOR -------------------- */
225 
227  : omega(omega)
228  {}
229 
230 
231 
234  {
235  initialize(matrix, additional_data);
236  }
237 
238 
239  void
241  const AdditionalData &additional_data_)
242  {
243  clear();
244 
245  matrix = static_cast<Mat>(matrix_);
246  additional_data = additional_data_;
247 
248  create_pc();
249 
250  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
251  AssertThrow(ierr == 0, ExcPETScError(ierr));
252 
253  // then set flags as given
254  ierr = PCSORSetOmega(pc, additional_data.omega);
255  AssertThrow(ierr == 0, ExcPETScError(ierr));
256 
257  ierr = PCSetFromOptions(pc);
258  AssertThrow(ierr == 0, ExcPETScError(ierr));
259 
260  ierr = PCSetUp(pc);
261  AssertThrow(ierr == 0, ExcPETScError(ierr));
262  }
263 
264 
265  /* ----------------- PreconditionSSOR -------------------- */
266 
268  : omega(omega)
269  {}
270 
271 
272 
275  {
276  initialize(matrix, additional_data);
277  }
278 
279 
280  void
282  const AdditionalData &additional_data_)
283  {
284  clear();
285 
286  matrix = static_cast<Mat>(matrix_);
287  additional_data = additional_data_;
288 
289  create_pc();
290 
291  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
292  AssertThrow(ierr == 0, ExcPETScError(ierr));
293 
294  // then set flags as given
295  ierr = PCSORSetOmega(pc, additional_data.omega);
296  AssertThrow(ierr == 0, ExcPETScError(ierr));
297 
298  // convert SOR to SSOR
299  ierr = PCSORSetSymmetric(pc, SOR_SYMMETRIC_SWEEP);
300  AssertThrow(ierr == 0, ExcPETScError(ierr));
301 
302  ierr = PCSetFromOptions(pc);
303  AssertThrow(ierr == 0, ExcPETScError(ierr));
304 
305  ierr = PCSetUp(pc);
306  AssertThrow(ierr == 0, ExcPETScError(ierr));
307  }
308 
309 
310  /* ----------------- PreconditionEisenstat -------------------- */
311 
313  : omega(omega)
314  {}
315 
316 
317 
319  const MatrixBase & matrix,
321  {
322  initialize(matrix, additional_data);
323  }
324 
325 
326  void
328  const AdditionalData &additional_data_)
329  {
330  clear();
331 
332  matrix = static_cast<Mat>(matrix_);
333  additional_data = additional_data_;
334 
335  create_pc();
336 
337  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCEISENSTAT));
338  AssertThrow(ierr == 0, ExcPETScError(ierr));
339 
340  // then set flags as given
341  ierr = PCEisenstatSetOmega(pc, additional_data.omega);
342  AssertThrow(ierr == 0, ExcPETScError(ierr));
343 
344  ierr = PCSetFromOptions(pc);
345  AssertThrow(ierr == 0, ExcPETScError(ierr));
346 
347  ierr = PCSetUp(pc);
348  AssertThrow(ierr == 0, ExcPETScError(ierr));
349  }
350 
351 
352  /* ----------------- PreconditionICC -------------------- */
353 
354 
356  : levels(levels)
357  {}
358 
359 
360 
363  {
364  initialize(matrix, additional_data);
365  }
366 
367 
368  void
370  const AdditionalData &additional_data_)
371  {
372  clear();
373 
374  matrix = static_cast<Mat>(matrix_);
375  additional_data = additional_data_;
376 
377  create_pc();
378 
379  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCICC));
380  AssertThrow(ierr == 0, ExcPETScError(ierr));
381 
382  // then set flags
383  ierr = PCFactorSetLevels(pc, additional_data.levels);
384  AssertThrow(ierr == 0, ExcPETScError(ierr));
385 
386  ierr = PCSetFromOptions(pc);
387  AssertThrow(ierr == 0, ExcPETScError(ierr));
388 
389  ierr = PCSetUp(pc);
390  AssertThrow(ierr == 0, ExcPETScError(ierr));
391  }
392 
393 
394  /* ----------------- PreconditionILU -------------------- */
395 
397  : levels(levels)
398  {}
399 
400 
401 
404  {
405  initialize(matrix, additional_data);
406  }
407 
408 
409  void
411  const AdditionalData &additional_data_)
412  {
413  clear();
414 
415  matrix = static_cast<Mat>(matrix_);
416  additional_data = additional_data_;
417 
418  create_pc();
419 
420  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCILU));
421  AssertThrow(ierr == 0, ExcPETScError(ierr));
422 
423  // then set flags
424  ierr = PCFactorSetLevels(pc, additional_data.levels);
425  AssertThrow(ierr == 0, ExcPETScError(ierr));
426 
427  ierr = PCSetFromOptions(pc);
428  AssertThrow(ierr == 0, ExcPETScError(ierr));
429 
430  ierr = PCSetUp(pc);
431  AssertThrow(ierr == 0, ExcPETScError(ierr));
432  }
433 
434 
435  /* ----------------- PreconditionBoomerAMG -------------------- */
436 
438  const bool symmetric_operator,
439  const double strong_threshold,
440  const double max_row_sum,
441  const unsigned int aggressive_coarsening_num_levels,
442  const bool output_details)
443  : symmetric_operator(symmetric_operator)
444  , strong_threshold(strong_threshold)
445  , max_row_sum(max_row_sum)
446  , aggressive_coarsening_num_levels(aggressive_coarsening_num_levels)
447  , output_details(output_details)
448  {}
449 
450 
451 
453  const MPI_Comm comm,
454  const AdditionalData &additional_data_)
455  {
456  additional_data = additional_data_;
457 
458  PetscErrorCode ierr = PCCreate(comm, &pc);
459  AssertThrow(ierr == 0, ExcPETScError(ierr));
460 
461 # ifdef DEAL_II_PETSC_WITH_HYPRE
462  initialize();
463 # else // DEAL_II_PETSC_WITH_HYPRE
464  (void)pc;
465  Assert(false,
466  ExcMessage("Your PETSc installation does not include a copy of "
467  "the hypre package necessary for this preconditioner."));
468 # endif
469  }
470 
471 
473  const MatrixBase & matrix,
475  {
476  initialize(matrix, additional_data);
477  }
478 
479  void
481  {
482 # ifdef DEAL_II_PETSC_WITH_HYPRE
483  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
484  AssertThrow(ierr == 0, ExcPETScError(ierr));
485 
486  ierr = PCHYPRESetType(pc, "boomeramg");
487  AssertThrow(ierr == 0, ExcPETScError(ierr));
488 
490  {
491  set_option_value("-pc_hypre_boomeramg_print_statistics", "1");
492  }
493 
494  set_option_value("-pc_hypre_boomeramg_agg_nl",
497 
498  std::stringstream ssStream;
499  ssStream << additional_data.max_row_sum;
500  set_option_value("-pc_hypre_boomeramg_max_row_sum", ssStream.str());
501 
502  ssStream.str(""); // empty the stringstream
503  ssStream << additional_data.strong_threshold;
504  set_option_value("-pc_hypre_boomeramg_strong_threshold", ssStream.str());
505 
507  {
508  set_option_value("-pc_hypre_boomeramg_relax_type_up",
509  "symmetric-SOR/Jacobi");
510  set_option_value("-pc_hypre_boomeramg_relax_type_down",
511  "symmetric-SOR/Jacobi");
512  set_option_value("-pc_hypre_boomeramg_relax_type_coarse",
513  "Gaussian-elimination");
514  }
515  else
516  {
517  set_option_value("-pc_hypre_boomeramg_relax_type_up", "SOR/Jacobi");
518  set_option_value("-pc_hypre_boomeramg_relax_type_down", "SOR/Jacobi");
519  set_option_value("-pc_hypre_boomeramg_relax_type_coarse",
520  "Gaussian-elimination");
521  }
522 
523  ierr = PCSetFromOptions(pc);
524  AssertThrow(ierr == 0, ExcPETScError(ierr));
525 # else
526  Assert(false,
527  ExcMessage("Your PETSc installation does not include a copy of "
528  "the hypre package necessary for this preconditioner."));
529 # endif
530  }
531 
532  void
534  const AdditionalData &additional_data_)
535  {
536 # ifdef DEAL_II_PETSC_WITH_HYPRE
537  clear();
538 
539  matrix = static_cast<Mat>(matrix_);
540  additional_data = additional_data_;
541 
542  create_pc();
543  initialize();
544 
545  PetscErrorCode ierr = PCSetUp(pc);
546  AssertThrow(ierr == 0, ExcPETScError(ierr));
547 
548 # else // DEAL_II_PETSC_WITH_HYPRE
549  (void)matrix_;
550  (void)additional_data_;
551  Assert(false,
552  ExcMessage("Your PETSc installation does not include a copy of "
553  "the hypre package necessary for this preconditioner."));
554 # endif
555  }
556 
557 
558  /* ----------------- PreconditionParaSails -------------------- */
559 
561  const unsigned int symmetric,
562  const unsigned int n_levels,
563  const double threshold,
564  const double filter,
565  const bool output_details)
566  : symmetric(symmetric)
567  , n_levels(n_levels)
568  , threshold(threshold)
569  , filter(filter)
570  , output_details(output_details)
571  {}
572 
573 
574 
576  const MatrixBase & matrix,
578  {
579  initialize(matrix, additional_data);
580  }
581 
582 
583  void
585  const AdditionalData &additional_data_)
586  {
587  clear();
588 
589  matrix = static_cast<Mat>(matrix_);
590  additional_data = additional_data_;
591 
592 # ifdef DEAL_II_PETSC_WITH_HYPRE
593  create_pc();
594 
595  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
596  AssertThrow(ierr == 0, ExcPETScError(ierr));
597 
598  ierr = PCHYPRESetType(pc, "parasails");
599  AssertThrow(ierr == 0, ExcPETScError(ierr));
600 
602  {
603  set_option_value("-pc_hypre_parasails_logging", "1");
604  }
605 
608  ExcMessage(
609  "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
610 
611  std::stringstream ssStream;
612 
613  switch (additional_data.symmetric)
614  {
615  case 0:
616  {
617  ssStream << "nonsymmetric";
618  break;
619  }
620 
621  case 1:
622  {
623  ssStream << "SPD";
624  break;
625  }
626 
627  case 2:
628  {
629  ssStream << "nonsymmetric,SPD";
630  break;
631  }
632 
633  default:
634  Assert(
635  false,
636  ExcMessage(
637  "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
638  }
639 
640  set_option_value("-pc_hypre_parasails_sym", ssStream.str());
641 
642  set_option_value("-pc_hypre_parasails_nlevels",
644 
645  ssStream.str(""); // empty the stringstream
646  ssStream << additional_data.threshold;
647  set_option_value("-pc_hypre_parasails_thresh", ssStream.str());
648 
649  ssStream.str(""); // empty the stringstream
650  ssStream << additional_data.filter;
651  set_option_value("-pc_hypre_parasails_filter", ssStream.str());
652 
653  ierr = PCSetFromOptions(pc);
654  AssertThrow(ierr == 0, ExcPETScError(ierr));
655 
656  ierr = PCSetUp(pc);
657  AssertThrow(ierr == 0, ExcPETScError(ierr));
658 
659 # else // DEAL_II_PETSC_WITH_HYPRE
660  (void)pc;
661  Assert(false,
662  ExcMessage("Your PETSc installation does not include a copy of "
663  "the hypre package necessary for this preconditioner."));
664 # endif
665  }
666 
667 
668  /* ----------------- PreconditionNone ------------------------- */
669 
672  {
673  initialize(matrix, additional_data);
674  }
675 
676 
677  void
679  const AdditionalData &additional_data_)
680  {
681  clear();
682 
683  matrix = static_cast<Mat>(matrix_);
684  additional_data = additional_data_;
685 
686  create_pc();
687 
688  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCNONE));
689  AssertThrow(ierr == 0, ExcPETScError(ierr));
690 
691  ierr = PCSetFromOptions(pc);
692  AssertThrow(ierr == 0, ExcPETScError(ierr));
693 
694  ierr = PCSetUp(pc);
695  AssertThrow(ierr == 0, ExcPETScError(ierr));
696  }
697 
698 
699  /* ----------------- PreconditionLU -------------------- */
700 
702  const double zero_pivot,
703  const double damping)
704  : pivoting(pivoting)
705  , zero_pivot(zero_pivot)
706  , damping(damping)
707  {}
708 
709 
710 
713  {
714  initialize(matrix, additional_data);
715  }
716 
717 
718  void
720  const AdditionalData &additional_data_)
721  {
722  clear();
723 
724  matrix = static_cast<Mat>(matrix_);
725  additional_data = additional_data_;
726 
727  create_pc();
728 
729  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCLU));
730  AssertThrow(ierr == 0, ExcPETScError(ierr));
731 
732  // set flags as given
733  ierr = PCFactorSetColumnPivot(pc, additional_data.pivoting);
734  AssertThrow(ierr == 0, ExcPETScError(ierr));
735 
736  ierr = PCFactorSetZeroPivot(pc, additional_data.zero_pivot);
737  AssertThrow(ierr == 0, ExcPETScError(ierr));
738 
739  ierr = PCFactorSetShiftAmount(pc, additional_data.damping);
740  AssertThrow(ierr == 0, ExcPETScError(ierr));
741 
742  ierr = PCSetFromOptions(pc);
743  AssertThrow(ierr == 0, ExcPETScError(ierr));
744 
745  ierr = PCSetUp(pc);
746  AssertThrow(ierr == 0, ExcPETScError(ierr));
747  }
748 
749 } // namespace PETScWrappers
750 
752 
753 #endif // DEAL_II_WITH_PETSC
Matrix is symmetric.
void vmult(VectorBase &dst, const VectorBase &src) const
AdditionalData(const double pivoting=1.e-6, const double zero_pivot=1.e-12, const double damping=0.0)
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
Contents is actually a matrix.
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
void set_option_value(const std::string &name, const std::string &value)
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1419
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
std::string to_string(const T &t)
Definition: patterns.h:2360
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
void Tvmult(VectorBase &dst, const VectorBase &src) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
AdditionalData(const unsigned int symmetric=1, const unsigned int n_levels=1, const double threshold=0.1, const double filter=0.05, const bool output_details=false)
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
AdditionalData(const bool symmetric_operator=false, const double strong_threshold=0.25, const double max_row_sum=0.9, const unsigned int aggressive_coarsening_num_levels=0, const bool output_details=false)