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_solver.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
16
18
20
21#ifdef DEAL_II_WITH_PETSC
22
28
29// Shorthand notation for PETSc error codes.
30# define AssertPETSc(code) \
31 do \
32 { \
33 PetscErrorCode ierr = (code); \
34 AssertThrow(ierr == 0, ExcPETScError(ierr)); \
35 } \
36 while (0)
37
39
40namespace PETScWrappers
41{
43 : ksp(nullptr)
44 , solver_control(nullptr)
45 {}
46
47
49 : ksp(nullptr)
50 , solver_control(&cn)
51 {}
52
53
54
55 void
57 {}
58
59
60
62 {
63 AssertPETSc(KSPDestroy(&ksp));
64 }
65
66
67
68 KSP
70 {
71 return ksp;
72 }
73
74
75
76 SolverBase::operator KSP() const
77 {
78 return ksp;
79 }
80
81
82
83 void
85 VectorBase & x,
86 const VectorBase & b,
87 const PreconditionBase &preconditioner)
88 {
89 // first create a solver object if this
90 // is necessary
91 if (ksp == nullptr)
92 {
94
95 // let derived classes set the solver
96 // type, and the preconditioning
97 // object set the type of
98 // preconditioner
100
101 AssertPETSc(KSPSetPC(ksp, preconditioner.get_pc()));
102
103 /*
104 * by default we set up the preconditioner only once.
105 * this can be overridden by command line.
106 */
107 AssertPETSc(KSPSetReusePreconditioner(ksp, PETSC_TRUE));
108 }
109
110 // setting the preconditioner overwrites the used matrices.
111 // hence, we need to set the matrices after the preconditioner.
112 Mat B;
113 AssertPETSc(KSPGetOperators(ksp, nullptr, &B));
114 AssertPETSc(KSPSetOperators(ksp, A, B));
115
116 // set the command line option prefix name
117 AssertPETSc(KSPSetOptionsPrefix(ksp, prefix_name.c_str()));
118
119 // set the command line options provided
120 // by the user to override the defaults
121 AssertPETSc(KSPSetFromOptions(ksp));
122
123 // then do the real work: set up solver
124 // internal data and solve the
125 // system.
126 AssertPETSc(KSPSetUp(ksp));
127
128 AssertPETSc(KSPSolve(ksp, b, x));
129
130 // in case of failure: throw
131 // exception
132 if (solver_control &&
133 solver_control->last_check() != SolverControl::success)
134 AssertThrow(false,
136 solver_control->last_value()));
137 // otherwise exit as normal
138 }
139
140
141 void
142 SolverBase::set_prefix(const std::string &prefix)
143 {
144 prefix_name = prefix;
145 }
146
147
148 void
150 {
151 AssertPETSc(KSPDestroy(&ksp));
152 }
153
154
157 {
161 "You need to create the solver with a SolverControl object if you want to call the function that returns it."));
162 return *solver_control;
163 }
164
165
166 int
168 const PetscInt iteration,
169 const PetscReal residual_norm,
170 KSPConvergedReason *reason,
171 void * solver_control_x)
172 {
174 *reinterpret_cast<SolverControl *>(solver_control_x);
175
176 const SolverControl::State state =
177 solver_control.check(iteration, residual_norm);
178
179 switch (state)
180 {
181 case ::SolverControl::iterate:
182 *reason = KSP_CONVERGED_ITERATING;
183 break;
184
185 case ::SolverControl::success:
186 *reason = KSP_CONVERGED_RTOL;
187 break;
188
189 case ::SolverControl::failure:
190 if (solver_control.last_step() > solver_control.max_steps())
191 *reason = KSP_DIVERGED_ITS;
192 else
193 *reason = KSP_DIVERGED_DTOL;
194 break;
195
196 default:
197 Assert(false, ExcNotImplemented());
198 }
199
200 // return without failure
201 return 0;
202 }
203
204
205
206 void
208 {
209 // Create the PETSc KSP object
210 AssertPETSc(KSPCreate(comm, &ksp));
211
212 // then a convergence monitor
213 // function that simply
214 // checks with the solver_control
215 // object we have in this object for
216 // convergence
218 }
219
220
221
222 void
224 {
225 if (ksp && solver_control)
227 KSPSetConvergenceTest(ksp, &convergence_test, solver_control, nullptr));
228 }
229
230
231 void
233 {
235
236 // let derived classes set the solver
237 // type, and the preconditioning
238 // object set the type of
239 // preconditioner
241
242 // set the command line options provided
243 // by the user to override the defaults
244 AssertPETSc(KSPSetFromOptions(ksp));
245 }
246
247
248
249 /* ---------------------- SolverRichardson ------------------------ */
250
252 : omega(omega)
253 {}
254
255
256
258 const AdditionalData &data)
259 : SolverBase(cn)
260 , additional_data(data)
261 {}
262
263
264
266 const MPI_Comm,
267 const AdditionalData &data)
268 : SolverRichardson(cn, data)
269 {}
270
271
272 void
274 {
275 AssertPETSc(KSPSetType(ksp, KSPRICHARDSON));
276
277 // set the damping factor from the data
278 AssertPETSc(KSPRichardsonSetScale(ksp, additional_data.omega));
279
280 // in the deal.II solvers, we always
281 // honor the initial guess in the
282 // solution vector. do so here as well:
283 AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
284
285 // Hand over the absolute
286 // tolerance and the maximum
287 // iteration number to the PETSc
288 // convergence criterion. The
289 // custom deal.II SolverControl
290 // object is ignored by the PETSc
291 // Richardson method (when no
292 // PETSc monitoring is present),
293 // since in this case PETSc
294 // uses a faster version of
295 // the Richardson iteration,
296 // where no residual is
297 // available.
298 AssertPETSc(KSPSetTolerances(ksp,
299 PETSC_DEFAULT,
300 this->solver_control->tolerance(),
301 PETSC_DEFAULT,
302 this->solver_control->max_steps() + 1));
303 }
304
305
306 /* ---------------------- SolverChebychev ------------------------ */
307
309 const AdditionalData &data)
310 : SolverBase(cn)
311 , additional_data(data)
312 {}
313
314
316 const MPI_Comm,
317 const AdditionalData &data)
318 : SolverChebychev(cn, data)
319 {}
320
321
322 void
324 {
325 AssertPETSc(KSPSetType(ksp, KSPCHEBYSHEV));
326
327 // in the deal.II solvers, we always
328 // honor the initial guess in the
329 // solution vector. do so here as well:
330 AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
331 }
332
333
334 /* ---------------------- SolverCG ------------------------ */
335
337 : SolverBase(cn)
338 , additional_data(data)
339 {}
340
341
343 const MPI_Comm,
344 const AdditionalData &data)
345 : SolverCG(cn, data)
346 {}
347
348
349 void
351 {
352 AssertPETSc(KSPSetType(ksp, KSPCG));
353
354 // in the deal.II solvers, we always
355 // honor the initial guess in the
356 // solution vector. do so here as well:
357 AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
358 }
359
360
361 /* ---------------------- SolverBiCG ------------------------ */
362
364 : SolverBase(cn)
365 , additional_data(data)
366 {}
367
368
370 const MPI_Comm,
371 const AdditionalData &data)
372 : SolverBiCG(cn, data)
373 {}
374
375
376 void
378 {
379 AssertPETSc(KSPSetType(ksp, KSPBICG));
380
381 // in the deal.II solvers, we always
382 // honor the initial guess in the
383 // solution vector. do so here as well:
384 AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
385 }
386
387
388 /* ---------------------- SolverGMRES ------------------------ */
389
391 const unsigned int restart_parameter,
392 const bool right_preconditioning)
393 : restart_parameter(restart_parameter)
394 , right_preconditioning(right_preconditioning)
395 {}
396
397
398
400 : SolverBase(cn)
401 , additional_data(data)
402 {}
403
404
406 const MPI_Comm,
407 const AdditionalData &data)
408 : SolverGMRES(cn, data)
409 {}
410
411
412 void
414 {
415 AssertPETSc(KSPSetType(ksp, KSPGMRES));
416
417 AssertPETSc(KSPGMRESSetRestart(ksp, additional_data.restart_parameter));
418
419 // Set preconditioning side to right
421 {
422 AssertPETSc(KSPSetPCSide(ksp, PC_RIGHT));
423 }
424
425 // in the deal.II solvers, we always
426 // honor the initial guess in the
427 // solution vector. do so here as well:
428 AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
429 }
430
431
432 /* ---------------------- SolverBicgstab ------------------------ */
433
435 : SolverBase(cn)
436 , additional_data(data)
437 {}
438
439
441 const MPI_Comm,
442 const AdditionalData &data)
443 : SolverBicgstab(cn, data)
444 {}
445
446
447 void
449 {
450 AssertPETSc(KSPSetType(ksp, KSPBCGS));
451
452 // in the deal.II solvers, we always
453 // honor the initial guess in the
454 // solution vector. do so here as well:
455 AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
456 }
457
458
459 /* ---------------------- SolverCGS ------------------------ */
460
462 : SolverBase(cn)
463 , additional_data(data)
464 {}
465
466
468 const MPI_Comm,
469 const AdditionalData &data)
470 : SolverCGS(cn, data)
471 {}
472
473
474 void
476 {
477 AssertPETSc(KSPSetType(ksp, KSPCGS));
478
479 // in the deal.II solvers, we always
480 // honor the initial guess in the
481 // solution vector. do so here as well:
482 AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
483 }
484
485
486 /* ---------------------- SolverTFQMR ------------------------ */
487
489 : SolverBase(cn)
490 , additional_data(data)
491 {}
492
493
495 const MPI_Comm,
496 const AdditionalData &data)
497 : SolverTFQMR(cn, data)
498 {}
499
500
501 void
503 {
504 AssertPETSc(KSPSetType(ksp, KSPTFQMR));
505
506 // in the deal.II solvers, we always
507 // honor the initial guess in the
508 // solution vector. do so here as well:
509 AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
510 }
511
512
513 /* ---------------------- SolverTCQMR ------------------------ */
514
516 : SolverBase(cn)
517 , additional_data(data)
518 {}
519
520
522 const MPI_Comm,
523 const AdditionalData &data)
524 : SolverTCQMR(cn, data)
525 {}
526
527
528 void
530 {
531 AssertPETSc(KSPSetType(ksp, KSPTCQMR));
532
533 // in the deal.II solvers, we always
534 // honor the initial guess in the
535 // solution vector. do so here as well:
536 AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
537 }
538
539
540 /* ---------------------- SolverCR ------------------------ */
541
543 : SolverBase(cn)
544 , additional_data(data)
545 {}
546
547
549 const MPI_Comm,
550 const AdditionalData &data)
551 : SolverCR(cn, data)
552 {}
553
554
555 void
557 {
558 AssertPETSc(KSPSetType(ksp, KSPCR));
559
560 // in the deal.II solvers, we always
561 // honor the initial guess in the
562 // solution vector. do so here as well:
563 AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
564 }
565
566
567 /* ---------------------- SolverLSQR ------------------------ */
568
570 : SolverBase(cn)
571 , additional_data(data)
572 {}
573
574
575
577 const MPI_Comm,
578 const AdditionalData &data)
579 : SolverLSQR(cn, data)
580 {}
581
582
583
584 void
586 {
587 AssertPETSc(KSPSetType(ksp, KSPLSQR));
588
589 // in the deal.II solvers, we always
590 // honor the initial guess in the
591 // solution vector. do so here as well:
592 AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
593
594 // The KSPLSQR implementation overwrites the user-defined
595 // convergence test at creation (i.e. KSPSetType) time.
596 // This is probably a bad design decision in PETSc.
597 // Anyway, here we make sure we use our own convergence
598 // test.
600 }
601
602
603 /* ---------------------- SolverPreOnly ------------------------ */
604
606 : SolverBase(cn)
607 , additional_data(data)
608 {}
609
610
612 const MPI_Comm,
613 const AdditionalData &data)
614 : SolverPreOnly(cn, data)
615 {}
616
617
618 void
620 {
621 AssertPETSc(KSPSetType(ksp, KSPPREONLY));
622
623 // The KSPPREONLY solver of
624 // PETSc never calls the convergence
625 // monitor, which leads to failure
626 // even when everything was ok.
627 // Therefore the SolverControl status
628 // is set to some nice values, which
629 // guarantee a nice result at the end
630 // of the solution process.
631 solver_control->check(1, 0.0);
632
633 // Using the PREONLY solver with
634 // a nonzero initial guess leads
635 // PETSc to produce some error messages.
636 AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_FALSE));
637 }
638
639
640 /* ---------------------- SparseDirectMUMPS------------------------ */
641
643 const AdditionalData &data)
644 : SolverBase(cn)
645 , additional_data(data)
646 , symmetric_mode(false)
647 {}
648
649
650
652 const MPI_Comm,
653 const AdditionalData &data)
654 : SparseDirectMUMPS(cn, data)
655 {}
656
657
658
659 void
661 {
662 /*
663 * KSPPREONLY implements a stub method that applies only the
664 * preconditioner. Its use is due to SparseDirectMUMPS being a direct
665 * (rather than iterative) solver
666 */
667 AssertPETSc(KSPSetType(ksp, KSPPREONLY));
668
669 /*
670 * The KSPPREONLY solver of PETSc never calls the convergence monitor,
671 * which leads to failure even when everything was ok. Therefore, the
672 * SolverControl status is set to some nice values, which guarantee a
673 * nice result at the end of the solution process.
674 */
675 solver_control->check(1, 0.0);
676
677 /*
678 * Using a PREONLY solver with a nonzero initial guess leads PETSc to
679 * produce some error messages.
680 */
681 AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_FALSE));
682 }
683
684 void
686 VectorBase & x,
687 const VectorBase &b)
688 {
689# ifdef DEAL_II_PETSC_WITH_MUMPS
690 /*
691 * creating a solver object if this is necessary
692 */
693 if (ksp == nullptr)
694 {
696
697 /*
698 * setting the solver type
699 */
701
702 /*
703 * set the matrices involved. the last argument is irrelevant here,
704 * since we use the solver only once anyway
705 */
706 AssertPETSc(KSPSetOperators(ksp, A, A));
707
708 /*
709 * getting the associated preconditioner context
710 */
711 PC pc;
712 AssertPETSc(KSPGetPC(ksp, &pc));
713
714 /*
715 * build PETSc PC for particular PCLU or PCCHOLESKY preconditioner
716 * depending on whether the symmetric mode has been set
717 */
718 if (symmetric_mode)
719 AssertPETSc(PCSetType(pc, PCCHOLESKY));
720 else
721 AssertPETSc(PCSetType(pc, PCLU));
722
723 /*
724 * set the software that is to be used to perform the lu
725 * factorization here we start to see differences with the base
726 * class solve function
727 */
728# if DEAL_II_PETSC_VERSION_LT(3, 9, 0)
729 AssertPETSc(PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS));
730# else
731 AssertPETSc(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS));
732# endif
733
734 /*
735 * set up the package to call for the factorization
736 */
737# if DEAL_II_PETSC_VERSION_LT(3, 9, 0)
738 AssertPETSc(PCFactorSetUpMatSolverPackage(pc));
739# else
740 AssertPETSc(PCFactorSetUpMatSolverType(pc));
741# endif
742
743 /*
744 * get the factored matrix F from the preconditioner context.
745 */
746 Mat F;
747 AssertPETSc(PCFactorGetMatrix(pc, &F));
748
749 /*
750 * pass control parameters to MUMPS.
751 * Setting entry 7 of MUMPS ICNTL array to a value
752 * of 2. This sets use of Approximate Minimum Fill (AMF)
753 */
754 AssertPETSc(MatMumpsSetIcntl(F, 7, 2));
755
756 /*
757 * by default we set up the preconditioner only once.
758 * this can be overridden by command line.
759 */
760 AssertPETSc(KSPSetReusePreconditioner(ksp, PETSC_TRUE));
761 }
762
763 /*
764 * set the matrices involved. the last argument is irrelevant here,
765 * since we use the solver only once anyway
766 */
767 AssertPETSc(KSPSetOperators(ksp, A, A));
768
769 /*
770 * set the command line option prefix name
771 */
772 AssertPETSc(KSPSetOptionsPrefix(ksp, prefix_name.c_str()));
773
774 /*
775 * set the command line options provided by the user to override
776 * the defaults
777 */
778 AssertPETSc(KSPSetFromOptions(ksp));
779
780 /*
781 * solve the linear system
782 */
783 AssertPETSc(KSPSolve(ksp, b, x));
784
785 /*
786 * in case of failure throw exception
787 */
788 if (solver_control &&
789 solver_control->last_check() != SolverControl::success)
790 {
791 AssertThrow(false,
793 solver_control->last_value()));
794 }
795
796# else // DEAL_II_PETSC_WITH_MUMPS
797 Assert(
798 false,
800 "Your PETSc installation does not include a copy of "
801 "the MUMPS package necessary for this solver. You will need to configure "
802 "PETSc so that it includes MUMPS, recompile it, and then re-configure "
803 "and recompile deal.II as well."));
804
805 // Cast to void to silence compiler warnings
806 (void)A;
807 (void)x;
808 (void)b;
809# endif
810 }
811
812
813
814 void
816 {
817 symmetric_mode = flag;
818 }
819
820} // namespace PETScWrappers
821
823
824#endif // DEAL_II_WITH_PETSC
MPI_Comm get_mpi_communicator() const
SolverControl & control() const
void initialize(const PreconditionBase &preconditioner)
void perhaps_set_convergence_test() const
static PetscErrorCode convergence_test(KSP ksp, const PetscInt iteration, const PetscReal residual_norm, KSPConvergedReason *reason, void *solver_control)
SmartPointer< SolverControl, SolverBase > solver_control
void set_prefix(const std::string &prefix)
void initialize_ksp_with_comm(const MPI_Comm comm)
void solve(const MatrixBase &A, VectorBase &x, const VectorBase &b, const PreconditionBase &preconditioner)
virtual void set_solver_type(KSP &ksp) const
SolverBiCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
virtual void set_solver_type(KSP &ksp) const override
virtual void set_solver_type(KSP &ksp) const override
SolverBicgstab(SolverControl &cn, const AdditionalData &data=AdditionalData())
virtual void set_solver_type(KSP &ksp) const override
SolverCGS(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
virtual void set_solver_type(KSP &ksp) const override
virtual void set_solver_type(KSP &ksp) const override
SolverCR(SolverControl &cn, const AdditionalData &data=AdditionalData())
virtual void set_solver_type(KSP &ksp) const override
SolverChebychev(SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
virtual void set_solver_type(KSP &ksp) const override
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
virtual void set_solver_type(KSP &ksp) const override
SolverLSQR(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverPreOnly(SolverControl &cn, const AdditionalData &data=AdditionalData())
virtual void set_solver_type(KSP &ksp) const override
virtual void set_solver_type(KSP &ksp) const override
const AdditionalData additional_data
SolverRichardson(SolverControl &cn, const AdditionalData &data=AdditionalData())
virtual void set_solver_type(KSP &ksp) const override
SolverTCQMR(SolverControl &cn, const AdditionalData &data=AdditionalData())
virtual void set_solver_type(KSP &ksp) const override
SolverTFQMR(SolverControl &cn, const AdditionalData &data=AdditionalData())
void solve(const MatrixBase &A, VectorBase &x, const VectorBase &b)
virtual void set_solver_type(KSP &ksp) const override
void set_symmetric_mode(const bool flag)
SparseDirectMUMPS(SolverControl &cn, const AdditionalData &data=AdditionalData())
@ success
Stop iteration, goal reached.
#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 & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define AssertPETSc(code)
AdditionalData(const unsigned int restart_parameter=30, const bool right_preconditioning=false)
const MPI_Comm comm