Reference documentation for deal.II version Git 42dd89c428 2021-07-27 06:40:55 +0200
\(\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 - 2021 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  PreconditionBase::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  /* ----------------- PreconditionICC -------------------- */
311 
312 
314  : levels(levels)
315  {}
316 
317 
318 
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 *>(PCICC));
338  AssertThrow(ierr == 0, ExcPETScError(ierr));
339 
340  // then set flags
341  ierr = PCFactorSetLevels(pc, additional_data.levels);
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  /* ----------------- PreconditionILU -------------------- */
353 
355  : levels(levels)
356  {}
357 
358 
359 
362  {
363  initialize(matrix, additional_data);
364  }
365 
366 
367  void
369  const AdditionalData &additional_data_)
370  {
371  clear();
372 
373  matrix = static_cast<Mat>(matrix_);
374  additional_data = additional_data_;
375 
376  create_pc();
377 
378  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCILU));
379  AssertThrow(ierr == 0, ExcPETScError(ierr));
380 
381  // then set flags
382  ierr = PCFactorSetLevels(pc, additional_data.levels);
383  AssertThrow(ierr == 0, ExcPETScError(ierr));
384 
385  ierr = PCSetFromOptions(pc);
386  AssertThrow(ierr == 0, ExcPETScError(ierr));
387 
388  ierr = PCSetUp(pc);
389  AssertThrow(ierr == 0, ExcPETScError(ierr));
390  }
391 
392 
393  /* ----------------- PreconditionBoomerAMG -------------------- */
394 
396  const bool symmetric_operator,
397  const double strong_threshold,
398  const double max_row_sum,
399  const unsigned int aggressive_coarsening_num_levels,
400  const bool output_details,
401  const RelaxationType relaxation_type_up,
402  const RelaxationType relaxation_type_down,
403  const RelaxationType relaxation_type_coarse,
404  const unsigned int n_sweeps_coarse,
405  const double tol,
406  const unsigned int max_iter,
407  const bool w_cycle)
408  : symmetric_operator(symmetric_operator)
409  , strong_threshold(strong_threshold)
410  , max_row_sum(max_row_sum)
411  , aggressive_coarsening_num_levels(aggressive_coarsening_num_levels)
412  , output_details(output_details)
413  , relaxation_type_up(relaxation_type_up)
414  , relaxation_type_down(relaxation_type_down)
415  , relaxation_type_coarse(relaxation_type_coarse)
416  , n_sweeps_coarse(n_sweeps_coarse)
417  , tol(tol)
418  , max_iter(max_iter)
419  , w_cycle(w_cycle)
420  {}
421 
422 
423 
424 # ifdef DEAL_II_PETSC_WITH_HYPRE
425  namespace
426  {
431  std::string
432  to_string(
434  {
435  std::string string_type;
436 
437  switch (relaxation_type)
438  {
440  string_type = "Jacobi";
441  break;
444  string_type = "sequential-Gauss-Seidel";
445  break;
448  string_type = "seqboundary-Gauss-Seidel";
449  break;
451  string_type = "SOR/Jacobi";
452  break;
455  string_type = "backward-SOR/Jacobi";
456  break;
459  string_type = "symmetric-SOR/Jacobi";
460  break;
463  string_type = " l1scaled-SOR/Jacobi";
464  break;
467  string_type = "Gaussian-elimination";
468  break;
471  string_type = "l1-Gauss-Seidel";
472  break;
475  string_type = "backward-l1-Gauss-Seidel";
476  break;
478  string_type = "CG";
479  break;
481  string_type = "Chebyshev";
482  break;
484  string_type = "FCF-Jacobi";
485  break;
488  string_type = "l1scaled-Jacobi";
489  break;
491  string_type = "None";
492  break;
493  default:
494  Assert(false, ExcNotImplemented());
495  }
496  return string_type;
497  }
498  } // namespace
499 # endif
500 
502  const MPI_Comm & comm,
503  const AdditionalData &additional_data_)
504  {
505  additional_data = additional_data_;
506 
507  PetscErrorCode ierr = PCCreate(comm, &pc);
508  AssertThrow(ierr == 0, ExcPETScError(ierr));
509 
510 # ifdef DEAL_II_PETSC_WITH_HYPRE
511  initialize();
512 # else // DEAL_II_PETSC_WITH_HYPRE
513  (void)pc;
514  Assert(false,
515  ExcMessage("Your PETSc installation does not include a copy of "
516  "the hypre package necessary for this preconditioner."));
517 # endif
518  }
519 
520 
522  const MatrixBase & matrix,
524  {
525  initialize(matrix, additional_data);
526  }
527 
528  void
530  {
531 # ifdef DEAL_II_PETSC_WITH_HYPRE
532  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
533  AssertThrow(ierr == 0, ExcPETScError(ierr));
534 
535  ierr = PCHYPRESetType(pc, "boomeramg");
536  AssertThrow(ierr == 0, ExcPETScError(ierr));
537 
539  {
540  set_option_value("-pc_hypre_boomeramg_print_statistics", "1");
541  }
542 
543  set_option_value("-pc_hypre_boomeramg_agg_nl",
546 
547  set_option_value("-pc_hypre_boomeramg_max_row_sum",
549 
550  set_option_value("-pc_hypre_boomeramg_strong_threshold",
552 
553  // change to symmetric SOR/Jacobi when using a symmetric operator for
554  // backward compatibility
558  {
561  }
562 
566  {
569  }
570 
574  {
577  }
578 
579  auto relaxation_type_is_symmetric =
580  [](AdditionalData::RelaxationType relaxation_type) {
581  return relaxation_type == AdditionalData::RelaxationType::Jacobi ||
582  relaxation_type ==
584  relaxation_type ==
586  relaxation_type == AdditionalData::RelaxationType::None ||
587  relaxation_type ==
589  relaxation_type == AdditionalData::RelaxationType::CG ||
591  };
592 
594  !relaxation_type_is_symmetric(additional_data.relaxation_type_up))
595  Assert(false,
596  ExcMessage("Use a symmetric smoother for relaxation_type_up"));
597 
599  !relaxation_type_is_symmetric(additional_data.relaxation_type_down))
600  Assert(false,
601  ExcMessage("Use a symmetric smoother for relaxation_type_down"));
602 
604  !relaxation_type_is_symmetric(additional_data.relaxation_type_coarse))
605  Assert(false,
606  ExcMessage("Use a symmetric smoother for relaxation_type_coarse"));
607 
608  set_option_value("-pc_hypre_boomeramg_relax_type_up",
610  set_option_value("-pc_hypre_boomeramg_relax_type_down",
612  set_option_value("-pc_hypre_boomeramg_relax_type_coarse",
614  set_option_value("-pc_hypre_boomeramg_grid_sweeps_coarse",
616 
617  set_option_value("-pc_hypre_boomeramg_tol",
619  set_option_value("-pc_hypre_boomeramg_max_iter",
621 
623  {
624  set_option_value("-pc_hypre_boomeramg_cycle_type", "W");
625  }
626 
627  ierr = PCSetFromOptions(pc);
628  AssertThrow(ierr == 0, ExcPETScError(ierr));
629 # else
630  Assert(false,
631  ExcMessage("Your PETSc installation does not include a copy of "
632  "the hypre package necessary for this preconditioner."));
633 # endif
634  }
635 
636  void
638  const AdditionalData &additional_data_)
639  {
640 # ifdef DEAL_II_PETSC_WITH_HYPRE
641  clear();
642 
643  matrix = static_cast<Mat>(matrix_);
644  additional_data = additional_data_;
645 
646  create_pc();
647  initialize();
648 
649  PetscErrorCode ierr = PCSetUp(pc);
650  AssertThrow(ierr == 0, ExcPETScError(ierr));
651 
652 # else // DEAL_II_PETSC_WITH_HYPRE
653  (void)matrix_;
654  (void)additional_data_;
655  Assert(false,
656  ExcMessage("Your PETSc installation does not include a copy of "
657  "the hypre package necessary for this preconditioner."));
658 # endif
659  }
660 
661 
662  /* ----------------- PreconditionParaSails -------------------- */
663 
665  const unsigned int symmetric,
666  const unsigned int n_levels,
667  const double threshold,
668  const double filter,
669  const bool output_details)
670  : symmetric(symmetric)
671  , n_levels(n_levels)
672  , threshold(threshold)
673  , filter(filter)
674  , output_details(output_details)
675  {}
676 
677 
678 
680  const MatrixBase & matrix,
682  {
683  initialize(matrix, additional_data);
684  }
685 
686 
687  void
689  const AdditionalData &additional_data_)
690  {
691  clear();
692 
693  matrix = static_cast<Mat>(matrix_);
694  additional_data = additional_data_;
695 
696 # ifdef DEAL_II_PETSC_WITH_HYPRE
697  create_pc();
698 
699  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
700  AssertThrow(ierr == 0, ExcPETScError(ierr));
701 
702  ierr = PCHYPRESetType(pc, "parasails");
703  AssertThrow(ierr == 0, ExcPETScError(ierr));
704 
706  {
707  set_option_value("-pc_hypre_parasails_logging", "1");
708  }
709 
712  ExcMessage(
713  "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
714 
715  std::stringstream ssStream;
716 
717  switch (additional_data.symmetric)
718  {
719  case 0:
720  {
721  ssStream << "nonsymmetric";
722  break;
723  }
724 
725  case 1:
726  {
727  ssStream << "SPD";
728  break;
729  }
730 
731  case 2:
732  {
733  ssStream << "nonsymmetric,SPD";
734  break;
735  }
736 
737  default:
738  Assert(
739  false,
740  ExcMessage(
741  "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
742  }
743 
744  set_option_value("-pc_hypre_parasails_sym", ssStream.str());
745 
746  set_option_value("-pc_hypre_parasails_nlevels",
748 
749  ssStream.str(""); // empty the stringstream
750  ssStream << additional_data.threshold;
751  set_option_value("-pc_hypre_parasails_thresh", ssStream.str());
752 
753  ssStream.str(""); // empty the stringstream
754  ssStream << additional_data.filter;
755  set_option_value("-pc_hypre_parasails_filter", ssStream.str());
756 
757  ierr = PCSetFromOptions(pc);
758  AssertThrow(ierr == 0, ExcPETScError(ierr));
759 
760  ierr = PCSetUp(pc);
761  AssertThrow(ierr == 0, ExcPETScError(ierr));
762 
763 # else // DEAL_II_PETSC_WITH_HYPRE
764  (void)pc;
765  Assert(false,
766  ExcMessage("Your PETSc installation does not include a copy of "
767  "the hypre package necessary for this preconditioner."));
768 # endif
769  }
770 
771 
772  /* ----------------- PreconditionNone ------------------------- */
773 
776  {
777  initialize(matrix, additional_data);
778  }
779 
780 
781  void
783  const AdditionalData &additional_data_)
784  {
785  clear();
786 
787  matrix = static_cast<Mat>(matrix_);
788  additional_data = additional_data_;
789 
790  create_pc();
791 
792  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCNONE));
793  AssertThrow(ierr == 0, ExcPETScError(ierr));
794 
795  ierr = PCSetFromOptions(pc);
796  AssertThrow(ierr == 0, ExcPETScError(ierr));
797 
798  ierr = PCSetUp(pc);
799  AssertThrow(ierr == 0, ExcPETScError(ierr));
800  }
801 
802 
803  /* ----------------- PreconditionLU -------------------- */
804 
806  const double zero_pivot,
807  const double damping)
808  : pivoting(pivoting)
809  , zero_pivot(zero_pivot)
810  , damping(damping)
811  {}
812 
813 
814 
817  {
818  initialize(matrix, additional_data);
819  }
820 
821 
822  void
824  const AdditionalData &additional_data_)
825  {
826  clear();
827 
828  matrix = static_cast<Mat>(matrix_);
829  additional_data = additional_data_;
830 
831  create_pc();
832 
833  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCLU));
834  AssertThrow(ierr == 0, ExcPETScError(ierr));
835 
836  // set flags as given
837  ierr = PCFactorSetColumnPivot(pc, additional_data.pivoting);
838  AssertThrow(ierr == 0, ExcPETScError(ierr));
839 
840  ierr = PCFactorSetZeroPivot(pc, additional_data.zero_pivot);
841  AssertThrow(ierr == 0, ExcPETScError(ierr));
842 
843  ierr = PCFactorSetShiftAmount(pc, additional_data.damping);
844  AssertThrow(ierr == 0, ExcPETScError(ierr));
845 
846  ierr = PCSetFromOptions(pc);
847  AssertThrow(ierr == 0, ExcPETScError(ierr));
848 
849  ierr = PCSetUp(pc);
850  AssertThrow(ierr == 0, ExcPETScError(ierr));
851  }
852 
853 } // namespace PETScWrappers
854 
856 
857 #endif // DEAL_II_WITH_PETSC
Matrix is symmetric.
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())
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, const RelaxationType relaxation_type_up=RelaxationType::SORJacobi, const RelaxationType relaxation_type_down=RelaxationType::SORJacobi, const RelaxationType relaxation_type_coarse=RelaxationType::GaussianElimination, const unsigned int n_sweeps_coarse=1, const double tol=0.0, const unsigned int max_iter=1, const bool w_cycle=false)
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
void set_option_value(const std::string &name, const std::string &value)
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
std::string to_string(const T &t)
Definition: patterns.h:2329
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
void vmult(VectorBase &dst, const VectorBase &src) const
*braid_SplitCommworld & comm
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
void Tvmult(VectorBase &dst, const VectorBase &src) const
static ::ExceptionBase & ExcNotImplemented()
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())