Reference documentation for deal.II version GIT 0980a66d4b 2023-03-23 20:20:03+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 - 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  {
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  const MPI_Comm &
95  {
96  static MPI_Comm comm = PETSC_COMM_SELF;
97  MPI_Comm pcomm = PetscObjectComm(reinterpret_cast<PetscObject>(pc));
98  if (pcomm != MPI_COMM_NULL)
99  comm = pcomm;
100  return comm;
101  }
102 
103  void
105  {
106  // only allow the creation of the
107  // preconditioner once
109 
110  MPI_Comm comm;
111  PetscErrorCode ierr = PetscObjectGetComm(
112  reinterpret_cast<PetscObject>(static_cast<const Mat &>(matrix)), &comm);
113  AssertThrow(ierr == 0, ExcPETScError(ierr));
114 
116 
117 # if DEAL_II_PETSC_VERSION_LT(3, 5, 0)
118  ierr = PCSetOperators(pc, matrix, matrix, SAME_PRECONDITIONER);
119 # else
120  ierr = PCSetOperators(pc, matrix, matrix);
121 # endif
122  AssertThrow(ierr == 0, ExcPETScError(ierr));
123  }
124 
125  void
127  {
128  clear();
129  PetscErrorCode ierr = PCCreate(comm, &pc);
130  AssertThrow(ierr == 0, ExcPETScError(ierr));
131  }
132 
133  const PC &
135  {
136  return pc;
137  }
138 
139 
140  /* ----------------- PreconditionJacobi -------------------- */
141 
143  : PreconditionBase()
144  {}
145 
146 
147 
149  const AdditionalData &additional_data_)
151  {
152  additional_data = additional_data_;
153 
154  PetscErrorCode ierr = PCCreate(comm, &pc);
155  AssertThrow(ierr == 0, ExcPETScError(ierr));
156 
157  initialize();
158  }
159 
160 
161 
163  const AdditionalData &additional_data)
164  : PreconditionBase(matrix.get_mpi_communicator())
165  {
167  }
168 
169 
170 
171  void
173  {
175 
176  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCJACOBI));
177  AssertThrow(ierr == 0, ExcPETScError(ierr));
178 
179  ierr = PCSetFromOptions(pc);
180  AssertThrow(ierr == 0, ExcPETScError(ierr));
181  }
182 
183 
184 
185  void
187  const AdditionalData &additional_data_)
188  {
189  clear();
190 
191  additional_data = additional_data_;
192 
193  create_pc_with_mat(matrix_);
194  initialize();
195  }
196 
197 
198  /* ----------------- PreconditionBlockJacobi -------------------- */
199 
201  : PreconditionBase()
202  {}
203 
205  const MPI_Comm & comm,
206  const AdditionalData &additional_data_)
208  {
209  additional_data = additional_data_;
210 
211  PetscErrorCode ierr = PCCreate(comm, &pc);
212  AssertThrow(ierr == 0, ExcPETScError(ierr));
213 
214  initialize();
215  }
216 
217 
218 
220  const MatrixBase & matrix,
221  const AdditionalData &additional_data)
222  : PreconditionBase(matrix.get_mpi_communicator())
223  {
225  }
226 
227 
228 
229  void
231  {
232  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCBJACOBI));
233  AssertThrow(ierr == 0, ExcPETScError(ierr));
234 
235  ierr = PCSetFromOptions(pc);
236  AssertThrow(ierr == 0, ExcPETScError(ierr));
237  }
238 
239 
240 
241  void
243  const AdditionalData &additional_data_)
244  {
245  clear();
246 
247  additional_data = additional_data_;
248 
249  create_pc_with_mat(matrix_);
250  initialize();
251  }
252 
253 
254  /* ----------------- PreconditionSOR -------------------- */
255 
257  : PreconditionBase()
258  {}
259 
260 
261 
263  : omega(omega)
264  {}
265 
266 
267 
271  {
273  }
274 
275 
276  void
278  const AdditionalData &additional_data_)
279  {
280  clear();
281 
282  additional_data = additional_data_;
283 
284  create_pc_with_mat(matrix_);
285 
286  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
287  AssertThrow(ierr == 0, ExcPETScError(ierr));
288 
289  // then set flags as given
290  ierr = PCSORSetOmega(pc, additional_data.omega);
291  AssertThrow(ierr == 0, ExcPETScError(ierr));
292 
293  ierr = PCSetFromOptions(pc);
294  AssertThrow(ierr == 0, ExcPETScError(ierr));
295  }
296 
297 
298  /* ----------------- PreconditionSSOR -------------------- */
299 
301  : PreconditionBase()
302  {}
303 
304 
305 
307  : omega(omega)
308  {}
309 
310 
311 
315  {
317  }
318 
319 
320  void
322  const AdditionalData &additional_data_)
323  {
324  clear();
325 
326  additional_data = additional_data_;
327 
328  create_pc_with_mat(matrix_);
329 
330  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
331  AssertThrow(ierr == 0, ExcPETScError(ierr));
332 
333  // then set flags as given
334  ierr = PCSORSetOmega(pc, additional_data.omega);
335  AssertThrow(ierr == 0, ExcPETScError(ierr));
336 
337  // convert SOR to SSOR
338  ierr = PCSORSetSymmetric(pc, SOR_SYMMETRIC_SWEEP);
339  AssertThrow(ierr == 0, ExcPETScError(ierr));
340 
341  ierr = PCSetFromOptions(pc);
342  AssertThrow(ierr == 0, ExcPETScError(ierr));
343  }
344 
345 
346  /* ----------------- PreconditionICC -------------------- */
347 
349  : PreconditionBase()
350  {}
351 
352 
353 
355  : levels(levels)
356  {}
357 
358 
359 
363  {
365  }
366 
367 
368  void
370  const AdditionalData &additional_data_)
371  {
372  clear();
373 
374  additional_data = additional_data_;
375 
376  create_pc_with_mat(matrix_);
377 
378  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCICC));
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 
389 
390  /* ----------------- PreconditionILU -------------------- */
391 
393  : PreconditionBase()
394  {}
395 
396 
397 
399  : levels(levels)
400  {}
401 
402 
403 
407  {
409  }
410 
411 
412  void
414  const AdditionalData &additional_data_)
415  {
416  clear();
417 
418  additional_data = additional_data_;
419 
420  create_pc_with_mat(matrix_);
421 
422  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCILU));
423  AssertThrow(ierr == 0, ExcPETScError(ierr));
424 
425  // then set flags
426  ierr = PCFactorSetLevels(pc, additional_data.levels);
427  AssertThrow(ierr == 0, ExcPETScError(ierr));
428 
429  ierr = PCSetFromOptions(pc);
430  AssertThrow(ierr == 0, ExcPETScError(ierr));
431  }
432 
433 
434  /* ----------------- PreconditionBoomerAMG -------------------- */
435 
437  const bool symmetric_operator,
438  const double strong_threshold,
439  const double max_row_sum,
440  const unsigned int aggressive_coarsening_num_levels,
441  const bool output_details,
442  const RelaxationType relaxation_type_up,
443  const RelaxationType relaxation_type_down,
444  const RelaxationType relaxation_type_coarse,
445  const unsigned int n_sweeps_coarse,
446  const double tol,
447  const unsigned int max_iter,
448  const bool w_cycle)
449  : symmetric_operator(symmetric_operator)
450  , strong_threshold(strong_threshold)
451  , max_row_sum(max_row_sum)
452  , aggressive_coarsening_num_levels(aggressive_coarsening_num_levels)
453  , output_details(output_details)
454  , relaxation_type_up(relaxation_type_up)
455  , relaxation_type_down(relaxation_type_down)
456  , relaxation_type_coarse(relaxation_type_coarse)
457  , n_sweeps_coarse(n_sweeps_coarse)
458  , tol(tol)
459  , max_iter(max_iter)
460  , w_cycle(w_cycle)
461  {}
462 
463 
464 
465 # ifdef DEAL_II_PETSC_WITH_HYPRE
466  namespace
467  {
472  std::string
473  to_string(
475  {
476  std::string string_type;
477 
478  switch (relaxation_type)
479  {
481  string_type = "Jacobi";
482  break;
485  string_type = "sequential-Gauss-Seidel";
486  break;
489  string_type = "seqboundary-Gauss-Seidel";
490  break;
492  string_type = "SOR/Jacobi";
493  break;
496  string_type = "backward-SOR/Jacobi";
497  break;
500  string_type = "symmetric-SOR/Jacobi";
501  break;
504  string_type = " l1scaled-SOR/Jacobi";
505  break;
508  string_type = "Gaussian-elimination";
509  break;
512  string_type = "l1-Gauss-Seidel";
513  break;
516  string_type = "backward-l1-Gauss-Seidel";
517  break;
519  string_type = "CG";
520  break;
522  string_type = "Chebyshev";
523  break;
525  string_type = "FCF-Jacobi";
526  break;
529  string_type = "l1scaled-Jacobi";
530  break;
532  string_type = "None";
533  break;
534  default:
535  Assert(false, ExcNotImplemented());
536  }
537  return string_type;
538  }
539  } // namespace
540 # endif
541 
542 
543 
545  : PreconditionBase()
546  {}
547 
548 
549 
551  const MPI_Comm & comm,
552  const AdditionalData &additional_data_)
554  {
555  additional_data = additional_data_;
556 
557  PetscErrorCode ierr = PCCreate(comm, &pc);
558  AssertThrow(ierr == 0, ExcPETScError(ierr));
559 
560 # ifdef DEAL_II_PETSC_WITH_HYPRE
561  initialize();
562 # else // DEAL_II_PETSC_WITH_HYPRE
563  (void)pc;
564  Assert(false,
565  ExcMessage("Your PETSc installation does not include a copy of "
566  "the hypre package necessary for this preconditioner."));
567 # endif
568  }
569 
570 
571 
573  const MatrixBase & matrix,
574  const AdditionalData &additional_data)
575  : PreconditionBase(matrix.get_mpi_communicator())
576  {
578  }
579 
580 
581 
582  void
584  {
585 # ifdef DEAL_II_PETSC_WITH_HYPRE
586  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
587  AssertThrow(ierr == 0, ExcPETScError(ierr));
588 
589  ierr = PCHYPRESetType(pc, "boomeramg");
590  AssertThrow(ierr == 0, ExcPETScError(ierr));
591 
593  {
594  set_option_value("-pc_hypre_boomeramg_print_statistics", "1");
595  }
596 
597  set_option_value("-pc_hypre_boomeramg_agg_nl",
600 
601  set_option_value("-pc_hypre_boomeramg_max_row_sum",
603 
604  set_option_value("-pc_hypre_boomeramg_strong_threshold",
606 
607  // change to symmetric SOR/Jacobi when using a symmetric operator for
608  // backward compatibility
612  {
615  }
616 
620  {
623  }
624 
628  {
631  }
632 
633  auto relaxation_type_is_symmetric =
634  [](AdditionalData::RelaxationType relaxation_type) {
635  return relaxation_type == AdditionalData::RelaxationType::Jacobi ||
636  relaxation_type ==
638  relaxation_type ==
640  relaxation_type == AdditionalData::RelaxationType::None ||
641  relaxation_type ==
643  relaxation_type == AdditionalData::RelaxationType::CG ||
645  };
646 
648  !relaxation_type_is_symmetric(additional_data.relaxation_type_up))
649  Assert(false,
650  ExcMessage("Use a symmetric smoother for relaxation_type_up"));
651 
653  !relaxation_type_is_symmetric(additional_data.relaxation_type_down))
654  Assert(false,
655  ExcMessage("Use a symmetric smoother for relaxation_type_down"));
656 
658  !relaxation_type_is_symmetric(additional_data.relaxation_type_coarse))
659  Assert(false,
660  ExcMessage("Use a symmetric smoother for relaxation_type_coarse"));
661 
662  set_option_value("-pc_hypre_boomeramg_relax_type_up",
664  set_option_value("-pc_hypre_boomeramg_relax_type_down",
666  set_option_value("-pc_hypre_boomeramg_relax_type_coarse",
668  set_option_value("-pc_hypre_boomeramg_grid_sweeps_coarse",
670 
671  set_option_value("-pc_hypre_boomeramg_tol",
673  set_option_value("-pc_hypre_boomeramg_max_iter",
675 
677  {
678  set_option_value("-pc_hypre_boomeramg_cycle_type", "W");
679  }
680 
681  ierr = PCSetFromOptions(pc);
682  AssertThrow(ierr == 0, ExcPETScError(ierr));
683 # else
684  Assert(false,
685  ExcMessage("Your PETSc installation does not include a copy of "
686  "the hypre package necessary for this preconditioner."));
687 # endif
688  }
689 
690 
691 
692  void
694  const AdditionalData &additional_data_)
695  {
696 # ifdef DEAL_II_PETSC_WITH_HYPRE
697  clear();
698 
699  additional_data = additional_data_;
700 
701  create_pc_with_mat(matrix_);
702  initialize();
703 
704 # else // DEAL_II_PETSC_WITH_HYPRE
705  (void)matrix_;
706  (void)additional_data_;
707  Assert(false,
708  ExcMessage("Your PETSc installation does not include a copy of "
709  "the hypre package necessary for this preconditioner."));
710 # endif
711  }
712 
713 
714  /* ----------------- PreconditionParaSails -------------------- */
715 
717  const unsigned int symmetric,
718  const unsigned int n_levels,
719  const double threshold,
720  const double filter,
721  const bool output_details)
723  , n_levels(n_levels)
724  , threshold(threshold)
725  , filter(filter)
726  , output_details(output_details)
727  {}
728 
729 
730 
732  : PreconditionBase()
733  {}
734 
735 
736 
738  const MatrixBase & matrix,
739  const AdditionalData &additional_data)
740  : PreconditionBase(matrix.get_mpi_communicator())
741  {
743  }
744 
745 
746  void
748  const AdditionalData &additional_data_)
749  {
750  clear();
751 
752  additional_data = additional_data_;
753 
754 # ifdef DEAL_II_PETSC_WITH_HYPRE
755  create_pc_with_mat(matrix_);
756 
757  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
758  AssertThrow(ierr == 0, ExcPETScError(ierr));
759 
760  ierr = PCHYPRESetType(pc, "parasails");
761  AssertThrow(ierr == 0, ExcPETScError(ierr));
762 
764  {
765  set_option_value("-pc_hypre_parasails_logging", "1");
766  }
767 
770  ExcMessage(
771  "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
772 
773  std::stringstream ssStream;
774 
775  switch (additional_data.symmetric)
776  {
777  case 0:
778  {
779  ssStream << "nonsymmetric";
780  break;
781  }
782 
783  case 1:
784  {
785  ssStream << "SPD";
786  break;
787  }
788 
789  case 2:
790  {
791  ssStream << "nonsymmetric,SPD";
792  break;
793  }
794 
795  default:
796  Assert(
797  false,
798  ExcMessage(
799  "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
800  }
801 
802  set_option_value("-pc_hypre_parasails_sym", ssStream.str());
803 
804  set_option_value("-pc_hypre_parasails_nlevels",
806 
807  ssStream.str(""); // empty the stringstream
808  ssStream << additional_data.threshold;
809  set_option_value("-pc_hypre_parasails_thresh", ssStream.str());
810 
811  ssStream.str(""); // empty the stringstream
812  ssStream << additional_data.filter;
813  set_option_value("-pc_hypre_parasails_filter", ssStream.str());
814 
815  ierr = PCSetFromOptions(pc);
816  AssertThrow(ierr == 0, ExcPETScError(ierr));
817 
818 # else // DEAL_II_PETSC_WITH_HYPRE
819  (void)matrix_;
820  Assert(false,
821  ExcMessage("Your PETSc installation does not include a copy of "
822  "the hypre package necessary for this preconditioner."));
823 # endif
824  }
825 
826 
827  /* ----------------- PreconditionNone ------------------------- */
828 
830  : PreconditionBase()
831  {}
832 
833 
834 
836  const AdditionalData &additional_data)
837  : PreconditionBase(matrix.get_mpi_communicator())
838  {
840  }
841 
842 
843  void
845  const AdditionalData &additional_data_)
846  {
847  clear();
848 
849  additional_data = additional_data_;
850 
851  create_pc_with_mat(matrix_);
852 
853  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCNONE));
854  AssertThrow(ierr == 0, ExcPETScError(ierr));
855 
856  ierr = PCSetFromOptions(pc);
857  AssertThrow(ierr == 0, ExcPETScError(ierr));
858  }
859 
860 
861  /* ----------------- PreconditionLU -------------------- */
862 
864  const double zero_pivot,
865  const double damping)
866  : pivoting(pivoting)
867  , zero_pivot(zero_pivot)
868  , damping(damping)
869  {}
870 
871 
872 
874  : PreconditionBase()
875  {}
876 
877 
878 
880  const AdditionalData &additional_data)
881  : PreconditionBase(matrix.get_mpi_communicator())
882  {
884  }
885 
886 
887  void
889  const AdditionalData &additional_data_)
890  {
891  clear();
892 
893  additional_data = additional_data_;
894 
895  create_pc_with_mat(matrix_);
896 
897  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCLU));
898  AssertThrow(ierr == 0, ExcPETScError(ierr));
899 
900  // set flags as given
901  ierr = PCFactorSetColumnPivot(pc, additional_data.pivoting);
902  AssertThrow(ierr == 0, ExcPETScError(ierr));
903 
904  ierr = PCFactorSetZeroPivot(pc, additional_data.zero_pivot);
905  AssertThrow(ierr == 0, ExcPETScError(ierr));
906 
907  ierr = PCFactorSetShiftAmount(pc, additional_data.damping);
908  AssertThrow(ierr == 0, ExcPETScError(ierr));
909 
910  ierr = PCSetFromOptions(pc);
911  AssertThrow(ierr == 0, ExcPETScError(ierr));
912  }
913 
914  /* ----------------- PreconditionBDDC -------------------- */
915 
916  template <int dim>
918  const bool use_vertices,
919  const bool use_edges,
920  const bool use_faces,
921  const bool symmetric,
922  const std::vector<Point<dim>> coords)
923  : use_vertices(use_vertices)
924  , use_edges(use_edges)
925  , use_faces(use_faces)
927  , coords(coords)
928  {}
929 
930 
931 
932  template <int dim>
934  : PreconditionBase()
935  {}
936 
937 
938 
939  template <int dim>
941  const MPI_Comm comm,
942  const AdditionalData &additional_data_)
944  {
945  additional_data = additional_data_;
946 
947  PetscErrorCode ierr = PCCreate(comm, &pc);
948  AssertThrow(ierr == 0, ExcPETScError(ierr));
949 
950  initialize();
951  }
952 
953 
954 
955  template <int dim>
957  const AdditionalData &additional_data)
958  : PreconditionBase(matrix.get_mpi_communicator())
959  {
961  }
962 
963 
964 
965  template <int dim>
966  void
968  {
969 # if DEAL_II_PETSC_VERSION_GTE(3, 10, 0)
970  PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCBDDC));
971  AssertThrow(ierr == 0, ExcPETScError(ierr));
972 
973  // The matrix must be of IS type. We check for this to avoid the PETSc error
974  // in order to suggest the correct matrix reinit method.
975  {
976  MatType current_type;
977  Mat A, P;
978  PetscBool flg;
979 
980  ierr = PCGetOperators(pc, &A, &P);
981  AssertThrow(ierr == 0, ExcPETScError(ierr));
982  ierr = PCGetUseAmat(pc, &flg);
983  AssertThrow(ierr == 0, ExcPETScError(ierr));
984 
985  ierr = MatGetType(flg ? A : P, &current_type);
986  AssertThrow(ierr == 0, ExcPETScError(ierr));
987  AssertThrow(
988  strcmp(current_type, MATIS) == 0,
989  ExcMessage(
990  "Matrix must be of IS type. For this, the variant of reinit that includes the active dofs must be used."));
991  }
992 
993 
994  std::stringstream ssStream;
995 
996  if (additional_data.use_vertices)
997  set_option_value("-pc_bddc_use_vertices", "true");
998  else
999  set_option_value("-pc_bddc_use_vertices", "false");
1000  if (additional_data.use_edges)
1001  set_option_value("-pc_bddc_use_edges", "true");
1002  else
1003  set_option_value("-pc_bddc_use_edges", "false");
1004  if (additional_data.use_faces)
1005  set_option_value("-pc_bddc_use_faces", "true");
1006  else
1007  set_option_value("-pc_bddc_use_faces", "false");
1008  if (additional_data.symmetric)
1009  set_option_value("-pc_bddc_symmetric", "true");
1010  else
1011  set_option_value("-pc_bddc_symmetric", "false");
1012  if (additional_data.coords.size() > 0)
1013  {
1014  set_option_value("-pc_bddc_corner_selection", "true");
1015  // Convert coords vector to PETSc data array
1016  std::vector<PetscReal> coords_petsc(additional_data.coords.size() *
1017  dim);
1018  for (unsigned int i = 0, j = 0; i < additional_data.coords.size(); ++i)
1019  {
1020  for (j = 0; j < dim; ++j)
1021  coords_petsc[dim * i + j] = additional_data.coords[i][j];
1022  }
1023 
1024  ierr = PCSetCoordinates(pc,
1025  dim,
1026  additional_data.coords.size(),
1027  coords_petsc.data());
1028  AssertThrow(ierr == 0, ExcPETScError(ierr));
1029  }
1030  else
1031  {
1032  set_option_value("-pc_bddc_corner_selection", "false");
1033  ierr = PCSetCoordinates(pc, 0, 0, nullptr);
1034  AssertThrow(ierr == 0, ExcPETScError(ierr));
1035  }
1036 
1037 
1038  ierr = PCSetFromOptions(pc);
1039  AssertThrow(ierr == 0, ExcPETScError(ierr));
1040 # else
1041  AssertThrow(
1042  false, ExcMessage("BDDC preconditioner requires PETSc 3.10.0 or newer"));
1043 # endif
1044  }
1045 
1046 
1047 
1048  template <int dim>
1049  void
1051  const AdditionalData &additional_data_)
1052  {
1053  clear();
1054 
1055  additional_data = additional_data_;
1056 
1057  create_pc_with_mat(matrix_);
1058  initialize();
1059  }
1060 
1061  /* ----------------- PreconditionShell -------------------- */
1062 
1064  {
1065  initialize(matrix);
1066  }
1067 
1069  {
1070  initialize(comm);
1071  }
1072 
1073  void
1075  {
1076  PetscErrorCode ierr;
1077  if (pc)
1078  {
1079  ierr = PCDestroy(&pc);
1080  AssertThrow(ierr == 0, ExcPETScError(ierr));
1081  }
1083 
1084  ierr = PCSetType(pc, PCSHELL);
1085  AssertThrow(ierr == 0, ExcPETScError(ierr));
1086  ierr = PCShellSetContext(pc, static_cast<void *>(this));
1087  AssertThrow(ierr == 0, ExcPETScError(ierr));
1088  ierr = PCShellSetApply(pc, PreconditionShell::pcapply);
1089  AssertThrow(ierr == 0, ExcPETScError(ierr));
1090  ierr = PCShellSetApplyTranspose(pc, PreconditionShell::pcapply_transpose);
1091  AssertThrow(ierr == 0, ExcPETScError(ierr));
1092  ierr = PCShellSetName(pc, "deal.II user solve");
1093  AssertThrow(ierr == 0, ExcPETScError(ierr));
1094  }
1095 
1096  void
1098  {
1099  initialize(matrix.get_mpi_communicator());
1100  PetscErrorCode ierr;
1101 # if DEAL_II_PETSC_VERSION_LT(3, 5, 0)
1102  ierr = PCSetOperators(pc, matrix, matrix, DIFFERENT_NONZERO_PATTERN);
1103 # else
1104  ierr = PCSetOperators(pc, matrix, matrix);
1105 # endif
1106  AssertThrow(ierr == 0, ExcPETScError(ierr));
1107  }
1108 
1109  int
1110  PreconditionShell::pcapply(PC ppc, Vec x, Vec y)
1111  {
1112  void *ctx;
1113 
1114  PetscFunctionBeginUser;
1115  PetscErrorCode ierr = PCShellGetContext(ppc, &ctx);
1116  AssertThrow(ierr == 0, ExcPETScError(ierr));
1117 
1118  auto user = static_cast<PreconditionShell *>(ctx);
1119  AssertThrow(user->vmult,
1121  "std::function vmult"));
1122 
1123  VectorBase src(x);
1124  VectorBase dst(y);
1125  user->vmult(dst, src);
1126  PetscFunctionReturn(0);
1127  }
1128 
1129  int
1131  {
1132  void *ctx;
1133 
1134  PetscFunctionBeginUser;
1135  PetscErrorCode ierr = PCShellGetContext(ppc, &ctx);
1136  AssertThrow(ierr == 0, ExcPETScError(ierr));
1137 
1138  auto user = static_cast<PreconditionShell *>(ctx);
1139  AssertThrow(user->vmultT,
1141  "std::function vmultT"));
1142 
1143  VectorBase src(x);
1144  VectorBase dst(y);
1145  user->vmult(dst, src);
1146  PetscFunctionReturn(0);
1147  }
1148 
1149 
1150 } // namespace PETScWrappers
1151 
1152 template class PETScWrappers::PreconditionBDDC<2>;
1153 template class PETScWrappers::PreconditionBDDC<3>;
1154 
1156 
1157 #endif // DEAL_II_WITH_PETSC
void Tvmult(VectorBase &dst, const VectorBase &src) const
const MPI_Comm & get_mpi_communicator() const
void vmult(VectorBase &dst, const VectorBase &src) const
void create_pc_with_comm(const MPI_Comm &)
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 int pcapply_transpose(PC pc, Vec src, Vec dst)
static int pcapply(PC pc, Vec src, Vec dst)
void initialize(const MPI_Comm &comm)
Definition: point.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
@ 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)
std::string to_string(const T &t)
Definition: patterns.h:2392
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