deal.II version GIT relicensing-3540-g7552a02177 2025-06-20 13:50:00+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\}}\)
Loading...
Searching...
No Matches
petsc_precondition.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2004 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16
17#ifdef DEAL_II_WITH_PETSC
18
20
26
27# include <petscconf.h>
28
29# include <cmath>
30
32
33namespace PETScWrappers
34{
40
42 : pc(nullptr)
43 {}
44
46 {
47 try
48 {
49 clear();
50 }
51 catch (...)
52 {}
53 }
54
55 void
57 {
58 if (pc)
59 {
60 PetscErrorCode ierr = PCDestroy(&pc);
61 AssertThrow(ierr == 0, ExcPETScError(ierr));
62 }
63 }
64
65 void
67 {
69
70 PetscErrorCode ierr = PCApply(pc, src, dst);
71 AssertThrow(ierr == 0, ExcPETScError(ierr));
72 }
73
74 void
76 {
78
79 PetscErrorCode ierr = PCApplyTranspose(pc, src, dst);
80 AssertThrow(ierr == 0, ExcPETScError(ierr));
81 }
82
83 void
85 {
87
88 PetscErrorCode ierr = PCSetUp(pc);
89 AssertThrow(ierr == 0, ExcPETScError(ierr));
90 }
91
94 {
95 return PetscObjectComm(reinterpret_cast<PetscObject>(pc));
96 }
97
98 void
100 {
101 // only allow the creation of the
102 // preconditioner once
104
106 PetscErrorCode ierr = PetscObjectGetComm(
107 reinterpret_cast<PetscObject>(static_cast<const Mat &>(matrix)), &comm);
108 AssertThrow(ierr == 0, ExcPETScError(ierr));
109
111
112 ierr = PCSetOperators(pc, matrix, matrix);
113 AssertThrow(ierr == 0, ExcPETScError(ierr));
114 }
115
116 void
118 {
119 clear();
120 PetscErrorCode ierr = PCCreate(comm, &pc);
121 AssertThrow(ierr == 0, ExcPETScError(ierr));
122 }
123
124 const PC &
126 {
127 return pc;
128 }
129
130
131 /* ----------------- PreconditionJacobi -------------------- */
132
136
137
138
140 const AdditionalData &additional_data_)
142 {
143 additional_data = additional_data_;
144
145 PetscErrorCode ierr = PCCreate(comm, &pc);
146 AssertThrow(ierr == 0, ExcPETScError(ierr));
147
148 initialize();
149 }
150
151
152
154 const AdditionalData &additional_data)
155 : PreconditionBase(matrix.get_mpi_communicator())
156 {
158 }
159
160
161
162 void
164 {
166
167 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCJACOBI));
168 AssertThrow(ierr == 0, ExcPETScError(ierr));
169
170 ierr = PCSetFromOptions(pc);
171 AssertThrow(ierr == 0, ExcPETScError(ierr));
172 }
173
174
175
176 void
178 const AdditionalData &additional_data_)
179 {
180 clear();
181
182 additional_data = additional_data_;
183
184 create_pc_with_mat(matrix_);
185 initialize();
186 }
187
188
189 /* ----------------- PreconditionBlockJacobi -------------------- */
190
194
196 const MPI_Comm comm,
197 const AdditionalData &additional_data_)
199 {
200 additional_data = additional_data_;
201
202 PetscErrorCode ierr = PCCreate(comm, &pc);
203 AssertThrow(ierr == 0, ExcPETScError(ierr));
204
205 initialize();
206 }
207
208
209
211 const MatrixBase &matrix,
212 const AdditionalData &additional_data)
213 : PreconditionBase(matrix.get_mpi_communicator())
214 {
216 }
217
218
219
220 void
222 {
223 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCBJACOBI));
224 AssertThrow(ierr == 0, ExcPETScError(ierr));
225
226 ierr = PCSetFromOptions(pc);
227 AssertThrow(ierr == 0, ExcPETScError(ierr));
228 }
229
230
231
232 void
234 const AdditionalData &additional_data_)
235 {
236 clear();
237
238 additional_data = additional_data_;
239
240 create_pc_with_mat(matrix_);
241 initialize();
242 }
243
244
245 /* ----------------- PreconditionSOR -------------------- */
246
250
251
252
254 : omega(omega)
255 {}
256
257
258
265
266
267 void
269 const AdditionalData &additional_data_)
270 {
271 clear();
272
273 additional_data = additional_data_;
274
275 create_pc_with_mat(matrix_);
276
277 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
278 AssertThrow(ierr == 0, ExcPETScError(ierr));
279
280 // then set flags as given
281 ierr = PCSORSetOmega(pc, additional_data.omega);
282 AssertThrow(ierr == 0, ExcPETScError(ierr));
283
284 ierr = PCSetFromOptions(pc);
285 AssertThrow(ierr == 0, ExcPETScError(ierr));
286 }
287
288
289 /* ----------------- PreconditionSSOR -------------------- */
290
294
295
296
298 : omega(omega)
299 {}
300
301
302
309
310
311 void
313 const AdditionalData &additional_data_)
314 {
315 clear();
316
317 additional_data = additional_data_;
318
319 create_pc_with_mat(matrix_);
320
321 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
322 AssertThrow(ierr == 0, ExcPETScError(ierr));
323
324 // then set flags as given
325 ierr = PCSORSetOmega(pc, additional_data.omega);
326 AssertThrow(ierr == 0, ExcPETScError(ierr));
327
328 // convert SOR to SSOR
329 ierr = PCSORSetSymmetric(pc, SOR_SYMMETRIC_SWEEP);
330 AssertThrow(ierr == 0, ExcPETScError(ierr));
331
332 ierr = PCSetFromOptions(pc);
333 AssertThrow(ierr == 0, ExcPETScError(ierr));
334 }
335
336
337 /* ----------------- PreconditionICC -------------------- */
338
342
343
344
346 : levels(levels)
347 {}
348
349
350
357
358
359 void
361 const AdditionalData &additional_data_)
362 {
363 clear();
364
365 additional_data = additional_data_;
366
367 create_pc_with_mat(matrix_);
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
380
381 /* ----------------- PreconditionILU -------------------- */
382
386
387
388
390 : levels(levels)
391 {}
392
393
394
401
402
403 void
405 const AdditionalData &additional_data_)
406 {
407 clear();
408
409 additional_data = additional_data_;
410
411 create_pc_with_mat(matrix_);
412
413 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCILU));
414 AssertThrow(ierr == 0, ExcPETScError(ierr));
415
416 // then set flags
417 ierr = PCFactorSetLevels(pc, additional_data.levels);
418 AssertThrow(ierr == 0, ExcPETScError(ierr));
419
420 ierr = PCSetFromOptions(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 const RelaxationType relaxation_type_up,
434 const RelaxationType relaxation_type_down,
435 const RelaxationType relaxation_type_coarse,
436 const unsigned int n_sweeps_coarse,
437 const double tol,
438 const unsigned int max_iter,
439 const bool w_cycle)
440 : symmetric_operator(symmetric_operator)
441 , strong_threshold(strong_threshold)
442 , max_row_sum(max_row_sum)
443 , aggressive_coarsening_num_levels(aggressive_coarsening_num_levels)
444 , output_details(output_details)
445 , relaxation_type_up(relaxation_type_up)
446 , relaxation_type_down(relaxation_type_down)
447 , relaxation_type_coarse(relaxation_type_coarse)
448 , n_sweeps_coarse(n_sweeps_coarse)
449 , tol(tol)
450 , max_iter(max_iter)
451 , w_cycle(w_cycle)
452 {}
453
454
455
456# ifdef DEAL_II_PETSC_WITH_HYPRE
457 namespace
458 {
463 std::string
464 to_string(
466 {
467 std::string string_type;
468
469 switch (relaxation_type)
470 {
472 string_type = "Jacobi";
473 break;
476 string_type = "sequential-Gauss-Seidel";
477 break;
480 string_type = "seqboundary-Gauss-Seidel";
481 break;
483 string_type = "SOR/Jacobi";
484 break;
487 string_type = "backward-SOR/Jacobi";
488 break;
491 string_type = "symmetric-SOR/Jacobi";
492 break;
495 string_type = " l1scaled-SOR/Jacobi";
496 break;
499 string_type = "Gaussian-elimination";
500 break;
503 string_type = "l1-Gauss-Seidel";
504 break;
507 string_type = "backward-l1-Gauss-Seidel";
508 break;
510 string_type = "CG";
511 break;
513 string_type = "Chebyshev";
514 break;
516 string_type = "FCF-Jacobi";
517 break;
520 string_type = "l1scaled-Jacobi";
521 break;
523 string_type = "None";
524 break;
525 default:
527 }
528 return string_type;
529 }
530 } // namespace
531# endif
532
533
534
538
539
540
542 const MPI_Comm comm,
543 const AdditionalData &additional_data_)
545 {
546 additional_data = additional_data_;
547
548 PetscErrorCode ierr = PCCreate(comm, &pc);
549 AssertThrow(ierr == 0, ExcPETScError(ierr));
550
551# ifdef DEAL_II_PETSC_WITH_HYPRE
552 initialize();
553# else // DEAL_II_PETSC_WITH_HYPRE
554 Assert(false,
555 ExcMessage("Your PETSc installation does not include a copy of "
556 "the hypre package necessary for this preconditioner."));
557# endif
558 }
559
560
561
563 const MatrixBase &matrix,
564 const AdditionalData &additional_data)
565 : PreconditionBase(matrix.get_mpi_communicator())
566 {
568 }
569
570
571
572 void
574 {
575# ifdef DEAL_II_PETSC_WITH_HYPRE
576 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
577 AssertThrow(ierr == 0, ExcPETScError(ierr));
578
579 ierr = PCHYPRESetType(pc, "boomeramg");
580 AssertThrow(ierr == 0, ExcPETScError(ierr));
581
583 {
584 set_option_value("-pc_hypre_boomeramg_print_statistics", "1");
585 }
586
587 set_option_value("-pc_hypre_boomeramg_agg_nl",
588 std::to_string(
590
591 set_option_value("-pc_hypre_boomeramg_max_row_sum",
592 std::to_string(additional_data.max_row_sum));
593
594 set_option_value("-pc_hypre_boomeramg_strong_threshold",
595 std::to_string(additional_data.strong_threshold));
596
597 // change to symmetric SOR/Jacobi when using a symmetric operator for
598 // backward compatibility
602 {
605 }
606
610 {
613 }
614
618 {
621 }
622
623 auto relaxation_type_is_symmetric =
624 [](AdditionalData::RelaxationType relaxation_type) {
625 return relaxation_type == AdditionalData::RelaxationType::Jacobi ||
626 relaxation_type ==
628 relaxation_type ==
630 relaxation_type == AdditionalData::RelaxationType::None ||
631 relaxation_type ==
633 relaxation_type == AdditionalData::RelaxationType::CG ||
635 };
636
638 !relaxation_type_is_symmetric(additional_data.relaxation_type_up))
639 Assert(false,
640 ExcMessage("Use a symmetric smoother for relaxation_type_up"));
641
643 !relaxation_type_is_symmetric(additional_data.relaxation_type_down))
644 Assert(false,
645 ExcMessage("Use a symmetric smoother for relaxation_type_down"));
646
648 !relaxation_type_is_symmetric(additional_data.relaxation_type_coarse))
649 Assert(false,
650 ExcMessage("Use a symmetric smoother for relaxation_type_coarse"));
651
652 set_option_value("-pc_hypre_boomeramg_relax_type_up",
654 set_option_value("-pc_hypre_boomeramg_relax_type_down",
656 set_option_value("-pc_hypre_boomeramg_relax_type_coarse",
658 set_option_value("-pc_hypre_boomeramg_grid_sweeps_coarse",
659 std::to_string(additional_data.n_sweeps_coarse));
660
661 set_option_value("-pc_hypre_boomeramg_tol",
662 std::to_string(additional_data.tol));
663 set_option_value("-pc_hypre_boomeramg_max_iter",
664 std::to_string(additional_data.max_iter));
665
667 {
668 set_option_value("-pc_hypre_boomeramg_cycle_type", "W");
669 }
670
671 ierr = PCSetFromOptions(pc);
672 AssertThrow(ierr == 0, ExcPETScError(ierr));
673# else
674 Assert(false,
675 ExcMessage("Your PETSc installation does not include a copy of "
676 "the hypre package necessary for this preconditioner."));
677# endif
678 }
679
680
681
682 void
684 const AdditionalData &additional_data_)
685 {
686# ifdef DEAL_II_PETSC_WITH_HYPRE
687 clear();
688
689 additional_data = additional_data_;
690
691 create_pc_with_mat(matrix_);
692 initialize();
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
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
724
725
726
728 const MatrixBase &matrix,
729 const AdditionalData &additional_data)
730 : PreconditionBase(matrix.get_mpi_communicator())
731 {
733 }
734
735
736 void
738 const AdditionalData &additional_data_)
739 {
740 clear();
741
742 additional_data = additional_data_;
743
744# ifdef DEAL_II_PETSC_WITH_HYPRE
745 create_pc_with_mat(matrix_);
746
747 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
748 AssertThrow(ierr == 0, ExcPETScError(ierr));
749
750 ierr = PCHYPRESetType(pc, "parasails");
751 AssertThrow(ierr == 0, ExcPETScError(ierr));
752
754 {
755 set_option_value("-pc_hypre_parasails_logging", "1");
756 }
757
761 "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
762
763 std::stringstream ssStream;
764
766 {
767 case 0:
768 {
769 ssStream << "nonsymmetric";
770 break;
771 }
772
773 case 1:
774 {
775 ssStream << "SPD";
776 break;
777 }
778
779 case 2:
780 {
781 ssStream << "nonsymmetric,SPD";
782 break;
783 }
784
785 default:
786 Assert(
787 false,
789 "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
790 }
791
792 set_option_value("-pc_hypre_parasails_sym", ssStream.str());
793
794 set_option_value("-pc_hypre_parasails_nlevels",
795 std::to_string(additional_data.n_levels));
796
797 ssStream.str(""); // empty the stringstream
798 ssStream << additional_data.threshold;
799 set_option_value("-pc_hypre_parasails_thresh", ssStream.str());
800
801 ssStream.str(""); // empty the stringstream
802 ssStream << additional_data.filter;
803 set_option_value("-pc_hypre_parasails_filter", ssStream.str());
804
805 ierr = PCSetFromOptions(pc);
806 AssertThrow(ierr == 0, ExcPETScError(ierr));
807
808# else // DEAL_II_PETSC_WITH_HYPRE
809 (void)matrix_;
810 Assert(false,
811 ExcMessage("Your PETSc installation does not include a copy of "
812 "the hypre package necessary for this preconditioner."));
813# endif
814 }
815
816
817 /* ----------------- PreconditionNone ------------------------- */
818
822
823
824
826 const AdditionalData &additional_data)
827 : PreconditionBase(matrix.get_mpi_communicator())
828 {
830 }
831
832
833 void
835 const AdditionalData &additional_data_)
836 {
837 clear();
838
839 additional_data = additional_data_;
840
841 create_pc_with_mat(matrix_);
842
843 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCNONE));
844 AssertThrow(ierr == 0, ExcPETScError(ierr));
845
846 ierr = PCSetFromOptions(pc);
847 AssertThrow(ierr == 0, ExcPETScError(ierr));
848 }
849
850
851 /* ----------------- PreconditionLU -------------------- */
852
854 const double zero_pivot,
855 const double damping)
856 : pivoting(pivoting)
857 , zero_pivot(zero_pivot)
858 , damping(damping)
859 {}
860
861
862
866
867
868
870 const AdditionalData &additional_data)
871 : PreconditionBase(matrix.get_mpi_communicator())
872 {
874 }
875
876
877 void
879 const AdditionalData &additional_data_)
880 {
881 clear();
882
883 additional_data = additional_data_;
884
885 create_pc_with_mat(matrix_);
886
887 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCLU));
888 AssertThrow(ierr == 0, ExcPETScError(ierr));
889
890 // set flags as given
891 ierr = PCFactorSetColumnPivot(pc, additional_data.pivoting);
892 AssertThrow(ierr == 0, ExcPETScError(ierr));
893
894 ierr = PCFactorSetZeroPivot(pc, additional_data.zero_pivot);
895 AssertThrow(ierr == 0, ExcPETScError(ierr));
896
897 ierr = PCFactorSetShiftAmount(pc, additional_data.damping);
898 AssertThrow(ierr == 0, ExcPETScError(ierr));
899
900 ierr = PCSetFromOptions(pc);
901 AssertThrow(ierr == 0, ExcPETScError(ierr));
902 }
903
904 /* ----------------- PreconditionBDDC -------------------- */
905
906 template <int dim>
908 const bool use_vertices,
909 const bool use_edges,
910 const bool use_faces,
911 const bool symmetric,
912 const std::vector<Point<dim>> coords)
913 : use_vertices(use_vertices)
914 , use_edges(use_edges)
915 , use_faces(use_faces)
916 , symmetric(symmetric)
917 , coords(coords)
918 {}
919
920
921
922 template <int dim>
926
927
928
929 template <int dim>
931 const MPI_Comm comm,
932 const AdditionalData &additional_data_)
934 {
935 additional_data = additional_data_;
936
937 PetscErrorCode ierr = PCCreate(comm, &pc);
938 AssertThrow(ierr == 0, ExcPETScError(ierr));
939
940 initialize();
941 }
942
943
944
945 template <int dim>
947 const AdditionalData &additional_data)
948 : PreconditionBase(matrix.get_mpi_communicator())
949 {
951 }
952
953
954
955 template <int dim>
956 void
958 {
959# if DEAL_II_PETSC_VERSION_GTE(3, 10, 0)
960 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCBDDC));
961 AssertThrow(ierr == 0, ExcPETScError(ierr));
962
963 // The matrix must be of IS type. We check for this to avoid the PETSc error
964 // in order to suggest the correct matrix reinit method.
965 {
966 MatType current_type;
967 Mat A, P;
968 PetscBool flg;
969
970 ierr = PCGetOperators(pc, &A, &P);
971 AssertThrow(ierr == 0, ExcPETScError(ierr));
972 ierr = PCGetUseAmat(pc, &flg);
973 AssertThrow(ierr == 0, ExcPETScError(ierr));
974
975 ierr = MatGetType(flg ? A : P, &current_type);
976 AssertThrow(ierr == 0, ExcPETScError(ierr));
978 strcmp(current_type, MATIS) == 0,
980 "Matrix must be of IS type. For this, the variant of reinit that includes the active dofs must be used."));
981 }
982
983
984 std::stringstream ssStream;
985
986 if (additional_data.use_vertices)
987 set_option_value("-pc_bddc_use_vertices", "true");
988 else
989 set_option_value("-pc_bddc_use_vertices", "false");
990 if (additional_data.use_edges)
991 set_option_value("-pc_bddc_use_edges", "true");
992 else
993 set_option_value("-pc_bddc_use_edges", "false");
994 if (additional_data.use_faces)
995 set_option_value("-pc_bddc_use_faces", "true");
996 else
997 set_option_value("-pc_bddc_use_faces", "false");
998 if (additional_data.symmetric)
999 set_option_value("-pc_bddc_symmetric", "true");
1000 else
1001 set_option_value("-pc_bddc_symmetric", "false");
1002 if (additional_data.coords.size() > 0)
1003 {
1004 set_option_value("-pc_bddc_corner_selection", "true");
1005 // Convert coords vector to PETSc data array
1006 std::vector<PetscReal> coords_petsc(additional_data.coords.size() *
1007 dim);
1008 for (unsigned int i = 0, j = 0; i < additional_data.coords.size(); ++i)
1009 {
1010 for (j = 0; j < dim; ++j)
1011 coords_petsc[dim * i + j] = additional_data.coords[i][j];
1012 }
1013
1014 ierr = PCSetCoordinates(pc,
1015 dim,
1016 additional_data.coords.size(),
1017 coords_petsc.data());
1018 AssertThrow(ierr == 0, ExcPETScError(ierr));
1019 }
1020 else
1021 {
1022 set_option_value("-pc_bddc_corner_selection", "false");
1023 ierr = PCSetCoordinates(pc, 0, 0, nullptr);
1024 AssertThrow(ierr == 0, ExcPETScError(ierr));
1025 }
1026
1027
1028 ierr = PCSetFromOptions(pc);
1029 AssertThrow(ierr == 0, ExcPETScError(ierr));
1030# else
1032 false, ExcMessage("BDDC preconditioner requires PETSc 3.10.0 or newer"));
1033# endif
1034 }
1035
1036
1037
1038 template <int dim>
1039 void
1041 const AdditionalData &additional_data_)
1042 {
1043 clear();
1044
1045 additional_data = additional_data_;
1046
1047 create_pc_with_mat(matrix_);
1048 initialize();
1049 }
1050
1051 /* ----------------- PreconditionShell -------------------- */
1052
1054 {
1055 initialize(matrix);
1056 }
1057
1062
1063 void
1065 {
1066 PetscErrorCode ierr;
1067 if (pc)
1068 {
1069 ierr = PCDestroy(&pc);
1070 AssertThrow(ierr == 0, ExcPETScError(ierr));
1071 }
1073
1074 ierr = PCSetType(pc, PCSHELL);
1075 AssertThrow(ierr == 0, ExcPETScError(ierr));
1076 ierr = PCShellSetContext(pc, static_cast<void *>(this));
1077 AssertThrow(ierr == 0, ExcPETScError(ierr));
1078 ierr = PCShellSetSetUp(pc, PreconditionShell::pcsetup);
1079 AssertThrow(ierr == 0, ExcPETScError(ierr));
1080 ierr = PCShellSetApply(pc, PreconditionShell::pcapply);
1081 AssertThrow(ierr == 0, ExcPETScError(ierr));
1082 ierr = PCShellSetApplyTranspose(pc, PreconditionShell::pcapply_transpose);
1083 AssertThrow(ierr == 0, ExcPETScError(ierr));
1084 ierr = PCShellSetName(pc, "deal.II user solve");
1085 AssertThrow(ierr == 0, ExcPETScError(ierr));
1086 }
1087
1088 void
1090 {
1091 initialize(matrix.get_mpi_communicator());
1092 PetscErrorCode ierr;
1093 ierr = PCSetOperators(pc, matrix, matrix);
1094 AssertThrow(ierr == 0, ExcPETScError(ierr));
1095 }
1096
1097# ifndef PetscCall
1098# define PetscCall(code) \
1099 do \
1100 { \
1101 PetscErrorCode ierr = (code); \
1102 CHKERRQ(ierr); \
1103 } \
1104 while (false)
1105# endif
1106
1107 PetscErrorCode
1109 {
1110 PetscFunctionBeginUser;
1111 // Failed reason is not reset uniformly within the
1112 // interface code of PCSetUp in PETSc.
1113 // We handle it here.
1114 PetscCall(pc_set_failed_reason(ppc, PC_NOERROR));
1115 PetscFunctionReturn(PETSC_SUCCESS);
1116 }
1117
1118 PetscErrorCode
1119 PreconditionShell::pcapply(PC ppc, Vec x, Vec y)
1120 {
1121 void *ctx;
1122
1123 PetscFunctionBeginUser;
1124 PetscCall(PCShellGetContext(ppc, &ctx));
1125
1126 auto *user = static_cast<PreconditionShell *>(ctx);
1127 if (!user->vmult)
1128 SETERRQ(
1129 PetscObjectComm((PetscObject)ppc),
1130 PETSC_ERR_LIB,
1131 "Failure in ::PETScWrappers::PreconditionShell::pcapply. Missing std::function vmult");
1132
1133 VectorBase src(x);
1134 VectorBase dst(y);
1135 const int lineno = __LINE__;
1136 try
1137 {
1138 user->vmult(dst, src);
1139 }
1140 catch (const RecoverableUserCallbackError &)
1141 {
1142 PetscCall(pc_set_failed_reason(ppc, PC_SUBPC_ERROR));
1143 }
1144 catch (...)
1145 {
1146 return PetscError(
1147 PetscObjectComm((PetscObject)ppc),
1148 lineno + 3,
1149 "vmult",
1150 __FILE__,
1151 PETSC_ERR_LIB,
1152 PETSC_ERROR_INITIAL,
1153 "Failure in pcapply from ::PETScWrappers::NonlinearSolver");
1154 }
1156 PetscFunctionReturn(PETSC_SUCCESS);
1157 }
1158
1159 PetscErrorCode
1161 {
1162 void *ctx;
1163
1164 PetscFunctionBeginUser;
1165 PetscCall(PCShellGetContext(ppc, &ctx));
1166
1167 auto *user = static_cast<PreconditionShell *>(ctx);
1168 if (!user->vmultT)
1169 SETERRQ(
1170 PetscObjectComm((PetscObject)ppc),
1171 PETSC_ERR_LIB,
1172 "Failure in ::PETScWrappers::PreconditionShell::pcapply_transpose. Missing std::function vmultT");
1173
1174 VectorBase src(x);
1175 VectorBase dst(y);
1176 const int lineno = __LINE__;
1177 try
1178 {
1179 user->vmultT(dst, src);
1180 }
1181 catch (const RecoverableUserCallbackError &)
1182 {
1183 PetscCall(pc_set_failed_reason(ppc, PC_SUBPC_ERROR));
1184 }
1185 catch (...)
1186 {
1187 return PetscError(
1188 PetscObjectComm((PetscObject)ppc),
1189 lineno + 3,
1190 "vmultT",
1191 __FILE__,
1192 PETSC_ERR_LIB,
1193 PETSC_ERROR_INITIAL,
1194 "Failure in pcapply_transpose from ::PETScWrappers::NonlinearSolver");
1195 }
1197 PetscFunctionReturn(PETSC_SUCCESS);
1198 }
1199
1200
1201} // namespace PETScWrappers
1202
1205
1207
1208#endif // DEAL_II_WITH_PETSC
void Tvmult(VectorBase &dst, const VectorBase &src) const
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:113
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_NOT_IMPLEMENTED()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & RecoverableUserCallbackError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
const MPI_Comm comm
Definition mpi.cc:928
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)
#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)