Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
petsc_precondition.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2019 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 
16 #include <deal.II/lac/petsc_precondition.h>
17 
18 #ifdef DEAL_II_WITH_PETSC
19 
20 # include <deal.II/base/utilities.h>
21 
22 # include <deal.II/lac/exceptions.h>
23 # include <deal.II/lac/petsc_compatibility.h>
24 # include <deal.II/lac/petsc_matrix_base.h>
25 # include <deal.II/lac/petsc_solver.h>
26 # include <deal.II/lac/petsc_vector_base.h>
27 
28 # include <petscconf.h>
29 
30 # include <cmath>
31 
32 DEAL_II_NAMESPACE_OPEN
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  {
78  // only allow the creation of the
79  // preconditioner once
81 
82  MPI_Comm comm;
83  // this ugly cast is necessary because the
84  // type Mat and PETScObject are
85  // unrelated.
86  PetscErrorCode ierr =
87  PetscObjectGetComm(reinterpret_cast<PetscObject>(matrix), &comm);
88  AssertThrow(ierr == 0, ExcPETScError(ierr));
89 
90  ierr = PCCreate(comm, &pc);
91  AssertThrow(ierr == 0, ExcPETScError(ierr));
92 
93 # if DEAL_II_PETSC_VERSION_LT(3, 5, 0)
94  ierr = PCSetOperators(pc, matrix, matrix, SAME_PRECONDITIONER);
95 # else
96  ierr = PCSetOperators(pc, matrix, matrix);
97 # endif
98  AssertThrow(ierr == 0, ExcPETScError(ierr));
99  }
100 
101 
102  const PC &
104  {
105  return pc;
106  }
107 
108 
109  PreconditionerBase::operator Mat() const
110  {
111  return matrix;
112  }
113 
114 
115  /* ----------------- PreconditionJacobi -------------------- */
117  const AdditionalData &additional_data_)
118  {
119  additional_data = additional_data_;
120 
121  PetscErrorCode ierr = PCCreate(comm, &pc);
122  AssertThrow(ierr == 0, ExcPETScError(ierr));
123 
124  initialize();
125  }
126 
127 
128 
130  const AdditionalData &additional_data)
131  {
132  initialize(matrix, additional_data);
133  }
134 
135  void
137  {
139 
140  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCJACOBI));
141  AssertThrow(ierr == 0, ExcPETScError(ierr));
142 
143  ierr = PCSetFromOptions(pc);
144  AssertThrow(ierr == 0, ExcPETScError(ierr));
145  }
146 
147  void
149  const AdditionalData &additional_data_)
150  {
151  clear();
152 
153  matrix = static_cast<Mat>(matrix_);
154  additional_data = additional_data_;
155 
156  create_pc();
157  initialize();
158 
159  PetscErrorCode ierr = PCSetUp(pc);
160  AssertThrow(ierr == 0, ExcPETScError(ierr));
161  }
162 
163 
164  /* ----------------- PreconditionBlockJacobi -------------------- */
166  const MPI_Comm comm,
167  const AdditionalData &additional_data_)
168  {
169  additional_data = additional_data_;
170 
171  PetscErrorCode ierr = PCCreate(comm, &pc);
172  AssertThrow(ierr == 0, ExcPETScError(ierr));
173 
174  initialize();
175  }
176 
177 
178 
180  const MatrixBase & matrix,
181  const AdditionalData &additional_data)
182  {
183  initialize(matrix, additional_data);
184  }
185 
186  void
188  {
189  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCBJACOBI));
190  AssertThrow(ierr == 0, ExcPETScError(ierr));
191 
192  ierr = PCSetFromOptions(pc);
193  AssertThrow(ierr == 0, ExcPETScError(ierr));
194  }
195 
196 
197  void
199  const AdditionalData &additional_data_)
200  {
201  clear();
202 
203  matrix = static_cast<Mat>(matrix_);
204  additional_data = additional_data_;
205 
206  create_pc();
207  initialize();
208 
209  PetscErrorCode ierr = PCSetUp(pc);
210  AssertThrow(ierr == 0, ExcPETScError(ierr));
211  }
212 
213 
214  /* ----------------- PreconditionSOR -------------------- */
215 
217  : omega(omega)
218  {}
219 
220 
221 
224  {
225  initialize(matrix, additional_data);
226  }
227 
228 
229  void
231  const AdditionalData &additional_data_)
232  {
233  clear();
234 
235  matrix = static_cast<Mat>(matrix_);
236  additional_data = additional_data_;
237 
238  create_pc();
239 
240  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
241  AssertThrow(ierr == 0, ExcPETScError(ierr));
242 
243  // then set flags as given
244  ierr = PCSORSetOmega(pc, additional_data.omega);
245  AssertThrow(ierr == 0, ExcPETScError(ierr));
246 
247  ierr = PCSetFromOptions(pc);
248  AssertThrow(ierr == 0, ExcPETScError(ierr));
249 
250  ierr = PCSetUp(pc);
251  AssertThrow(ierr == 0, ExcPETScError(ierr));
252  }
253 
254 
255  /* ----------------- PreconditionSSOR -------------------- */
256 
258  : omega(omega)
259  {}
260 
261 
262 
265  {
266  initialize(matrix, additional_data);
267  }
268 
269 
270  void
272  const AdditionalData &additional_data_)
273  {
274  clear();
275 
276  matrix = static_cast<Mat>(matrix_);
277  additional_data = additional_data_;
278 
279  create_pc();
280 
281  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
282  AssertThrow(ierr == 0, ExcPETScError(ierr));
283 
284  // then set flags as given
285  ierr = PCSORSetOmega(pc, additional_data.omega);
286  AssertThrow(ierr == 0, ExcPETScError(ierr));
287 
288  // convert SOR to SSOR
289  ierr = PCSORSetSymmetric(pc, SOR_SYMMETRIC_SWEEP);
290  AssertThrow(ierr == 0, ExcPETScError(ierr));
291 
292  ierr = PCSetFromOptions(pc);
293  AssertThrow(ierr == 0, ExcPETScError(ierr));
294 
295  ierr = PCSetUp(pc);
296  AssertThrow(ierr == 0, ExcPETScError(ierr));
297  }
298 
299 
300  /* ----------------- PreconditionEisenstat -------------------- */
301 
303  : omega(omega)
304  {}
305 
306 
307 
309  const MatrixBase & matrix,
311  {
312  initialize(matrix, additional_data);
313  }
314 
315 
316  void
318  const AdditionalData &additional_data_)
319  {
320  clear();
321 
322  matrix = static_cast<Mat>(matrix_);
323  additional_data = additional_data_;
324 
325  create_pc();
326 
327  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCEISENSTAT));
328  AssertThrow(ierr == 0, ExcPETScError(ierr));
329 
330  // then set flags as given
331  ierr = PCEisenstatSetOmega(pc, additional_data.omega);
332  AssertThrow(ierr == 0, ExcPETScError(ierr));
333 
334  ierr = PCSetFromOptions(pc);
335  AssertThrow(ierr == 0, ExcPETScError(ierr));
336 
337  ierr = PCSetUp(pc);
338  AssertThrow(ierr == 0, ExcPETScError(ierr));
339  }
340 
341 
342  /* ----------------- PreconditionICC -------------------- */
343 
344 
346  : levels(levels)
347  {}
348 
349 
350 
353  {
354  initialize(matrix, additional_data);
355  }
356 
357 
358  void
360  const AdditionalData &additional_data_)
361  {
362  clear();
363 
364  matrix = static_cast<Mat>(matrix_);
365  additional_data = additional_data_;
366 
367  create_pc();
368 
369  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCICC));
370  AssertThrow(ierr == 0, ExcPETScError(ierr));
371 
372  // then set flags
373  ierr = PCFactorSetLevels(pc, additional_data.levels);
374  AssertThrow(ierr == 0, ExcPETScError(ierr));
375 
376  ierr = PCSetFromOptions(pc);
377  AssertThrow(ierr == 0, ExcPETScError(ierr));
378 
379  ierr = PCSetUp(pc);
380  AssertThrow(ierr == 0, ExcPETScError(ierr));
381  }
382 
383 
384  /* ----------------- PreconditionILU -------------------- */
385 
387  : levels(levels)
388  {}
389 
390 
391 
394  {
395  initialize(matrix, additional_data);
396  }
397 
398 
399  void
401  const AdditionalData &additional_data_)
402  {
403  clear();
404 
405  matrix = static_cast<Mat>(matrix_);
406  additional_data = additional_data_;
407 
408  create_pc();
409 
410  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCILU));
411  AssertThrow(ierr == 0, ExcPETScError(ierr));
412 
413  // then set flags
414  ierr = PCFactorSetLevels(pc, additional_data.levels);
415  AssertThrow(ierr == 0, ExcPETScError(ierr));
416 
417  ierr = PCSetFromOptions(pc);
418  AssertThrow(ierr == 0, ExcPETScError(ierr));
419 
420  ierr = PCSetUp(pc);
421  AssertThrow(ierr == 0, ExcPETScError(ierr));
422  }
423 
424 
425  /* ----------------- PreconditionBoomerAMG -------------------- */
426 
428  const bool symmetric_operator,
429  const double strong_threshold,
430  const double max_row_sum,
431  const unsigned int aggressive_coarsening_num_levels,
432  const bool output_details)
433  : symmetric_operator(symmetric_operator)
434  , strong_threshold(strong_threshold)
435  , max_row_sum(max_row_sum)
436  , aggressive_coarsening_num_levels(aggressive_coarsening_num_levels)
437  , output_details(output_details)
438  {}
439 
440 
441 
443  const MPI_Comm comm,
444  const AdditionalData &additional_data_)
445  {
446  additional_data = additional_data_;
447 
448  PetscErrorCode ierr = PCCreate(comm, &pc);
449  AssertThrow(ierr == 0, ExcPETScError(ierr));
450 
451 # ifdef DEAL_II_PETSC_WITH_HYPRE
452  initialize();
453 # else // DEAL_II_PETSC_WITH_HYPRE
454  (void)pc;
455  Assert(false,
456  ExcMessage("Your PETSc installation does not include a copy of "
457  "the hypre package necessary for this preconditioner."));
458 # endif
459  }
460 
461 
463  const MatrixBase & matrix,
465  {
466  initialize(matrix, additional_data);
467  }
468 
469  void
471  {
472 # ifdef DEAL_II_PETSC_WITH_HYPRE
473  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
474  AssertThrow(ierr == 0, ExcPETScError(ierr));
475 
476  ierr = PCHYPRESetType(pc, "boomeramg");
477  AssertThrow(ierr == 0, ExcPETScError(ierr));
478 
480  {
481  set_option_value("-pc_hypre_boomeramg_print_statistics", "1");
482  }
483 
484  set_option_value("-pc_hypre_boomeramg_agg_nl",
487 
488  std::stringstream ssStream;
489  ssStream << additional_data.max_row_sum;
490  set_option_value("-pc_hypre_boomeramg_max_row_sum", ssStream.str());
491 
492  ssStream.str(""); // empty the stringstream
493  ssStream << additional_data.strong_threshold;
494  set_option_value("-pc_hypre_boomeramg_strong_threshold", ssStream.str());
495 
497  {
498  set_option_value("-pc_hypre_boomeramg_relax_type_up",
499  "symmetric-SOR/Jacobi");
500  set_option_value("-pc_hypre_boomeramg_relax_type_down",
501  "symmetric-SOR/Jacobi");
502  set_option_value("-pc_hypre_boomeramg_relax_type_coarse",
503  "Gaussian-elimination");
504  }
505  else
506  {
507  set_option_value("-pc_hypre_boomeramg_relax_type_up", "SOR/Jacobi");
508  set_option_value("-pc_hypre_boomeramg_relax_type_down", "SOR/Jacobi");
509  set_option_value("-pc_hypre_boomeramg_relax_type_coarse",
510  "Gaussian-elimination");
511  }
512 
513  ierr = PCSetFromOptions(pc);
514  AssertThrow(ierr == 0, ExcPETScError(ierr));
515 # else
516  Assert(false,
517  ExcMessage("Your PETSc installation does not include a copy of "
518  "the hypre package necessary for this preconditioner."));
519 # endif
520  }
521 
522  void
524  const AdditionalData &additional_data_)
525  {
526 # ifdef DEAL_II_PETSC_WITH_HYPRE
527  clear();
528 
529  matrix = static_cast<Mat>(matrix_);
530  additional_data = additional_data_;
531 
532  create_pc();
533  initialize();
534 
535  PetscErrorCode ierr = PCSetUp(pc);
536  AssertThrow(ierr == 0, ExcPETScError(ierr));
537 
538 # else // DEAL_II_PETSC_WITH_HYPRE
539  (void)matrix_;
540  (void)additional_data_;
541  Assert(false,
542  ExcMessage("Your PETSc installation does not include a copy of "
543  "the hypre package necessary for this preconditioner."));
544 # endif
545  }
546 
547 
548  /* ----------------- PreconditionParaSails -------------------- */
549 
551  const unsigned int symmetric,
552  const unsigned int n_levels,
553  const double threshold,
554  const double filter,
555  const bool output_details)
556  : symmetric(symmetric)
557  , n_levels(n_levels)
558  , threshold(threshold)
559  , filter(filter)
560  , output_details(output_details)
561  {}
562 
563 
564 
566  const MatrixBase & matrix,
568  {
569  initialize(matrix, additional_data);
570  }
571 
572 
573  void
575  const AdditionalData &additional_data_)
576  {
577  clear();
578 
579  matrix = static_cast<Mat>(matrix_);
580  additional_data = additional_data_;
581 
582 # ifdef DEAL_II_PETSC_WITH_HYPRE
583  create_pc();
584 
585  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
586  AssertThrow(ierr == 0, ExcPETScError(ierr));
587 
588  ierr = PCHYPRESetType(pc, "parasails");
589  AssertThrow(ierr == 0, ExcPETScError(ierr));
590 
592  {
593  set_option_value("-pc_hypre_parasails_logging", "1");
594  }
595 
598  ExcMessage(
599  "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
600 
601  std::stringstream ssStream;
602 
603  switch (additional_data.symmetric)
604  {
605  case 0:
606  {
607  ssStream << "nonsymmetric";
608  break;
609  }
610 
611  case 1:
612  {
613  ssStream << "SPD";
614  break;
615  }
616 
617  case 2:
618  {
619  ssStream << "nonsymmetric,SPD";
620  break;
621  }
622 
623  default:
624  Assert(
625  false,
626  ExcMessage(
627  "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
628  }
629 
630  set_option_value("-pc_hypre_parasails_sym", ssStream.str());
631 
632  set_option_value("-pc_hypre_parasails_nlevels",
634 
635  ssStream.str(""); // empty the stringstream
636  ssStream << additional_data.threshold;
637  set_option_value("-pc_hypre_parasails_thresh", ssStream.str());
638 
639  ssStream.str(""); // empty the stringstream
640  ssStream << additional_data.filter;
641  set_option_value("-pc_hypre_parasails_filter", ssStream.str());
642 
643  ierr = PCSetFromOptions(pc);
644  AssertThrow(ierr == 0, ExcPETScError(ierr));
645 
646  ierr = PCSetUp(pc);
647  AssertThrow(ierr == 0, ExcPETScError(ierr));
648 
649 # else // DEAL_II_PETSC_WITH_HYPRE
650  (void)pc;
651  Assert(false,
652  ExcMessage("Your PETSc installation does not include a copy of "
653  "the hypre package necessary for this preconditioner."));
654 # endif
655  }
656 
657 
658  /* ----------------- PreconditionNone ------------------------- */
659 
662  {
663  initialize(matrix, additional_data);
664  }
665 
666 
667  void
669  const AdditionalData &additional_data_)
670  {
671  clear();
672 
673  matrix = static_cast<Mat>(matrix_);
674  additional_data = additional_data_;
675 
676  create_pc();
677 
678  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCNONE));
679  AssertThrow(ierr == 0, ExcPETScError(ierr));
680 
681  ierr = PCSetFromOptions(pc);
682  AssertThrow(ierr == 0, ExcPETScError(ierr));
683 
684  ierr = PCSetUp(pc);
685  AssertThrow(ierr == 0, ExcPETScError(ierr));
686  }
687 
688 
689  /* ----------------- PreconditionLU -------------------- */
690 
692  const double zero_pivot,
693  const double damping)
694  : pivoting(pivoting)
695  , zero_pivot(zero_pivot)
696  , damping(damping)
697  {}
698 
699 
700 
703  {
704  initialize(matrix, additional_data);
705  }
706 
707 
708  void
710  const AdditionalData &additional_data_)
711  {
712  clear();
713 
714  matrix = static_cast<Mat>(matrix_);
715  additional_data = additional_data_;
716 
717  create_pc();
718 
719  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCLU));
720  AssertThrow(ierr == 0, ExcPETScError(ierr));
721 
722  // set flags as given
723  ierr = PCFactorSetColumnPivot(pc, additional_data.pivoting);
724  AssertThrow(ierr == 0, ExcPETScError(ierr));
725 
726  ierr = PCFactorSetZeroPivot(pc, additional_data.zero_pivot);
727  AssertThrow(ierr == 0, ExcPETScError(ierr));
728 
729  ierr = PCFactorSetShiftAmount(pc, additional_data.damping);
730  AssertThrow(ierr == 0, ExcPETScError(ierr));
731 
732  ierr = PCSetFromOptions(pc);
733  AssertThrow(ierr == 0, ExcPETScError(ierr));
734 
735  ierr = PCSetUp(pc);
736  AssertThrow(ierr == 0, ExcPETScError(ierr));
737  }
738 
739 } // namespace PETScWrappers
740 
741 DEAL_II_NAMESPACE_CLOSE
742 
743 #endif // DEAL_II_WITH_PETSC
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())
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:1519
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:392
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:1407
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
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)