Reference documentation for deal.II version GIT 01a9543f1b 2023-12-05 20:40:02+00:00
\(\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 - 2023 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  {
40  }
41 
43  : pc(nullptr)
44  {}
45 
47  {
48  try
49  {
50  clear();
51  }
52  catch (...)
53  {}
54  }
55 
56  void
58  {
59  if (pc)
60  {
61  PetscErrorCode ierr = PCDestroy(&pc);
62  AssertThrow(ierr == 0, ExcPETScError(ierr));
63  }
64  }
65 
66  void
68  {
70 
71  PetscErrorCode ierr = PCApply(pc, src, dst);
72  AssertThrow(ierr == 0, ExcPETScError(ierr));
73  }
74 
75  void
77  {
79 
80  PetscErrorCode ierr = PCApplyTranspose(pc, src, dst);
81  AssertThrow(ierr == 0, ExcPETScError(ierr));
82  }
83 
84  void
86  {
88 
89  PetscErrorCode ierr = PCSetUp(pc);
90  AssertThrow(ierr == 0, ExcPETScError(ierr));
91  }
92 
93  MPI_Comm
95  {
96  return PetscObjectComm(reinterpret_cast<PetscObject>(pc));
97  }
98 
99  void
101  {
102  // only allow the creation of the
103  // preconditioner once
105 
106  MPI_Comm comm;
107  PetscErrorCode ierr = PetscObjectGetComm(
108  reinterpret_cast<PetscObject>(static_cast<const Mat &>(matrix)), &comm);
109  AssertThrow(ierr == 0, ExcPETScError(ierr));
110 
112 
113  ierr = PCSetOperators(pc, matrix, matrix);
114  AssertThrow(ierr == 0, ExcPETScError(ierr));
115  }
116 
117  void
119  {
120  clear();
121  PetscErrorCode ierr = PCCreate(comm, &pc);
122  AssertThrow(ierr == 0, ExcPETScError(ierr));
123  }
124 
125  const PC &
127  {
128  return pc;
129  }
130 
131 
132  /* ----------------- PreconditionJacobi -------------------- */
133 
135  : PreconditionBase()
136  {}
137 
138 
139 
141  const AdditionalData &additional_data_)
143  {
144  additional_data = additional_data_;
145 
146  PetscErrorCode ierr = PCCreate(comm, &pc);
147  AssertThrow(ierr == 0, ExcPETScError(ierr));
148 
149  initialize();
150  }
151 
152 
153 
155  const AdditionalData &additional_data)
156  : PreconditionBase(matrix.get_mpi_communicator())
157  {
159  }
160 
161 
162 
163  void
165  {
167 
168  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCJACOBI));
169  AssertThrow(ierr == 0, ExcPETScError(ierr));
170 
171  ierr = PCSetFromOptions(pc);
172  AssertThrow(ierr == 0, ExcPETScError(ierr));
173  }
174 
175 
176 
177  void
179  const AdditionalData &additional_data_)
180  {
181  clear();
182 
183  additional_data = additional_data_;
184 
185  create_pc_with_mat(matrix_);
186  initialize();
187  }
188 
189 
190  /* ----------------- PreconditionBlockJacobi -------------------- */
191 
193  : PreconditionBase()
194  {}
195 
197  const MPI_Comm comm,
198  const AdditionalData &additional_data_)
200  {
201  additional_data = additional_data_;
202 
203  PetscErrorCode ierr = PCCreate(comm, &pc);
204  AssertThrow(ierr == 0, ExcPETScError(ierr));
205 
206  initialize();
207  }
208 
209 
210 
212  const MatrixBase &matrix,
213  const AdditionalData &additional_data)
214  : PreconditionBase(matrix.get_mpi_communicator())
215  {
217  }
218 
219 
220 
221  void
223  {
224  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCBJACOBI));
225  AssertThrow(ierr == 0, ExcPETScError(ierr));
226 
227  ierr = PCSetFromOptions(pc);
228  AssertThrow(ierr == 0, ExcPETScError(ierr));
229  }
230 
231 
232 
233  void
235  const AdditionalData &additional_data_)
236  {
237  clear();
238 
239  additional_data = additional_data_;
240 
241  create_pc_with_mat(matrix_);
242  initialize();
243  }
244 
245 
246  /* ----------------- PreconditionSOR -------------------- */
247 
249  : PreconditionBase()
250  {}
251 
252 
253 
255  : omega(omega)
256  {}
257 
258 
259 
263  {
265  }
266 
267 
268  void
270  const AdditionalData &additional_data_)
271  {
272  clear();
273 
274  additional_data = additional_data_;
275 
276  create_pc_with_mat(matrix_);
277 
278  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
279  AssertThrow(ierr == 0, ExcPETScError(ierr));
280 
281  // then set flags as given
282  ierr = PCSORSetOmega(pc, additional_data.omega);
283  AssertThrow(ierr == 0, ExcPETScError(ierr));
284 
285  ierr = PCSetFromOptions(pc);
286  AssertThrow(ierr == 0, ExcPETScError(ierr));
287  }
288 
289 
290  /* ----------------- PreconditionSSOR -------------------- */
291 
293  : PreconditionBase()
294  {}
295 
296 
297 
299  : omega(omega)
300  {}
301 
302 
303 
307  {
309  }
310 
311 
312  void
314  const AdditionalData &additional_data_)
315  {
316  clear();
317 
318  additional_data = additional_data_;
319 
320  create_pc_with_mat(matrix_);
321 
322  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
323  AssertThrow(ierr == 0, ExcPETScError(ierr));
324 
325  // then set flags as given
326  ierr = PCSORSetOmega(pc, additional_data.omega);
327  AssertThrow(ierr == 0, ExcPETScError(ierr));
328 
329  // convert SOR to SSOR
330  ierr = PCSORSetSymmetric(pc, SOR_SYMMETRIC_SWEEP);
331  AssertThrow(ierr == 0, ExcPETScError(ierr));
332 
333  ierr = PCSetFromOptions(pc);
334  AssertThrow(ierr == 0, ExcPETScError(ierr));
335  }
336 
337 
338  /* ----------------- PreconditionICC -------------------- */
339 
341  : PreconditionBase()
342  {}
343 
344 
345 
347  : levels(levels)
348  {}
349 
350 
351 
355  {
357  }
358 
359 
360  void
362  const AdditionalData &additional_data_)
363  {
364  clear();
365 
366  additional_data = additional_data_;
367 
368  create_pc_with_mat(matrix_);
369 
370  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCICC));
371  AssertThrow(ierr == 0, ExcPETScError(ierr));
372 
373  // then set flags
374  ierr = PCFactorSetLevels(pc, additional_data.levels);
375  AssertThrow(ierr == 0, ExcPETScError(ierr));
376 
377  ierr = PCSetFromOptions(pc);
378  AssertThrow(ierr == 0, ExcPETScError(ierr));
379  }
380 
381 
382  /* ----------------- PreconditionILU -------------------- */
383 
385  : PreconditionBase()
386  {}
387 
388 
389 
391  : levels(levels)
392  {}
393 
394 
395 
399  {
401  }
402 
403 
404  void
406  const AdditionalData &additional_data_)
407  {
408  clear();
409 
410  additional_data = additional_data_;
411 
412  create_pc_with_mat(matrix_);
413 
414  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCILU));
415  AssertThrow(ierr == 0, ExcPETScError(ierr));
416 
417  // then set flags
418  ierr = PCFactorSetLevels(pc, additional_data.levels);
419  AssertThrow(ierr == 0, ExcPETScError(ierr));
420 
421  ierr = PCSetFromOptions(pc);
422  AssertThrow(ierr == 0, ExcPETScError(ierr));
423  }
424 
425 
426  /* ----------------- PreconditionBoomerAMG -------------------- */
427 
429  const bool symmetric_operator,
430  const double strong_threshold,
431  const double max_row_sum,
432  const unsigned int aggressive_coarsening_num_levels,
433  const bool output_details,
434  const RelaxationType relaxation_type_up,
435  const RelaxationType relaxation_type_down,
436  const RelaxationType relaxation_type_coarse,
437  const unsigned int n_sweeps_coarse,
438  const double tol,
439  const unsigned int max_iter,
440  const bool w_cycle)
441  : symmetric_operator(symmetric_operator)
442  , strong_threshold(strong_threshold)
443  , max_row_sum(max_row_sum)
444  , aggressive_coarsening_num_levels(aggressive_coarsening_num_levels)
445  , output_details(output_details)
446  , relaxation_type_up(relaxation_type_up)
447  , relaxation_type_down(relaxation_type_down)
448  , relaxation_type_coarse(relaxation_type_coarse)
449  , n_sweeps_coarse(n_sweeps_coarse)
450  , tol(tol)
451  , max_iter(max_iter)
452  , w_cycle(w_cycle)
453  {}
454 
455 
456 
457 # ifdef DEAL_II_PETSC_WITH_HYPRE
458  namespace
459  {
464  std::string
465  to_string(
467  {
468  std::string string_type;
469 
470  switch (relaxation_type)
471  {
473  string_type = "Jacobi";
474  break;
477  string_type = "sequential-Gauss-Seidel";
478  break;
481  string_type = "seqboundary-Gauss-Seidel";
482  break;
484  string_type = "SOR/Jacobi";
485  break;
488  string_type = "backward-SOR/Jacobi";
489  break;
492  string_type = "symmetric-SOR/Jacobi";
493  break;
496  string_type = " l1scaled-SOR/Jacobi";
497  break;
500  string_type = "Gaussian-elimination";
501  break;
504  string_type = "l1-Gauss-Seidel";
505  break;
508  string_type = "backward-l1-Gauss-Seidel";
509  break;
511  string_type = "CG";
512  break;
514  string_type = "Chebyshev";
515  break;
517  string_type = "FCF-Jacobi";
518  break;
521  string_type = "l1scaled-Jacobi";
522  break;
524  string_type = "None";
525  break;
526  default:
527  Assert(false, ExcNotImplemented());
528  }
529  return string_type;
530  }
531  } // namespace
532 # endif
533 
534 
535 
537  : PreconditionBase()
538  {}
539 
540 
541 
543  const MPI_Comm comm,
544  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 
565  const MatrixBase &matrix,
566  const AdditionalData &additional_data)
567  : PreconditionBase(matrix.get_mpi_communicator())
568  {
570  }
571 
572 
573 
574  void
576  {
577 # ifdef DEAL_II_PETSC_WITH_HYPRE
578  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
579  AssertThrow(ierr == 0, ExcPETScError(ierr));
580 
581  ierr = PCHYPRESetType(pc, "boomeramg");
582  AssertThrow(ierr == 0, ExcPETScError(ierr));
583 
585  {
586  set_option_value("-pc_hypre_boomeramg_print_statistics", "1");
587  }
588 
589  set_option_value("-pc_hypre_boomeramg_agg_nl",
592 
593  set_option_value("-pc_hypre_boomeramg_max_row_sum",
595 
596  set_option_value("-pc_hypre_boomeramg_strong_threshold",
598 
599  // change to symmetric SOR/Jacobi when using a symmetric operator for
600  // backward compatibility
604  {
607  }
608 
612  {
615  }
616 
620  {
623  }
624 
625  auto relaxation_type_is_symmetric =
626  [](AdditionalData::RelaxationType relaxation_type) {
627  return relaxation_type == AdditionalData::RelaxationType::Jacobi ||
628  relaxation_type ==
630  relaxation_type ==
632  relaxation_type == AdditionalData::RelaxationType::None ||
633  relaxation_type ==
635  relaxation_type == AdditionalData::RelaxationType::CG ||
637  };
638 
640  !relaxation_type_is_symmetric(additional_data.relaxation_type_up))
641  Assert(false,
642  ExcMessage("Use a symmetric smoother for relaxation_type_up"));
643 
645  !relaxation_type_is_symmetric(additional_data.relaxation_type_down))
646  Assert(false,
647  ExcMessage("Use a symmetric smoother for relaxation_type_down"));
648 
650  !relaxation_type_is_symmetric(additional_data.relaxation_type_coarse))
651  Assert(false,
652  ExcMessage("Use a symmetric smoother for relaxation_type_coarse"));
653 
654  set_option_value("-pc_hypre_boomeramg_relax_type_up",
656  set_option_value("-pc_hypre_boomeramg_relax_type_down",
658  set_option_value("-pc_hypre_boomeramg_relax_type_coarse",
660  set_option_value("-pc_hypre_boomeramg_grid_sweeps_coarse",
662 
663  set_option_value("-pc_hypre_boomeramg_tol",
665  set_option_value("-pc_hypre_boomeramg_max_iter",
667 
669  {
670  set_option_value("-pc_hypre_boomeramg_cycle_type", "W");
671  }
672 
673  ierr = PCSetFromOptions(pc);
674  AssertThrow(ierr == 0, ExcPETScError(ierr));
675 # else
676  Assert(false,
677  ExcMessage("Your PETSc installation does not include a copy of "
678  "the hypre package necessary for this preconditioner."));
679 # endif
680  }
681 
682 
683 
684  void
686  const AdditionalData &additional_data_)
687  {
688 # ifdef DEAL_II_PETSC_WITH_HYPRE
689  clear();
690 
691  additional_data = additional_data_;
692 
693  create_pc_with_mat(matrix_);
694  initialize();
695 
696 # else // DEAL_II_PETSC_WITH_HYPRE
697  (void)matrix_;
698  (void)additional_data_;
699  Assert(false,
700  ExcMessage("Your PETSc installation does not include a copy of "
701  "the hypre package necessary for this preconditioner."));
702 # endif
703  }
704 
705 
706  /* ----------------- PreconditionParaSails -------------------- */
707 
709  const unsigned int symmetric,
710  const unsigned int n_levels,
711  const double threshold,
712  const double filter,
713  const bool output_details)
715  , n_levels(n_levels)
716  , threshold(threshold)
717  , filter(filter)
718  , output_details(output_details)
719  {}
720 
721 
722 
724  : PreconditionBase()
725  {}
726 
727 
728 
730  const MatrixBase &matrix,
731  const AdditionalData &additional_data)
732  : PreconditionBase(matrix.get_mpi_communicator())
733  {
735  }
736 
737 
738  void
740  const AdditionalData &additional_data_)
741  {
742  clear();
743 
744  additional_data = additional_data_;
745 
746 # ifdef DEAL_II_PETSC_WITH_HYPRE
747  create_pc_with_mat(matrix_);
748 
749  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
750  AssertThrow(ierr == 0, ExcPETScError(ierr));
751 
752  ierr = PCHYPRESetType(pc, "parasails");
753  AssertThrow(ierr == 0, ExcPETScError(ierr));
754 
756  {
757  set_option_value("-pc_hypre_parasails_logging", "1");
758  }
759 
762  ExcMessage(
763  "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
764 
765  std::stringstream ssStream;
766 
767  switch (additional_data.symmetric)
768  {
769  case 0:
770  {
771  ssStream << "nonsymmetric";
772  break;
773  }
774 
775  case 1:
776  {
777  ssStream << "SPD";
778  break;
779  }
780 
781  case 2:
782  {
783  ssStream << "nonsymmetric,SPD";
784  break;
785  }
786 
787  default:
788  Assert(
789  false,
790  ExcMessage(
791  "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
792  }
793 
794  set_option_value("-pc_hypre_parasails_sym", ssStream.str());
795 
796  set_option_value("-pc_hypre_parasails_nlevels",
798 
799  ssStream.str(""); // empty the stringstream
800  ssStream << additional_data.threshold;
801  set_option_value("-pc_hypre_parasails_thresh", ssStream.str());
802 
803  ssStream.str(""); // empty the stringstream
804  ssStream << additional_data.filter;
805  set_option_value("-pc_hypre_parasails_filter", ssStream.str());
806 
807  ierr = PCSetFromOptions(pc);
808  AssertThrow(ierr == 0, ExcPETScError(ierr));
809 
810 # else // DEAL_II_PETSC_WITH_HYPRE
811  (void)matrix_;
812  Assert(false,
813  ExcMessage("Your PETSc installation does not include a copy of "
814  "the hypre package necessary for this preconditioner."));
815 # endif
816  }
817 
818 
819  /* ----------------- PreconditionNone ------------------------- */
820 
822  : PreconditionBase()
823  {}
824 
825 
826 
828  const AdditionalData &additional_data)
829  : PreconditionBase(matrix.get_mpi_communicator())
830  {
832  }
833 
834 
835  void
837  const AdditionalData &additional_data_)
838  {
839  clear();
840 
841  additional_data = additional_data_;
842 
843  create_pc_with_mat(matrix_);
844 
845  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCNONE));
846  AssertThrow(ierr == 0, ExcPETScError(ierr));
847 
848  ierr = PCSetFromOptions(pc);
849  AssertThrow(ierr == 0, ExcPETScError(ierr));
850  }
851 
852 
853  /* ----------------- PreconditionLU -------------------- */
854 
856  const double zero_pivot,
857  const double damping)
858  : pivoting(pivoting)
859  , zero_pivot(zero_pivot)
860  , damping(damping)
861  {}
862 
863 
864 
866  : PreconditionBase()
867  {}
868 
869 
870 
872  const AdditionalData &additional_data)
873  : PreconditionBase(matrix.get_mpi_communicator())
874  {
876  }
877 
878 
879  void
881  const AdditionalData &additional_data_)
882  {
883  clear();
884 
885  additional_data = additional_data_;
886 
887  create_pc_with_mat(matrix_);
888 
889  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCLU));
890  AssertThrow(ierr == 0, ExcPETScError(ierr));
891 
892  // set flags as given
893  ierr = PCFactorSetColumnPivot(pc, additional_data.pivoting);
894  AssertThrow(ierr == 0, ExcPETScError(ierr));
895 
896  ierr = PCFactorSetZeroPivot(pc, additional_data.zero_pivot);
897  AssertThrow(ierr == 0, ExcPETScError(ierr));
898 
899  ierr = PCFactorSetShiftAmount(pc, additional_data.damping);
900  AssertThrow(ierr == 0, ExcPETScError(ierr));
901 
902  ierr = PCSetFromOptions(pc);
903  AssertThrow(ierr == 0, ExcPETScError(ierr));
904  }
905 
906  /* ----------------- PreconditionBDDC -------------------- */
907 
908  template <int dim>
910  const bool use_vertices,
911  const bool use_edges,
912  const bool use_faces,
913  const bool symmetric,
914  const std::vector<Point<dim>> coords)
915  : use_vertices(use_vertices)
916  , use_edges(use_edges)
917  , use_faces(use_faces)
919  , coords(coords)
920  {}
921 
922 
923 
924  template <int dim>
926  : PreconditionBase()
927  {}
928 
929 
930 
931  template <int dim>
933  const MPI_Comm comm,
934  const AdditionalData &additional_data_)
936  {
937  additional_data = additional_data_;
938 
939  PetscErrorCode ierr = PCCreate(comm, &pc);
940  AssertThrow(ierr == 0, ExcPETScError(ierr));
941 
942  initialize();
943  }
944 
945 
946 
947  template <int dim>
949  const AdditionalData &additional_data)
950  : PreconditionBase(matrix.get_mpi_communicator())
951  {
953  }
954 
955 
956 
957  template <int dim>
958  void
960  {
961 # if DEAL_II_PETSC_VERSION_GTE(3, 10, 0)
962  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCBDDC));
963  AssertThrow(ierr == 0, ExcPETScError(ierr));
964 
965  // The matrix must be of IS type. We check for this to avoid the PETSc error
966  // in order to suggest the correct matrix reinit method.
967  {
968  MatType current_type;
969  Mat A, P;
970  PetscBool flg;
971 
972  ierr = PCGetOperators(pc, &A, &P);
973  AssertThrow(ierr == 0, ExcPETScError(ierr));
974  ierr = PCGetUseAmat(pc, &flg);
975  AssertThrow(ierr == 0, ExcPETScError(ierr));
976 
977  ierr = MatGetType(flg ? A : P, &current_type);
978  AssertThrow(ierr == 0, ExcPETScError(ierr));
979  AssertThrow(
980  strcmp(current_type, MATIS) == 0,
981  ExcMessage(
982  "Matrix must be of IS type. For this, the variant of reinit that includes the active dofs must be used."));
983  }
984 
985 
986  std::stringstream ssStream;
987 
988  if (additional_data.use_vertices)
989  set_option_value("-pc_bddc_use_vertices", "true");
990  else
991  set_option_value("-pc_bddc_use_vertices", "false");
992  if (additional_data.use_edges)
993  set_option_value("-pc_bddc_use_edges", "true");
994  else
995  set_option_value("-pc_bddc_use_edges", "false");
996  if (additional_data.use_faces)
997  set_option_value("-pc_bddc_use_faces", "true");
998  else
999  set_option_value("-pc_bddc_use_faces", "false");
1000  if (additional_data.symmetric)
1001  set_option_value("-pc_bddc_symmetric", "true");
1002  else
1003  set_option_value("-pc_bddc_symmetric", "false");
1004  if (additional_data.coords.size() > 0)
1005  {
1006  set_option_value("-pc_bddc_corner_selection", "true");
1007  // Convert coords vector to PETSc data array
1008  std::vector<PetscReal> coords_petsc(additional_data.coords.size() *
1009  dim);
1010  for (unsigned int i = 0, j = 0; i < additional_data.coords.size(); ++i)
1011  {
1012  for (j = 0; j < dim; ++j)
1013  coords_petsc[dim * i + j] = additional_data.coords[i][j];
1014  }
1015 
1016  ierr = PCSetCoordinates(pc,
1017  dim,
1018  additional_data.coords.size(),
1019  coords_petsc.data());
1020  AssertThrow(ierr == 0, ExcPETScError(ierr));
1021  }
1022  else
1023  {
1024  set_option_value("-pc_bddc_corner_selection", "false");
1025  ierr = PCSetCoordinates(pc, 0, 0, nullptr);
1026  AssertThrow(ierr == 0, ExcPETScError(ierr));
1027  }
1028 
1029 
1030  ierr = PCSetFromOptions(pc);
1031  AssertThrow(ierr == 0, ExcPETScError(ierr));
1032 # else
1033  AssertThrow(
1034  false, ExcMessage("BDDC preconditioner requires PETSc 3.10.0 or newer"));
1035 # endif
1036  }
1037 
1038 
1039 
1040  template <int dim>
1041  void
1043  const AdditionalData &additional_data_)
1044  {
1045  clear();
1046 
1047  additional_data = additional_data_;
1048 
1049  create_pc_with_mat(matrix_);
1050  initialize();
1051  }
1052 
1053  /* ----------------- PreconditionShell -------------------- */
1054 
1056  {
1057  initialize(matrix);
1058  }
1059 
1061  {
1062  initialize(comm);
1063  }
1064 
1065  void
1067  {
1068  PetscErrorCode ierr;
1069  if (pc)
1070  {
1071  ierr = PCDestroy(&pc);
1072  AssertThrow(ierr == 0, ExcPETScError(ierr));
1073  }
1075 
1076  ierr = PCSetType(pc, PCSHELL);
1077  AssertThrow(ierr == 0, ExcPETScError(ierr));
1078  ierr = PCShellSetContext(pc, static_cast<void *>(this));
1079  AssertThrow(ierr == 0, ExcPETScError(ierr));
1080  ierr = PCShellSetSetUp(pc, PreconditionShell::pcsetup);
1081  AssertThrow(ierr == 0, ExcPETScError(ierr));
1082  ierr = PCShellSetApply(pc, PreconditionShell::pcapply);
1083  AssertThrow(ierr == 0, ExcPETScError(ierr));
1084  ierr = PCShellSetApplyTranspose(pc, PreconditionShell::pcapply_transpose);
1085  AssertThrow(ierr == 0, ExcPETScError(ierr));
1086  ierr = PCShellSetName(pc, "deal.II user solve");
1087  AssertThrow(ierr == 0, ExcPETScError(ierr));
1088  }
1089 
1090  void
1092  {
1093  initialize(matrix.get_mpi_communicator());
1094  PetscErrorCode ierr;
1095  ierr = PCSetOperators(pc, matrix, matrix);
1096  AssertThrow(ierr == 0, ExcPETScError(ierr));
1097  }
1098 
1099 # ifndef PetscCall
1100 # define PetscCall(code) \
1101  do \
1102  { \
1103  PetscErrorCode ierr = (code); \
1104  CHKERRQ(ierr); \
1105  } \
1106  while (false)
1107 # endif
1108 
1109  PetscErrorCode
1111  {
1112  PetscFunctionBeginUser;
1113  // Failed reason is not reset uniformly within the
1114  // interface code of PCSetUp in PETSc.
1115  // We handle it here.
1116  PetscCall(pc_set_failed_reason(ppc, PC_NOERROR));
1117  PetscFunctionReturn(PETSC_SUCCESS);
1118  }
1119 
1120  PetscErrorCode
1121  PreconditionShell::pcapply(PC ppc, Vec x, Vec y)
1122  {
1123  void *ctx;
1124 
1125  PetscFunctionBeginUser;
1126  PetscCall(PCShellGetContext(ppc, &ctx));
1127 
1128  auto *user = static_cast<PreconditionShell *>(ctx);
1129  if (!user->vmult)
1130  SETERRQ(
1131  PetscObjectComm((PetscObject)ppc),
1132  PETSC_ERR_LIB,
1133  "Failure in ::PETScWrappers::PreconditionShell::pcapply. Missing std::function vmult");
1134 
1135  VectorBase src(x);
1136  VectorBase dst(y);
1137  const int lineno = __LINE__;
1138  try
1139  {
1140  user->vmult(dst, src);
1141  }
1142  catch (const RecoverableUserCallbackError &)
1143  {
1144  PetscCall(pc_set_failed_reason(ppc, PC_SUBPC_ERROR));
1145  }
1146  catch (...)
1147  {
1148  return PetscError(
1149  PetscObjectComm((PetscObject)ppc),
1150  lineno + 3,
1151  "vmult",
1152  __FILE__,
1153  PETSC_ERR_LIB,
1154  PETSC_ERROR_INITIAL,
1155  "Failure in pcapply from ::PETScWrappers::NonlinearSolver");
1156  }
1158  PetscFunctionReturn(PETSC_SUCCESS);
1159  }
1160 
1161  PetscErrorCode
1163  {
1164  void *ctx;
1165 
1166  PetscFunctionBeginUser;
1167  PetscCall(PCShellGetContext(ppc, &ctx));
1168 
1169  auto *user = static_cast<PreconditionShell *>(ctx);
1170  if (!user->vmultT)
1171  SETERRQ(
1172  PetscObjectComm((PetscObject)ppc),
1173  PETSC_ERR_LIB,
1174  "Failure in ::PETScWrappers::PreconditionShell::pcapply_transpose. Missing std::function vmultT");
1175 
1176  VectorBase src(x);
1177  VectorBase dst(y);
1178  const int lineno = __LINE__;
1179  try
1180  {
1181  user->vmultT(dst, src);
1182  }
1183  catch (const RecoverableUserCallbackError &)
1184  {
1185  PetscCall(pc_set_failed_reason(ppc, PC_SUBPC_ERROR));
1186  }
1187  catch (...)
1188  {
1189  return PetscError(
1190  PetscObjectComm((PetscObject)ppc),
1191  lineno + 3,
1192  "vmultT",
1193  __FILE__,
1194  PETSC_ERR_LIB,
1195  PETSC_ERROR_INITIAL,
1196  "Failure in pcapply_transpose from ::PETScWrappers::NonlinearSolver");
1197  }
1199  PetscFunctionReturn(PETSC_SUCCESS);
1200  }
1201 
1202 
1203 } // namespace PETScWrappers
1204 
1205 template class PETScWrappers::PreconditionBDDC<2>;
1206 template class PETScWrappers::PreconditionBDDC<3>;
1207 
1209 
1210 #endif // DEAL_II_WITH_PETSC
void Tvmult(VectorBase &dst, const VectorBase &src) const
void create_pc_with_comm(const MPI_Comm)
void vmult(VectorBase &dst, const VectorBase &src) const
void create_pc_with_mat(const MatrixBase &)
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())
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())
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
static PetscErrorCode pcapply_transpose(PC pc, Vec src, Vec dst)
static PetscErrorCode pcsetup(PC pc)
void initialize(const MPI_Comm comm)
static PetscErrorCode pcapply(PC pc, Vec src, Vec dst)
Definition: point.h:112
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & RecoverableUserCallbackError()
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
@ matrix
Contents is actually a matrix.
static const char A
@ symmetric
Matrix is symmetric.
void set_option_value(const std::string &name, const std::string &value)
PetscErrorCode pc_set_failed_reason(PC pc, PCFailedReason reason)
void petsc_increment_state_counter(Vec v)
std::string to_string(const T &t)
Definition: patterns.h:2391
#define PetscCall(code)
AdditionalData(const bool use_vertices=true, const bool use_edges=false, const bool use_faces=false, const bool symmetric=false, const std::vector< Point< dim >> coords={})
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)
AdditionalData(const double pivoting=1.e-6, const double zero_pivot=1.e-12, const double damping=0.0)
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)
const MPI_Comm comm