Reference documentation for deal.II version GIT relicensing-822-gf1642d4ea5 2024-06-12 14:10: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
44
46 {
47 try
48 {
49 clear();
50 }
51 catch (...)
52 {}
53 }
54
55 void
57 {
58 if (pc)
59 {
62 }
63 }
64
65 void
73
74 void
82
83 void
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
107 reinterpret_cast<PetscObject>(static_cast<const Mat &>(matrix)), &comm);
109
111
112 ierr = PCSetOperators(pc, matrix, matrix);
114 }
115
116 void
123
124 const PC &
126 {
127 return pc;
128 }
129
130
131 /* ----------------- PreconditionJacobi -------------------- */
132
136
137
138
150
151
152
154 const AdditionalData &additional_data)
155 : PreconditionBase(matrix.get_mpi_communicator())
156 {
158 }
159
160
161
162 void
173
174
175
176 void
187
188
189 /* ----------------- PreconditionBlockJacobi -------------------- */
190
194
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
229
230
231
232 void
243
244
245 /* ----------------- PreconditionSOR -------------------- */
246
250
251
252
254 : omega(omega)
255 {}
256
257
258
265
266
267 void
270 {
271 clear();
272
274
276
277 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
279
280 // then set flags as given
283
286 }
287
288
289 /* ----------------- PreconditionSSOR -------------------- */
290
294
295
296
298 : omega(omega)
299 {}
300
301
302
309
310
311 void
314 {
315 clear();
316
318
320
321 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCSOR));
323
324 // then set flags as given
327
328 // convert SOR to SSOR
331
334 }
335
336
337 /* ----------------- PreconditionICC -------------------- */
338
342
343
344
346 : levels(levels)
347 {}
348
349
350
357
358
359 void
379
380
381 /* ----------------- PreconditionILU -------------------- */
382
386
387
388
390 : levels(levels)
391 {}
392
393
394
401
402
403 void
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,
545 {
547
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));
579
580 ierr = PCHYPRESetType(pc, "boomeramg");
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
636 };
637
640 Assert(false,
641 ExcMessage("Use a symmetric smoother for relaxation_type_up"));
642
645 Assert(false,
646 ExcMessage("Use a symmetric smoother for relaxation_type_down"));
647
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
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
686 {
687# ifdef DEAL_II_PETSC_WITH_HYPRE
688 clear();
689
691
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
740 {
741 clear();
742
744
745# ifdef DEAL_II_PETSC_WITH_HYPRE
747
748 PetscErrorCode ierr = PCSetType(pc, const_cast<char *>(PCHYPRE));
750
751 ierr = PCHYPRESetType(pc, "parasails");
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
800 set_option_value("-pc_hypre_parasails_thresh", ssStream.str());
801
802 ssStream.str(""); // empty the stringstream
804 set_option_value("-pc_hypre_parasails_filter", ssStream.str());
805
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
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
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>
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));
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 {
968 Mat A, P;
970
971 ierr = PCGetOperators(pc, &A, &P);
973 ierr = PCGetUseAmat(pc, &flg);
975
976 ierr = MatGetType(flg ? A : P, &current_type);
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
1016 dim,
1017 additional_data.coords.size(),
1018 coords_petsc.data());
1020 }
1021 else
1022 {
1023 set_option_value("-pc_bddc_corner_selection", "false");
1024 ierr = PCSetCoordinates(pc, 0, 0, nullptr);
1026 }
1027
1028
1029 ierr = PCSetFromOptions(pc);
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
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
1088
1089 void
1091 {
1092 initialize(matrix.get_mpi_communicator());
1094 ierr = PCSetOperators(pc, matrix, matrix);
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
1110 {
1112 // Failed reason is not reset uniformly within the
1113 // interface code of PCSetUp in PETSc.
1114 // We handle it here.
1117 }
1118
1121 {
1122 void *ctx;
1123
1126
1127 auto *user = static_cast<PreconditionShell *>(ctx);
1128 if (!user->vmult)
1129 SETERRQ(
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 {
1144 }
1145 catch (...)
1146 {
1147 return PetscError(
1149 lineno + 3,
1150 "vmult",
1151 __FILE__,
1154 "Failure in pcapply from ::PETScWrappers::NonlinearSolver");
1155 }
1158 }
1159
1162 {
1163 void *ctx;
1164
1167
1168 auto *user = static_cast<PreconditionShell *>(ctx);
1169 if (!user->vmultT)
1170 SETERRQ(
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 {
1185 }
1186 catch (...)
1187 {
1188 return PetscError(
1190 lineno + 3,
1191 "vmultT",
1192 __FILE__,
1195 "Failure in pcapply_transpose from ::PETScWrappers::NonlinearSolver");
1196 }
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)
#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)