Reference documentation for deal.II version 9.3.2
\(\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 
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 {
36  PreconditionBase::PreconditionBase()
37  : pc(nullptr)
38  , matrix(nullptr)
39  {}
40 
41  PreconditionBase::~PreconditionBase()
42  {
43  try
44  {
45  clear();
46  }
47  catch (...)
48  {}
49  }
50 
51  void
52  PreconditionBase::clear()
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
66  PreconditionBase::vmult(VectorBase &dst, const VectorBase &src) const
67  {
68  AssertThrow(pc != nullptr, StandardExceptions::ExcInvalidState());
69 
70  const PetscErrorCode ierr = PCApply(pc, src, dst);
71  AssertThrow(ierr == 0, ExcPETScError(ierr));
72  }
73 
74 
75  void
76  PreconditionBase::Tvmult(VectorBase &dst, const VectorBase &src) const
77  {
78  AssertThrow(pc != nullptr, StandardExceptions::ExcInvalidState());
79 
80  const PetscErrorCode ierr = PCApplyTranspose(pc, src, dst);
81  AssertThrow(ierr == 0, ExcPETScError(ierr));
82  }
83 
84 
85  void
86  PreconditionBase::create_pc()
87  {
88  // only allow the creation of the
89  // preconditioner once
90  AssertThrow(pc == nullptr, StandardExceptions::ExcInvalidState());
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 &
113  PreconditionBase::get_pc() const
114  {
115  return pc;
116  }
117 
118 
119  PreconditionBase::operator Mat() const
120  {
121  return matrix;
122  }
123 
124 
125  /* ----------------- PreconditionJacobi -------------------- */
126  PreconditionJacobi::PreconditionJacobi(const MPI_Comm & comm,
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 
139  PreconditionJacobi::PreconditionJacobi(const MatrixBase & matrix,
140  const AdditionalData &additional_data)
141  {
142  initialize(matrix, additional_data);
143  }
144 
145  void
146  PreconditionJacobi::initialize()
147  {
148  AssertThrow(pc != nullptr, StandardExceptions::ExcInvalidState());
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
158  PreconditionJacobi::initialize(const MatrixBase & matrix_,
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 -------------------- */
175  PreconditionBlockJacobi::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 
189  PreconditionBlockJacobi::PreconditionBlockJacobi(
190  const MatrixBase & matrix,
191  const AdditionalData &additional_data)
192  {
193  initialize(matrix, additional_data);
194  }
195 
196  void
197  PreconditionBlockJacobi::initialize()
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
208  PreconditionBlockJacobi::initialize(const MatrixBase & matrix_,
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 
226  PreconditionSOR::AdditionalData::AdditionalData(const double omega)
227  : omega(omega)
228  {}
229 
230 
231 
232  PreconditionSOR::PreconditionSOR(const MatrixBase & matrix,
233  const AdditionalData &additional_data)
234  {
235  initialize(matrix, additional_data);
236  }
237 
238 
239  void
240  PreconditionSOR::initialize(const MatrixBase & matrix_,
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 
267  PreconditionSSOR::AdditionalData::AdditionalData(const double omega)
268  : omega(omega)
269  {}
270 
271 
272 
273  PreconditionSSOR::PreconditionSSOR(const MatrixBase & matrix,
274  const AdditionalData &additional_data)
275  {
276  initialize(matrix, additional_data);
277  }
278 
279 
280  void
281  PreconditionSSOR::initialize(const MatrixBase & matrix_,
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 
312  PreconditionEisenstat::AdditionalData::AdditionalData(const double omega)
313  : omega(omega)
314  {}
315 
316 
317 
318  PreconditionEisenstat::PreconditionEisenstat(
319  const MatrixBase & matrix,
320  const AdditionalData &additional_data)
321  {
322  initialize(matrix, additional_data);
323  }
324 
325 
326  void
327  PreconditionEisenstat::initialize(const MatrixBase & matrix_,
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 
355  PreconditionICC::AdditionalData::AdditionalData(const unsigned int levels)
356  : levels(levels)
357  {}
358 
359 
360 
361  PreconditionICC::PreconditionICC(const MatrixBase & matrix,
362  const AdditionalData &additional_data)
363  {
364  initialize(matrix, additional_data);
365  }
366 
367 
368  void
369  PreconditionICC::initialize(const MatrixBase & matrix_,
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 
396  PreconditionILU::AdditionalData::AdditionalData(const unsigned int levels)
397  : levels(levels)
398  {}
399 
400 
401 
402  PreconditionILU::PreconditionILU(const MatrixBase & matrix,
403  const AdditionalData &additional_data)
404  {
405  initialize(matrix, additional_data);
406  }
407 
408 
409  void
410  PreconditionILU::initialize(const MatrixBase & matrix_,
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 
437  PreconditionBoomerAMG::AdditionalData::AdditionalData(
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  const RelaxationType relaxation_type_up,
444  const RelaxationType relaxation_type_down,
445  const RelaxationType relaxation_type_coarse,
446  const unsigned int n_sweeps_coarse,
447  const double tol,
448  const unsigned int max_iter,
449  const bool w_cycle)
450  : symmetric_operator(symmetric_operator)
451  , strong_threshold(strong_threshold)
452  , max_row_sum(max_row_sum)
453  , aggressive_coarsening_num_levels(aggressive_coarsening_num_levels)
454  , output_details(output_details)
455  , relaxation_type_up(relaxation_type_up)
456  , relaxation_type_down(relaxation_type_down)
457  , relaxation_type_coarse(relaxation_type_coarse)
458  , n_sweeps_coarse(n_sweeps_coarse)
459  , tol(tol)
460  , max_iter(max_iter)
461  , w_cycle(w_cycle)
462  {}
463 
464 
465 
466 # ifdef DEAL_II_PETSC_WITH_HYPRE
467  namespace
468  {
473  std::string
474  to_string(
475  PreconditionBoomerAMG::AdditionalData::RelaxationType relaxation_type)
476  {
477  std::string string_type;
478 
479  switch (relaxation_type)
480  {
481  case PreconditionBoomerAMG::AdditionalData::RelaxationType::Jacobi:
482  string_type = "Jacobi";
483  break;
484  case PreconditionBoomerAMG::AdditionalData::RelaxationType::
485  sequentialGaussSeidel:
486  string_type = "sequential-Gauss-Seidel";
487  break;
488  case PreconditionBoomerAMG::AdditionalData::RelaxationType::
489  seqboundaryGaussSeidel:
490  string_type = "seqboundary-Gauss-Seidel";
491  break;
492  case PreconditionBoomerAMG::AdditionalData::RelaxationType::SORJacobi:
493  string_type = "SOR/Jacobi";
494  break;
495  case PreconditionBoomerAMG::AdditionalData::RelaxationType::
496  backwardSORJacobi:
497  string_type = "backward-SOR/Jacobi";
498  break;
499  case PreconditionBoomerAMG::AdditionalData::RelaxationType::
500  symmetricSORJacobi:
501  string_type = "symmetric-SOR/Jacobi";
502  break;
503  case PreconditionBoomerAMG::AdditionalData::RelaxationType::
504  l1scaledSORJacobi:
505  string_type = " l1scaled-SOR/Jacobi";
506  break;
507  case PreconditionBoomerAMG::AdditionalData::RelaxationType::
508  GaussianElimination:
509  string_type = "Gaussian-elimination";
510  break;
511  case PreconditionBoomerAMG::AdditionalData::RelaxationType::
512  l1GaussSeidel:
513  string_type = "l1-Gauss-Seidel";
514  break;
515  case PreconditionBoomerAMG::AdditionalData::RelaxationType::
516  backwardl1GaussSeidel:
517  string_type = "backward-l1-Gauss-Seidel";
518  break;
519  case PreconditionBoomerAMG::AdditionalData::RelaxationType::CG:
520  string_type = "CG";
521  break;
522  case PreconditionBoomerAMG::AdditionalData::RelaxationType::Chebyshev:
523  string_type = "Chebyshev";
524  break;
525  case PreconditionBoomerAMG::AdditionalData::RelaxationType::FCFJacobi:
526  string_type = "FCF-Jacobi";
527  break;
528  case PreconditionBoomerAMG::AdditionalData::RelaxationType::
529  l1scaledJacobi:
530  string_type = "l1scaled-Jacobi";
531  break;
532  case PreconditionBoomerAMG::AdditionalData::RelaxationType::None:
533  string_type = "None";
534  break;
535  default:
536  Assert(false, ExcNotImplemented());
537  }
538  return string_type;
539  }
540  } // namespace
541 # endif
542 
543  PreconditionBoomerAMG::PreconditionBoomerAMG(
544  const MPI_Comm & comm,
545  const AdditionalData &additional_data_)
546  {
547  additional_data = additional_data_;
548 
549  PetscErrorCode ierr = PCCreate(comm, &pc);
550  AssertThrow(ierr == 0, ExcPETScError(ierr));
551 
552 # ifdef DEAL_II_PETSC_WITH_HYPRE
553  initialize();
554 # else // DEAL_II_PETSC_WITH_HYPRE
555  (void)pc;
556  Assert(false,
557  ExcMessage("Your PETSc installation does not include a copy of "
558  "the hypre package necessary for this preconditioner."));
559 # endif
560  }
561 
562 
563  PreconditionBoomerAMG::PreconditionBoomerAMG(
564  const MatrixBase & matrix,
565  const AdditionalData &additional_data)
566  {
567  initialize(matrix, additional_data);
568  }
569 
570  void
571  PreconditionBoomerAMG::initialize()
572  {
573 # ifdef DEAL_II_PETSC_WITH_HYPRE
574  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
575  AssertThrow(ierr == 0, ExcPETScError(ierr));
576 
577  ierr = PCHYPRESetType(pc, "boomeramg");
578  AssertThrow(ierr == 0, ExcPETScError(ierr));
579 
580  if (additional_data.output_details)
581  {
582  set_option_value("-pc_hypre_boomeramg_print_statistics", "1");
583  }
584 
585  set_option_value("-pc_hypre_boomeramg_agg_nl",
586  std::to_string(
587  additional_data.aggressive_coarsening_num_levels));
588 
589  set_option_value("-pc_hypre_boomeramg_max_row_sum",
590  std::to_string(additional_data.max_row_sum));
591 
592  set_option_value("-pc_hypre_boomeramg_strong_threshold",
593  std::to_string(additional_data.strong_threshold));
594 
595  // change to symmetric SOR/Jacobi when using a symmetric operator for
596  // backward compatibility
597  if (additional_data.relaxation_type_up ==
598  AdditionalData::RelaxationType::SORJacobi &&
599  additional_data.symmetric_operator)
600  {
601  additional_data.relaxation_type_up =
602  AdditionalData::RelaxationType::symmetricSORJacobi;
603  }
604 
605  if (additional_data.relaxation_type_down ==
606  AdditionalData::RelaxationType::SORJacobi &&
607  additional_data.symmetric_operator)
608  {
609  additional_data.relaxation_type_down =
610  AdditionalData::RelaxationType::symmetricSORJacobi;
611  }
612 
613  if (additional_data.relaxation_type_coarse ==
614  AdditionalData::RelaxationType::SORJacobi &&
615  additional_data.symmetric_operator)
616  {
617  additional_data.relaxation_type_coarse =
618  AdditionalData::RelaxationType::symmetricSORJacobi;
619  }
620 
621  auto relaxation_type_is_symmetric =
622  [](AdditionalData::RelaxationType relaxation_type) {
623  return relaxation_type == AdditionalData::RelaxationType::Jacobi ||
624  relaxation_type ==
625  AdditionalData::RelaxationType::symmetricSORJacobi ||
626  relaxation_type ==
627  AdditionalData::RelaxationType::GaussianElimination ||
628  relaxation_type == AdditionalData::RelaxationType::None ||
629  relaxation_type ==
630  AdditionalData::RelaxationType::l1scaledJacobi ||
631  relaxation_type == AdditionalData::RelaxationType::CG ||
632  relaxation_type == AdditionalData::RelaxationType::Chebyshev;
633  };
634 
635  if (additional_data.symmetric_operator &&
636  !relaxation_type_is_symmetric(additional_data.relaxation_type_up))
637  Assert(false,
638  ExcMessage("Use a symmetric smoother for relaxation_type_up"));
639 
640  if (additional_data.symmetric_operator &&
641  !relaxation_type_is_symmetric(additional_data.relaxation_type_down))
642  Assert(false,
643  ExcMessage("Use a symmetric smoother for relaxation_type_down"));
644 
645  if (additional_data.symmetric_operator &&
646  !relaxation_type_is_symmetric(additional_data.relaxation_type_coarse))
647  Assert(false,
648  ExcMessage("Use a symmetric smoother for relaxation_type_coarse"));
649 
650  set_option_value("-pc_hypre_boomeramg_relax_type_up",
651  to_string(additional_data.relaxation_type_up));
652  set_option_value("-pc_hypre_boomeramg_relax_type_down",
653  to_string(additional_data.relaxation_type_down));
654  set_option_value("-pc_hypre_boomeramg_relax_type_coarse",
655  to_string(additional_data.relaxation_type_coarse));
656  set_option_value("-pc_hypre_boomeramg_grid_sweeps_coarse",
657  std::to_string(additional_data.n_sweeps_coarse));
658 
659  set_option_value("-pc_hypre_boomeramg_tol",
660  std::to_string(additional_data.tol));
661  set_option_value("-pc_hypre_boomeramg_max_iter",
662  std::to_string(additional_data.max_iter));
663 
664  if (additional_data.w_cycle)
665  {
666  set_option_value("-pc_hypre_boomeramg_cycle_type", "W");
667  }
668 
669  ierr = PCSetFromOptions(pc);
670  AssertThrow(ierr == 0, ExcPETScError(ierr));
671 # else
672  Assert(false,
673  ExcMessage("Your PETSc installation does not include a copy of "
674  "the hypre package necessary for this preconditioner."));
675 # endif
676  }
677 
678  void
679  PreconditionBoomerAMG::initialize(const MatrixBase & matrix_,
680  const AdditionalData &additional_data_)
681  {
682 # ifdef DEAL_II_PETSC_WITH_HYPRE
683  clear();
684 
685  matrix = static_cast<Mat>(matrix_);
686  additional_data = additional_data_;
687 
688  create_pc();
689  initialize();
690 
691  PetscErrorCode ierr = PCSetUp(pc);
692  AssertThrow(ierr == 0, ExcPETScError(ierr));
693 
694 # else // DEAL_II_PETSC_WITH_HYPRE
695  (void)matrix_;
696  (void)additional_data_;
697  Assert(false,
698  ExcMessage("Your PETSc installation does not include a copy of "
699  "the hypre package necessary for this preconditioner."));
700 # endif
701  }
702 
703 
704  /* ----------------- PreconditionParaSails -------------------- */
705 
706  PreconditionParaSails::AdditionalData::AdditionalData(
707  const unsigned int symmetric,
708  const unsigned int n_levels,
709  const double threshold,
710  const double filter,
711  const bool output_details)
712  : symmetric(symmetric)
713  , n_levels(n_levels)
714  , threshold(threshold)
715  , filter(filter)
716  , output_details(output_details)
717  {}
718 
719 
720 
721  PreconditionParaSails::PreconditionParaSails(
722  const MatrixBase & matrix,
723  const AdditionalData &additional_data)
724  {
725  initialize(matrix, additional_data);
726  }
727 
728 
729  void
730  PreconditionParaSails::initialize(const MatrixBase & matrix_,
731  const AdditionalData &additional_data_)
732  {
733  clear();
734 
735  matrix = static_cast<Mat>(matrix_);
736  additional_data = additional_data_;
737 
738 # ifdef DEAL_II_PETSC_WITH_HYPRE
739  create_pc();
740 
741  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
742  AssertThrow(ierr == 0, ExcPETScError(ierr));
743 
744  ierr = PCHYPRESetType(pc, "parasails");
745  AssertThrow(ierr == 0, ExcPETScError(ierr));
746 
747  if (additional_data.output_details)
748  {
749  set_option_value("-pc_hypre_parasails_logging", "1");
750  }
751 
752  Assert((additional_data.symmetric == 0 || additional_data.symmetric == 1 ||
753  additional_data.symmetric == 2),
754  ExcMessage(
755  "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
756 
757  std::stringstream ssStream;
758 
759  switch (additional_data.symmetric)
760  {
761  case 0:
762  {
763  ssStream << "nonsymmetric";
764  break;
765  }
766 
767  case 1:
768  {
769  ssStream << "SPD";
770  break;
771  }
772 
773  case 2:
774  {
775  ssStream << "nonsymmetric,SPD";
776  break;
777  }
778 
779  default:
780  Assert(
781  false,
782  ExcMessage(
783  "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
784  }
785 
786  set_option_value("-pc_hypre_parasails_sym", ssStream.str());
787 
788  set_option_value("-pc_hypre_parasails_nlevels",
789  std::to_string(additional_data.n_levels));
790 
791  ssStream.str(""); // empty the stringstream
792  ssStream << additional_data.threshold;
793  set_option_value("-pc_hypre_parasails_thresh", ssStream.str());
794 
795  ssStream.str(""); // empty the stringstream
796  ssStream << additional_data.filter;
797  set_option_value("-pc_hypre_parasails_filter", ssStream.str());
798 
799  ierr = PCSetFromOptions(pc);
800  AssertThrow(ierr == 0, ExcPETScError(ierr));
801 
802  ierr = PCSetUp(pc);
803  AssertThrow(ierr == 0, ExcPETScError(ierr));
804 
805 # else // DEAL_II_PETSC_WITH_HYPRE
806  (void)pc;
807  Assert(false,
808  ExcMessage("Your PETSc installation does not include a copy of "
809  "the hypre package necessary for this preconditioner."));
810 # endif
811  }
812 
813 
814  /* ----------------- PreconditionNone ------------------------- */
815 
816  PreconditionNone::PreconditionNone(const MatrixBase & matrix,
817  const AdditionalData &additional_data)
818  {
819  initialize(matrix, additional_data);
820  }
821 
822 
823  void
824  PreconditionNone::initialize(const MatrixBase & matrix_,
825  const AdditionalData &additional_data_)
826  {
827  clear();
828 
829  matrix = static_cast<Mat>(matrix_);
830  additional_data = additional_data_;
831 
832  create_pc();
833 
834  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCNONE));
835  AssertThrow(ierr == 0, ExcPETScError(ierr));
836 
837  ierr = PCSetFromOptions(pc);
838  AssertThrow(ierr == 0, ExcPETScError(ierr));
839 
840  ierr = PCSetUp(pc);
841  AssertThrow(ierr == 0, ExcPETScError(ierr));
842  }
843 
844 
845  /* ----------------- PreconditionLU -------------------- */
846 
847  PreconditionLU::AdditionalData::AdditionalData(const double pivoting,
848  const double zero_pivot,
849  const double damping)
850  : pivoting(pivoting)
851  , zero_pivot(zero_pivot)
852  , damping(damping)
853  {}
854 
855 
856 
857  PreconditionLU::PreconditionLU(const MatrixBase & matrix,
858  const AdditionalData &additional_data)
859  {
860  initialize(matrix, additional_data);
861  }
862 
863 
864  void
865  PreconditionLU::initialize(const MatrixBase & matrix_,
866  const AdditionalData &additional_data_)
867  {
868  clear();
869 
870  matrix = static_cast<Mat>(matrix_);
871  additional_data = additional_data_;
872 
873  create_pc();
874 
875  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCLU));
876  AssertThrow(ierr == 0, ExcPETScError(ierr));
877 
878  // set flags as given
879  ierr = PCFactorSetColumnPivot(pc, additional_data.pivoting);
880  AssertThrow(ierr == 0, ExcPETScError(ierr));
881 
882  ierr = PCFactorSetZeroPivot(pc, additional_data.zero_pivot);
883  AssertThrow(ierr == 0, ExcPETScError(ierr));
884 
885  ierr = PCFactorSetShiftAmount(pc, additional_data.damping);
886  AssertThrow(ierr == 0, ExcPETScError(ierr));
887 
888  ierr = PCSetFromOptions(pc);
889  AssertThrow(ierr == 0, ExcPETScError(ierr));
890 
891  ierr = PCSetUp(pc);
892  AssertThrow(ierr == 0, ExcPETScError(ierr));
893  }
894 
895 } // namespace PETScWrappers
896 
897 DEAL_II_NAMESPACE_CLOSE
898 
899 #endif // DEAL_II_WITH_PETSC
PETScWrappers::PreconditionBase::Tvmult
void Tvmult(VectorBase &dst, const VectorBase &src) const
Definition: petsc_precondition.cc:76
PETScWrappers::PreconditionBase::vmult
void vmult(VectorBase &dst, const VectorBase &src) const
Definition: petsc_precondition.cc:66
PETScWrappers::PreconditionEisenstat::initialize
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
Definition: petsc_precondition.cc:327
PETScWrappers::PreconditionICC::initialize
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
Definition: petsc_precondition.cc:369
PETScWrappers::PreconditionILU::initialize
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
Definition: petsc_precondition.cc:410
PETScWrappers::PreconditionLU::initialize
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
Definition: petsc_precondition.cc:865
PETScWrappers::PreconditionNone::initialize
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
Definition: petsc_precondition.cc:824
PETScWrappers::PreconditionParaSails::initialize
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
Definition: petsc_precondition.cc:730
PETScWrappers::PreconditionSOR::initialize
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
Definition: petsc_precondition.cc:240
PETScWrappers::PreconditionSSOR::initialize
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
Definition: petsc_precondition.cc:281
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:396
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:397
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1465
Patterns::Tools::to_string
std::string to_string(const T &t)
Definition: patterns.h:2329
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
StandardExceptions::ExcInvalidState
static ::ExceptionBase & ExcInvalidState()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:58
LAPACKSupport::symmetric
@ symmetric
Matrix is symmetric.
Definition: lapack_support.h:113
PETScWrappers::set_option_value
void set_option_value(const std::string &name, const std::string &value)
Definition: petsc_compatibility.h:46
petsc_compatibility.h
petsc_matrix_base.h
petsc_precondition.h
petsc_solver.h
petsc_vector_base.h
PETScWrappers::PreconditionBoomerAMG::AdditionalData::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)
Definition: petsc_precondition.cc:437
PETScWrappers::PreconditionLU::AdditionalData::AdditionalData
AdditionalData(const double pivoting=1.e-6, const double zero_pivot=1.e-12, const double damping=0.0)
Definition: petsc_precondition.cc:847
PETScWrappers::PreconditionParaSails::AdditionalData::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)
Definition: petsc_precondition.cc:706
Generated by doxygen 1.9.1