Reference documentation for deal.II version GIT relicensing-220-g40aaea9664 2024-03-28 18:00: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\}}\)
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 (void)pc;
555 Assert(false,
556 ExcMessage("Your PETSc installation does not include a copy of "
557 "the hypre package necessary for this preconditioner."));
558# endif
559 }
560
561
562
564 const MatrixBase &matrix,
565 const AdditionalData &additional_data)
566 : PreconditionBase(matrix.get_mpi_communicator())
567 {
569 }
570
571
572
573 void
575 {
576# ifdef DEAL_II_PETSC_WITH_HYPRE
577 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
578 AssertThrow(ierr == 0, ExcPETScError(ierr));
579
580 ierr = PCHYPRESetType(pc, "boomeramg");
581 AssertThrow(ierr == 0, ExcPETScError(ierr));
582
584 {
585 set_option_value("-pc_hypre_boomeramg_print_statistics", "1");
586 }
587
588 set_option_value("-pc_hypre_boomeramg_agg_nl",
589 std::to_string(
591
592 set_option_value("-pc_hypre_boomeramg_max_row_sum",
593 std::to_string(additional_data.max_row_sum));
594
595 set_option_value("-pc_hypre_boomeramg_strong_threshold",
596 std::to_string(additional_data.strong_threshold));
597
598 // change to symmetric SOR/Jacobi when using a symmetric operator for
599 // backward compatibility
603 {
606 }
607
611 {
614 }
615
619 {
622 }
623
624 auto relaxation_type_is_symmetric =
625 [](AdditionalData::RelaxationType relaxation_type) {
626 return relaxation_type == AdditionalData::RelaxationType::Jacobi ||
627 relaxation_type ==
629 relaxation_type ==
631 relaxation_type == AdditionalData::RelaxationType::None ||
632 relaxation_type ==
634 relaxation_type == AdditionalData::RelaxationType::CG ||
636 };
637
639 !relaxation_type_is_symmetric(additional_data.relaxation_type_up))
640 Assert(false,
641 ExcMessage("Use a symmetric smoother for relaxation_type_up"));
642
644 !relaxation_type_is_symmetric(additional_data.relaxation_type_down))
645 Assert(false,
646 ExcMessage("Use a symmetric smoother for relaxation_type_down"));
647
649 !relaxation_type_is_symmetric(additional_data.relaxation_type_coarse))
650 Assert(false,
651 ExcMessage("Use a symmetric smoother for relaxation_type_coarse"));
652
653 set_option_value("-pc_hypre_boomeramg_relax_type_up",
655 set_option_value("-pc_hypre_boomeramg_relax_type_down",
657 set_option_value("-pc_hypre_boomeramg_relax_type_coarse",
659 set_option_value("-pc_hypre_boomeramg_grid_sweeps_coarse",
660 std::to_string(additional_data.n_sweeps_coarse));
661
662 set_option_value("-pc_hypre_boomeramg_tol",
663 std::to_string(additional_data.tol));
664 set_option_value("-pc_hypre_boomeramg_max_iter",
665 std::to_string(additional_data.max_iter));
666
668 {
669 set_option_value("-pc_hypre_boomeramg_cycle_type", "W");
670 }
671
672 ierr = PCSetFromOptions(pc);
673 AssertThrow(ierr == 0, ExcPETScError(ierr));
674# else
675 Assert(false,
676 ExcMessage("Your PETSc installation does not include a copy of "
677 "the hypre package necessary for this preconditioner."));
678# endif
679 }
680
681
682
683 void
685 const AdditionalData &additional_data_)
686 {
687# ifdef DEAL_II_PETSC_WITH_HYPRE
688 clear();
689
690 additional_data = additional_data_;
691
692 create_pc_with_mat(matrix_);
693 initialize();
694
695# else // DEAL_II_PETSC_WITH_HYPRE
696 (void)matrix_;
697 (void)additional_data_;
698 Assert(false,
699 ExcMessage("Your PETSc installation does not include a copy of "
700 "the hypre package necessary for this preconditioner."));
701# endif
702 }
703
704
705 /* ----------------- PreconditionParaSails -------------------- */
706
708 const unsigned int symmetric,
709 const unsigned int n_levels,
710 const double threshold,
711 const double filter,
712 const bool output_details)
713 : symmetric(symmetric)
714 , n_levels(n_levels)
715 , threshold(threshold)
716 , filter(filter)
717 , output_details(output_details)
718 {}
719
720
721
725
726
727
729 const MatrixBase &matrix,
730 const AdditionalData &additional_data)
731 : PreconditionBase(matrix.get_mpi_communicator())
732 {
734 }
735
736
737 void
739 const AdditionalData &additional_data_)
740 {
741 clear();
742
743 additional_data = additional_data_;
744
745# ifdef DEAL_II_PETSC_WITH_HYPRE
746 create_pc_with_mat(matrix_);
747
748 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
749 AssertThrow(ierr == 0, ExcPETScError(ierr));
750
751 ierr = PCHYPRESetType(pc, "parasails");
752 AssertThrow(ierr == 0, ExcPETScError(ierr));
753
755 {
756 set_option_value("-pc_hypre_parasails_logging", "1");
757 }
758
762 "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
763
764 std::stringstream ssStream;
765
767 {
768 case 0:
769 {
770 ssStream << "nonsymmetric";
771 break;
772 }
773
774 case 1:
775 {
776 ssStream << "SPD";
777 break;
778 }
779
780 case 2:
781 {
782 ssStream << "nonsymmetric,SPD";
783 break;
784 }
785
786 default:
787 Assert(
788 false,
790 "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
791 }
792
793 set_option_value("-pc_hypre_parasails_sym", ssStream.str());
794
795 set_option_value("-pc_hypre_parasails_nlevels",
796 std::to_string(additional_data.n_levels));
797
798 ssStream.str(""); // empty the stringstream
799 ssStream << additional_data.threshold;
800 set_option_value("-pc_hypre_parasails_thresh", ssStream.str());
801
802 ssStream.str(""); // empty the stringstream
803 ssStream << additional_data.filter;
804 set_option_value("-pc_hypre_parasails_filter", ssStream.str());
805
806 ierr = PCSetFromOptions(pc);
807 AssertThrow(ierr == 0, ExcPETScError(ierr));
808
809# else // DEAL_II_PETSC_WITH_HYPRE
810 (void)matrix_;
811 Assert(false,
812 ExcMessage("Your PETSc installation does not include a copy of "
813 "the hypre package necessary for this preconditioner."));
814# endif
815 }
816
817
818 /* ----------------- PreconditionNone ------------------------- */
819
823
824
825
827 const AdditionalData &additional_data)
828 : PreconditionBase(matrix.get_mpi_communicator())
829 {
831 }
832
833
834 void
836 const AdditionalData &additional_data_)
837 {
838 clear();
839
840 additional_data = additional_data_;
841
842 create_pc_with_mat(matrix_);
843
844 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCNONE));
845 AssertThrow(ierr == 0, ExcPETScError(ierr));
846
847 ierr = PCSetFromOptions(pc);
848 AssertThrow(ierr == 0, ExcPETScError(ierr));
849 }
850
851
852 /* ----------------- PreconditionLU -------------------- */
853
855 const double zero_pivot,
856 const double damping)
857 : pivoting(pivoting)
858 , zero_pivot(zero_pivot)
859 , damping(damping)
860 {}
861
862
863
867
868
869
871 const AdditionalData &additional_data)
872 : PreconditionBase(matrix.get_mpi_communicator())
873 {
875 }
876
877
878 void
880 const AdditionalData &additional_data_)
881 {
882 clear();
883
884 additional_data = additional_data_;
885
886 create_pc_with_mat(matrix_);
887
888 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCLU));
889 AssertThrow(ierr == 0, ExcPETScError(ierr));
890
891 // set flags as given
892 ierr = PCFactorSetColumnPivot(pc, additional_data.pivoting);
893 AssertThrow(ierr == 0, ExcPETScError(ierr));
894
895 ierr = PCFactorSetZeroPivot(pc, additional_data.zero_pivot);
896 AssertThrow(ierr == 0, ExcPETScError(ierr));
897
898 ierr = PCFactorSetShiftAmount(pc, additional_data.damping);
899 AssertThrow(ierr == 0, ExcPETScError(ierr));
900
901 ierr = PCSetFromOptions(pc);
902 AssertThrow(ierr == 0, ExcPETScError(ierr));
903 }
904
905 /* ----------------- PreconditionBDDC -------------------- */
906
907 template <int dim>
909 const bool use_vertices,
910 const bool use_edges,
911 const bool use_faces,
912 const bool symmetric,
913 const std::vector<Point<dim>> coords)
914 : use_vertices(use_vertices)
915 , use_edges(use_edges)
916 , use_faces(use_faces)
917 , symmetric(symmetric)
918 , coords(coords)
919 {}
920
921
922
923 template <int dim>
927
928
929
930 template <int dim>
932 const MPI_Comm comm,
933 const AdditionalData &additional_data_)
935 {
936 additional_data = additional_data_;
937
938 PetscErrorCode ierr = PCCreate(comm, &pc);
939 AssertThrow(ierr == 0, ExcPETScError(ierr));
940
941 initialize();
942 }
943
944
945
946 template <int dim>
948 const AdditionalData &additional_data)
949 : PreconditionBase(matrix.get_mpi_communicator())
950 {
952 }
953
954
955
956 template <int dim>
957 void
959 {
960# if DEAL_II_PETSC_VERSION_GTE(3, 10, 0)
961 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCBDDC));
962 AssertThrow(ierr == 0, ExcPETScError(ierr));
963
964 // The matrix must be of IS type. We check for this to avoid the PETSc error
965 // in order to suggest the correct matrix reinit method.
966 {
967 MatType current_type;
968 Mat A, P;
969 PetscBool flg;
970
971 ierr = PCGetOperators(pc, &A, &P);
972 AssertThrow(ierr == 0, ExcPETScError(ierr));
973 ierr = PCGetUseAmat(pc, &flg);
974 AssertThrow(ierr == 0, ExcPETScError(ierr));
975
976 ierr = MatGetType(flg ? A : P, &current_type);
977 AssertThrow(ierr == 0, ExcPETScError(ierr));
979 strcmp(current_type, MATIS) == 0,
981 "Matrix must be of IS type. For this, the variant of reinit that includes the active dofs must be used."));
982 }
983
984
985 std::stringstream ssStream;
986
987 if (additional_data.use_vertices)
988 set_option_value("-pc_bddc_use_vertices", "true");
989 else
990 set_option_value("-pc_bddc_use_vertices", "false");
991 if (additional_data.use_edges)
992 set_option_value("-pc_bddc_use_edges", "true");
993 else
994 set_option_value("-pc_bddc_use_edges", "false");
995 if (additional_data.use_faces)
996 set_option_value("-pc_bddc_use_faces", "true");
997 else
998 set_option_value("-pc_bddc_use_faces", "false");
999 if (additional_data.symmetric)
1000 set_option_value("-pc_bddc_symmetric", "true");
1001 else
1002 set_option_value("-pc_bddc_symmetric", "false");
1003 if (additional_data.coords.size() > 0)
1004 {
1005 set_option_value("-pc_bddc_corner_selection", "true");
1006 // Convert coords vector to PETSc data array
1007 std::vector<PetscReal> coords_petsc(additional_data.coords.size() *
1008 dim);
1009 for (unsigned int i = 0, j = 0; i < additional_data.coords.size(); ++i)
1010 {
1011 for (j = 0; j < dim; ++j)
1012 coords_petsc[dim * i + j] = additional_data.coords[i][j];
1013 }
1014
1015 ierr = PCSetCoordinates(pc,
1016 dim,
1017 additional_data.coords.size(),
1018 coords_petsc.data());
1019 AssertThrow(ierr == 0, ExcPETScError(ierr));
1020 }
1021 else
1022 {
1023 set_option_value("-pc_bddc_corner_selection", "false");
1024 ierr = PCSetCoordinates(pc, 0, 0, nullptr);
1025 AssertThrow(ierr == 0, ExcPETScError(ierr));
1026 }
1027
1028
1029 ierr = PCSetFromOptions(pc);
1030 AssertThrow(ierr == 0, ExcPETScError(ierr));
1031# else
1033 false, ExcMessage("BDDC preconditioner requires PETSc 3.10.0 or newer"));
1034# endif
1035 }
1036
1037
1038
1039 template <int dim>
1040 void
1042 const AdditionalData &additional_data_)
1043 {
1044 clear();
1045
1046 additional_data = additional_data_;
1047
1048 create_pc_with_mat(matrix_);
1049 initialize();
1050 }
1051
1052 /* ----------------- PreconditionShell -------------------- */
1053
1055 {
1056 initialize(matrix);
1057 }
1058
1063
1064 void
1066 {
1067 PetscErrorCode ierr;
1068 if (pc)
1069 {
1070 ierr = PCDestroy(&pc);
1071 AssertThrow(ierr == 0, ExcPETScError(ierr));
1072 }
1074
1075 ierr = PCSetType(pc, PCSHELL);
1076 AssertThrow(ierr == 0, ExcPETScError(ierr));
1077 ierr = PCShellSetContext(pc, static_cast<void *>(this));
1078 AssertThrow(ierr == 0, ExcPETScError(ierr));
1079 ierr = PCShellSetSetUp(pc, PreconditionShell::pcsetup);
1080 AssertThrow(ierr == 0, ExcPETScError(ierr));
1081 ierr = PCShellSetApply(pc, PreconditionShell::pcapply);
1082 AssertThrow(ierr == 0, ExcPETScError(ierr));
1083 ierr = PCShellSetApplyTranspose(pc, PreconditionShell::pcapply_transpose);
1084 AssertThrow(ierr == 0, ExcPETScError(ierr));
1085 ierr = PCShellSetName(pc, "deal.II user solve");
1086 AssertThrow(ierr == 0, ExcPETScError(ierr));
1087 }
1088
1089 void
1091 {
1092 initialize(matrix.get_mpi_communicator());
1093 PetscErrorCode ierr;
1094 ierr = PCSetOperators(pc, matrix, matrix);
1095 AssertThrow(ierr == 0, ExcPETScError(ierr));
1096 }
1097
1098# ifndef PetscCall
1099# define PetscCall(code) \
1100 do \
1101 { \
1102 PetscErrorCode ierr = (code); \
1103 CHKERRQ(ierr); \
1104 } \
1105 while (false)
1106# endif
1107
1108 PetscErrorCode
1110 {
1111 PetscFunctionBeginUser;
1112 // Failed reason is not reset uniformly within the
1113 // interface code of PCSetUp in PETSc.
1114 // We handle it here.
1115 PetscCall(pc_set_failed_reason(ppc, PC_NOERROR));
1116 PetscFunctionReturn(PETSC_SUCCESS);
1117 }
1118
1119 PetscErrorCode
1120 PreconditionShell::pcapply(PC ppc, Vec x, Vec y)
1121 {
1122 void *ctx;
1123
1124 PetscFunctionBeginUser;
1125 PetscCall(PCShellGetContext(ppc, &ctx));
1126
1127 auto *user = static_cast<PreconditionShell *>(ctx);
1128 if (!user->vmult)
1129 SETERRQ(
1130 PetscObjectComm((PetscObject)ppc),
1131 PETSC_ERR_LIB,
1132 "Failure in ::PETScWrappers::PreconditionShell::pcapply. Missing std::function vmult");
1133
1134 VectorBase src(x);
1135 VectorBase dst(y);
1136 const int lineno = __LINE__;
1137 try
1138 {
1139 user->vmult(dst, src);
1140 }
1141 catch (const RecoverableUserCallbackError &)
1142 {
1143 PetscCall(pc_set_failed_reason(ppc, PC_SUBPC_ERROR));
1144 }
1145 catch (...)
1146 {
1147 return PetscError(
1148 PetscObjectComm((PetscObject)ppc),
1149 lineno + 3,
1150 "vmult",
1151 __FILE__,
1152 PETSC_ERR_LIB,
1153 PETSC_ERROR_INITIAL,
1154 "Failure in pcapply from ::PETScWrappers::NonlinearSolver");
1155 }
1157 PetscFunctionReturn(PETSC_SUCCESS);
1158 }
1159
1160 PetscErrorCode
1162 {
1163 void *ctx;
1164
1165 PetscFunctionBeginUser;
1166 PetscCall(PCShellGetContext(ppc, &ctx));
1167
1168 auto *user = static_cast<PreconditionShell *>(ctx);
1169 if (!user->vmultT)
1170 SETERRQ(
1171 PetscObjectComm((PetscObject)ppc),
1172 PETSC_ERR_LIB,
1173 "Failure in ::PETScWrappers::PreconditionShell::pcapply_transpose. Missing std::function vmultT");
1174
1175 VectorBase src(x);
1176 VectorBase dst(y);
1177 const int lineno = __LINE__;
1178 try
1179 {
1180 user->vmultT(dst, src);
1181 }
1182 catch (const RecoverableUserCallbackError &)
1183 {
1184 PetscCall(pc_set_failed_reason(ppc, PC_SUBPC_ERROR));
1185 }
1186 catch (...)
1187 {
1188 return PetscError(
1189 PetscObjectComm((PetscObject)ppc),
1190 lineno + 3,
1191 "vmultT",
1192 __FILE__,
1193 PETSC_ERR_LIB,
1194 PETSC_ERROR_INITIAL,
1195 "Failure in pcapply_transpose from ::PETScWrappers::NonlinearSolver");
1196 }
1198 PetscFunctionReturn(PETSC_SUCCESS);
1199 }
1200
1201
1202} // namespace PETScWrappers
1203
1206
1208
1209#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:111
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & RecoverableUserCallbackError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
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)
*braid_SplitCommworld & comm
#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)