LinearSolver.hpp 10.2 KB
Newer Older
1
2
3
4
5
6
#pragma once

#include <petscpc.h>
#include <petscksp.h>

#include <amdis/Initfile.hpp>
7
#include <amdis/linearalgebra/LinearSolverInterface.hpp>
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
#include <amdis/linearalgebra/petsc/Interface.hpp>

#ifndef KSPDGGMRES
#define KSPDGGMRES "dggmres"
#endif


namespace AMDiS
{
  /// \brief Wrapper around PETSc KSP and PC objects to solve a linear system
  /**
   * Configure the KSP and PC types using initfile parameters, based on an initfile prefix `[prefix]`:
   * All parameters are names as the PETSc command-line parameters with underscore replaced by space
   * and removing the *type* postfix, e.g. instead of `ksp_type` write `ksp`, instead of `mat_solver_type`
   * write `mat solver`.
   *
   * **Examples:**
   * ```
   * [prefix]->max it: 100    % maximum solver iterations
   * [prefix]->rtol: 1.e-10   % relative residual tolerance
   * [prefix]->pc: bjacobi    % preconditioner type
   * [prefix]->pc->sub ksp: preonly   % sub KSP type used on the blocks
   * [prefix]->pc->sub ksp->pc: ilu   % preconditioner of the sub KSP type
   * ```
   *
   * For the configuration of sub PC types and sub KSP types, one can use the current parameter
   * as new prefix, e.g. `[sub-prefix] := [prefix]->pc`, or `[sub-sub-prefix] := [prefix]->pc->sub ksp`.
   * Those sub PC and sub KSP types can be configured with all the possible parameters.
   *
   * For the configuration using command-line arguments possible a ksp-prefix or pc-prefix must be
   * assigned to distinguish different KSP and PC objects. Therefore, for each initfile prefix, a
   * PETSc prefix can be set: `[prefix]->prefix: name`, where `name` can be used on the commandline, e.g.
   * `-name_ksp_type fgmres`.
   **/
42
43
44
  template <class Matrix, class VectorX, class VectorY = VectorX>
  class LinearSolver
      : public LinearSolverInterface<Matrix,VectorX,VectorY>
45
  {
46
    using Vector = VectorX;
Praetorius, Simon's avatar
Praetorius, Simon committed
47

48
49
50
51
52
53
54
55
56
57
  public:
    /// Constructor.
    /**
     * Stores the initfile prefix to configure the ksp and pc solver.
     *
     * Reads an info flag `[prefix]->info`:
     * 0 ... No solver convergence information
     * 1 ... Minimal output, print residual every 10th iteration
     * 2 ... Full convergence output. Print residual, true residual and rel. residual every iteration.
     **/
58
    LinearSolver(std::string const& /*name*/, std::string const& prefix)
59
60
61
62
63
      : prefix_(prefix)
    {
      Parameters::get(prefix + "->info", info_);
    }

64
    ~LinearSolver()
65
    {
66
      finish();
67
68
    }

69
70
    /// Implements \ref LinearSolverInterface::init()
    void init(Matrix const& A) override
71
    {
72
      finish();
73
74
75
76
77
78
#if HAVE_MPI
      KSPCreate(Environment::comm(), &ksp_);
#else
      KSPCreate(PETSC_COMM_SELF, &ksp_);
#endif

79
      PETSc::KSPSetOperators(ksp_, A.matrix(), A.matrix());
80
81
82
83
      initKSP(ksp_, prefix_);
      initialized_ = true;
    }

84
85
    /// Implements \ref LinearSolverInterface::finish()
    void finish() override
86
87
88
89
90
91
    {
      if (initialized_)
        KSPDestroy(&ksp_);
      initialized_ = false;
    }

92
    /// Implements \ref Dune::InverOperator::apply()
Praetorius, Simon's avatar
Praetorius, Simon committed
93
    void apply(Vector& x, Vector const& b, Dune::InverseOperatorResult& stat) override
94
    {
95
      KSPSolve(ksp_, b.vector(), x.vector());
96
97
98
99
100
101

      if (info_ > 0)
        KSPReasonView(ksp_, PETSC_VIEWER_STDOUT_WORLD);

      KSPConvergedReason reason;
      KSPGetConvergedReason(ksp_, &reason);
102
103

      stat.converged = (reason > 0);
104
105
106
107
108
109
    }

    KSP ksp() { return ksp_; }

  protected:
    // initialize the KSP solver from the initfile
110
    virtual void initKSP(KSP ksp, std::string prefix) const
111
112
    {
      // see https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPType.html
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
      auto kspType = Parameters::get<std::string>(prefix + "->ksp");
      std::string kspTypeStr = kspType.value_or("default");

      if (!kspType)
        Parameters::get(prefix, kspTypeStr);

      if (kspTypeStr == "direct")
        initDirectSolver(ksp, prefix);
      else if (kspTypeStr != "default") {
        KSPSetType(ksp, kspTypeStr.c_str());

        // initialize some KSP specific parameters
        initKSPParameters(ksp, kspTypeStr.c_str(), prefix);

        // set initial guess to nonzero only for non-preonly ksp type
        if (kspTypeStr != "preonly")
          KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
      }

      // set a KSPMonitor if info > 0
      int info = 0;
      Parameters::get(prefix + "->info", info);
      if (info == 1) {
        KSPMonitorSet(ksp, PETSc::KSPMonitorCyclic, nullptr, nullptr);
      } else if (info > 1) {
        KSPMonitorSet(ksp, PETSc::KSPMonitorNoisy, nullptr, nullptr);
      }

      // set solver tolerances
      auto maxIt = Parameters::get<PetscInt>(prefix + "->max it");
      auto rTol = Parameters::get<PetscReal>(prefix + "->rtol");
      auto aTol = Parameters::get<PetscReal>(prefix + "->atol");
      auto dTol = Parameters::get<PetscReal>(prefix + "->divtol");
      if (maxIt || rTol || aTol || dTol)
        KSPSetTolerances(ksp,
          rTol.value_or(PETSC_DEFAULT), aTol.value_or(PETSC_DEFAULT),
          dTol.value_or(PETSC_DEFAULT), maxIt.value_or(PETSC_DEFAULT));

      auto kspPrefixParam = Parameters::get<std::string>(prefix + "->prefix");
      if (kspPrefixParam) {
        std::string kspPrefix = kspPrefixParam.value() + "_";
        KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
        KSPSetFromOptions(ksp);
      }

      PC pc;
      KSPGetPC(ksp, &pc);
      initPC(pc, prefix + "->pc");

      // show details of ksp object
      if (info >= 10) {
        KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD);
      }
    }


