Reference documentation for deal.II version 9.5.0
\(\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// 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
21
27
28# include <petscconf.h>
29
30# include <cmath>
31
33
34namespace 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
95 {
96 return PetscObjectComm(reinterpret_cast<PetscObject>(pc));
97 }
98
99 void
101 {
102 // only allow the creation of the
103 // preconditioner once
105
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
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
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
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
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
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
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
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",
590 std::to_string(
592
593 set_option_value("-pc_hypre_boomeramg_max_row_sum",
594 std::to_string(additional_data.max_row_sum));
595
596 set_option_value("-pc_hypre_boomeramg_strong_threshold",
597 std::to_string(additional_data.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",
661 std::to_string(additional_data.n_sweeps_coarse));
662
663 set_option_value("-pc_hypre_boomeramg_tol",
664 std::to_string(additional_data.tol));
665 set_option_value("-pc_hypre_boomeramg_max_iter",
666 std::to_string(additional_data.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)
714 : symmetric(symmetric)
715 , n_levels(n_levels)
716 , threshold(threshold)
717 , filter(filter)
718 , output_details(output_details)
719 {}
720
721
722
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
763 "ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
764
765 std::stringstream ssStream;
766
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,
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",
797 std::to_string(additional_data.n_levels));
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
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
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)
918 , symmetric(symmetric)
919 , coords(coords)
920 {}
921
922
923
924 template <int dim>
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));
980 strcmp(current_type, MATIS) == 0,
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
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 {
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 (0)
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(0);
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(0);
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(0);
1200 }
1201
1202
1203} // namespace PETScWrappers
1204
1207
1209
1210#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:112
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & RecoverableUserCallbackError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
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)
const MPI_Comm comm