PetscRunner.hpp 10.3 KB
Newer Older
1
2
3
4
5
6
7
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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#pragma once

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

#include <dune/common/unused.hh>

#include <amdis/Initfile.hpp>
#include <amdis/linearalgebra/RunnerInterface.hpp>
#include <amdis/linearalgebra/SolverInfo.hpp>
#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`.
   **/
  template <class Traits>
  class PetscRunner
      : public RunnerInterface<Traits>
  {
  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.
     **/
    explicit PetscRunner(std::string const& prefix)
      : prefix_(prefix)
    {
      Parameters::get(prefix + "->info", info_);
    }

    ~PetscRunner()
    {
      exit();
    }

    /// Implements \ref RunnerInterface::init()
71
    void init(typename Traits::Mat const& mat, typename Traits::Comm const& comm) override
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
    {
      exit();
#if HAVE_MPI
      KSPCreate(Environment::comm(), &ksp_);
#else
      KSPCreate(PETSC_COMM_SELF, &ksp_);
#endif

      PETSc::KSPSetOperators(ksp_, mat, mat);
      initKSP(ksp_, prefix_);
      initialized_ = true;
    }

    /// Implements \ref RunnerInterface::exit()
    void exit() override
    {
      if (initialized_)
        KSPDestroy(&ksp_);
      initialized_ = false;
    }

    /// Implements \ref RunnerInterface::solve()
    int solve(typename Traits::Mat const& A, typename Traits::Vec& x, typename Traits::Vec const& b, SolverInfo& solverInfo) override
    {
      assert(initialized_);
      DUNE_UNUSED_PARAMETER(A);

      KSPSolve(ksp_, b, x);

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

      KSPConvergedReason reason;
      KSPGetConvergedReason(ksp_, &reason);
      return reason > 0 ? 0 : int(reason);
    }

    KSP ksp() { return ksp_; }

  protected:
    // initialize the KSP solver from the initfile
113
    virtual void initKSP(KSP ksp, std::string prefix) const
114
115
    {
      // see https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPType.html
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
      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);
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
170
171
172
      }

      // 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
173
    virtual void initDirectSolver(KSP ksp, std::string prefix) const
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
    {
      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
198
    virtual void initPC(PC pc, std::string prefix) const
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
237
238
239
    {
      // 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
240
    virtual void initKSPParameters(KSP ksp, char const* ksptype, std::string prefix) const
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
291
292
293
    {
      // 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());
      }
    }

294
  protected:
295
296
297
298
299
300
301
302
    std::string prefix_;
    int info_ = 0;

    KSP ksp_;
    bool initialized_ = false;
  };

} // end namespace AMDiS