    // initialize a direct solver from the initfile
170
    virtual void initDirectSolver(KSP ksp, std::string prefix) const
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
    {
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
      KSPSetType(ksp, KSPRICHARDSON);

      PC pc;
      KSPGetPC(ksp, &pc);
      PCSetType(pc, PCLU);
      auto matSolverType = Parameters::get<std::string>(prefix + "->mat solver");
      if (matSolverType)
        PETSc::PCFactorSetMatSolverType(pc, matSolverType.value().c_str());
      else {
        Mat Amat, Pmat;
        KSPGetOperators(ksp, &Amat, &Pmat);
        MatType type;
        MatGetType(Pmat, &type);
        if (std::strcmp(type, MATSEQAIJ) == 0)
          PETSc::PCFactorSetMatSolverType(pc, MATSOLVERUMFPACK);
        else if (std::strcmp(type, MATAIJ) == 0)
          PETSc::PCFactorSetMatSolverType(pc, MATSOLVERMUMPS);
      }
    }


    // initialize the preconditioner pc from the initfile
195
    virtual void initPC(PC pc, std::string prefix) const
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
    {
      // see https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCType.html
      auto pcType = Parameters::get<std::string>(prefix);
      std::string pcTypeStr = "default";
      if (pcType) {
        pcTypeStr = pcType.value();
        PCSetType(pc, pcTypeStr.c_str());
      }

      if (pcTypeStr == "lu") {
        // configure matsolverpackage
        // see https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSolverType.html
        auto matSolverType = Parameters::get<std::string>(prefix + "->mat solver");
        if (matSolverType)
          PETSc::PCFactorSetMatSolverType(pc, matSolverType.value().c_str());
      } else if (pcTypeStr == "bjacobi") {
        // configure sub solver
        PetscInt n_local, first;
        KSP* ksps;
        PCSetUp(pc);
        PCBJacobiGetSubKSP(pc, &n_local, &first, &ksps);
        for (PetscInt i = 0; i < n_local; ++i)
          initKSP(ksps[i], prefix + "->sub ksp");
      } else if (pcTypeStr == "ksp") {
        // configure sub ksp type
        KSP ksp;
        PCSetUp(pc);
        PCKSPGetKSP(pc, &ksp);
        initKSP(ksp, prefix + "->ksp");
      }

      auto pcPrefixParam = Parameters::get<std::string>(prefix + "->prefix");
      if (pcPrefixParam) {
        std::string pcPrefix = pcPrefixParam.value() + "_";
        PCSetOptionsPrefix(pc, pcPrefix.c_str());
        PCSetFromOptions(pc);
      }
    }


    // provide initfile parameters for some PETSc KSP parameters
237
    virtual void initKSPParameters(KSP ksp, char const* ksptype, std::string prefix) const
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
    {
      // parameters for the Richardson solver
      if (std::strcmp(ksptype, KSPRICHARDSON) == 0)
      {
        auto scale = Parameters::get<PetscReal>(prefix + "->scale");
        if (scale)
          KSPRichardsonSetScale(ksp, scale.value());
      }
      // parameters for the gmres solver
      else if (std::strcmp(ksptype, KSPGMRES) == 0 ||
               std::strcmp(ksptype, KSPFGMRES) == 0 ||
               std::strcmp(ksptype, KSPLGMRES) == 0 ||
               std::strcmp(ksptype, KSPDGGMRES) == 0 ||
               std::strcmp(ksptype, KSPPGMRES) == 0)
      {
        auto restart = Parameters::get<PetscInt>(prefix + "->restart");
        if (restart)
          KSPGMRESSetRestart(ksp, restart.value());

        auto gramschmidt = Parameters::get<std::string>(prefix + "->gramschmidt");
        if (gramschmidt) {
          if (gramschmidt.value() == "classical") {
            KSPGMRESSetOrthogonalization(ksp, KSPGMRESClassicalGramSchmidtOrthogonalization);
            auto refinementType = Parameters::get<std::string>(prefix + "->refinement type");
            if (refinementType) {
              if (refinementType.value() == "never")
                KSPGMRESSetCGSRefinementType(ksp, KSP_GMRES_CGS_REFINE_NEVER);
              else if (refinementType.value() == "ifneeded")
                KSPGMRESSetCGSRefinementType(ksp, KSP_GMRES_CGS_REFINE_IFNEEDED);
              else if (refinementType.value() == "always")
                KSPGMRESSetCGSRefinementType(ksp, KSP_GMRES_CGS_REFINE_ALWAYS);
            }
          } else if (gramschmidt.value() == "modified") {
            KSPGMRESSetOrthogonalization(ksp, KSPGMRESModifiedGramSchmidtOrthogonalization);
          }
        }
      }
      // parameters for the BiCGStab_ell solver
      else if (std::strcmp(ksptype, KSPBCGSL) == 0)
      {
        auto ell = Parameters::get<PetscInt>(prefix + "->ell");
        if (ell)
           KSPBCGSLSetEll(ksp, ell.value());
      }
      // parameters for the GCR solver
      else if (std::strcmp(ksptype, KSPGCR) == 0)
      {
        auto restart = Parameters::get<PetscInt>(prefix + "->restart");
        if (restart)
           KSPGCRSetRestart(ksp, restart.value());
      }
    }

291
  protected:
292
293
294
295
296
297
298
299
    std::string prefix_;
    int info_ = 0;

    KSP ksp_;
    bool initialized_ = false;
  };

} // end namespace AMDiS