PetscSolver.h 5.12 KB
Newer Older
1 2 3 4 5 6 7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9 10 11 12 13 14 15 16 17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19 20 21 22 23 24 25 26
 ******************************************************************************/


/** \file PetscSolver.h */

#ifndef AMDIS_PETSC_SOLVER_SEQ_H
#define AMDIS_PETSC_SOLVER_SEQ_H

27
#ifdef HAVE_SEQ_PETSC
28

29
#include "solver/LinearSolver.h"
30
#include "solver/PetscTypes.h"
31
#include "solver/MatrixStreams.h"
Praetorius, Simon's avatar
Praetorius, Simon committed
32
#include "utility/PetscWrapper.h"
33 34 35 36 37 38 39
#include "Timer.h"
#include <vector>
#include <iostream>
#include <boost/mpl/bool.hpp>

#include <petsc.h>
#include <petscksp.h>
40
#include <petscmat.h>
41 42 43 44 45
#include <petscvec.h>
#include <petscsys.h>
#include <petscao.h>

namespace AMDiS {
46

47 48
  /**
   * \ingroup Solver
49
   *
50 51 52
   * \brief Common base class for wrappers to use Petsc preconditioners in AMDiS.
   */
  template< class MatrixType, class VectorType >
53
  struct PetscPreconditionerInterface : public PreconditionerInterface
54
  {
55
    PetscPreconditionerInterface(std::string prefix, std::string name)
56
      : prefix(prefix), name(name) {}
57

58
    virtual ~PetscPreconditionerInterface() {}
59 60

    virtual void init(PC pc, const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix)
61
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
62
      petsc::options_insert_string(("-" + prefix + "pc_type " + name).c_str());
63 64
      PCSetFromOptions(pc);
    }
65

66
    virtual void exit() {}
67

68 69 70 71 72
  protected:
    std::string prefix;
    std::string name;
  };

73

74 75
  /// PETSc preconditioners using MATSEQAIJ and VECSEQ data structures
  typedef PetscPreconditionerInterface<PetscMatrix, PetscVector> PetscPreconditioner;
76

77 78
  /// PETSc preconditioner using MATNEST and VECNEST data structures
  typedef PetscPreconditionerInterface<PetscMatrixNested, PetscVectorNested> PetscPreconditionerNested;
79 80


81
  /// Runner that can be passed to LinearSolver to redirect the solution to PETSc
82
  template<typename MatrixType, typename VectorType>
83
  class PetscRunner : public RunnerInterface
84 85 86
  {
  public:
    /// Constructor of standard PETSc runner. Reads ksp and pc parameters from initfile.
87
    PetscRunner(LinearSolverInterface* oem);
88

89
    /// Destructor, deletes preconditioner object \ref preconditioner
90
    ~PetscRunner()
91 92 93 94 95 96
    {
      if (preconditioner) {
	delete preconditioner;
	preconditioner = NULL;
      }
    }
97

98 99
    /// Initialize the solver \ref ksp and preconditioner \ref pc
    void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix);
100

101 102
    /// Solve the linear equation \f$ A\cdot x = b \f$ by applying the PETSc solver \ref ksp
    int solve(const MatrixType& A, VectorType& x, const VectorType& b);
103

104
    /// Destroy solver \ref ksp and preconditioner \ref pc
105 106
    virtual void exit()
    {
107
      preconditioner->exit();
108 109 110
      KSPDestroy(&ksp);
    }

111
    /// Get the PETSc solver \ref ksp
112 113 114
    KSP getSolver()
    {
      return ksp;
115 116
    }

117
    /// Get the PETSc preconditioner \ref pc
118 119 120
    PC getPc()
    {
      return pc;
121 122
    }

123
    /// Returns preconditioner object \ref preconditioner
124
    PreconditionerInterface* getLeftPrecon()
125 126 127 128 129
    {
      return preconditioner;
    }

    /// Returns preconditioner object \ref preconditioner
130
    PreconditionerInterface* getRightPrecon()
131 132 133
    {
      return preconditioner;
    }
134 135


136
    static void createSubSolver(KSP &ksp_, Mat m, std::string kspPrefix_);
137

138
    static void setSolver(KSP ksp_, std::string kspPrefix_,
139
			  KSPType kspType, PCType pcType,
140 141 142
			  PetscReal rtol = PETSC_DEFAULT,
			  PetscReal atol = PETSC_DEFAULT,
			  PetscInt maxIt = PETSC_DEFAULT);
143 144 145


  protected:
146 147
    /// Creates a preconditioner object \ref preconditioner
    void setPrecon();
148

149
  protected:
150
    LinearSolverInterface& oem;
151

152 153 154 155 156 157 158
    /// PETSc solver object
    KSP ksp;

    /// PETSc preconditioner object
    PC pc;

    /// KSP database prefix
159
    std::string kspPrefix;
160

161 162 163
    bool zeroStartVector;	///< initialize startsolution for iterative solver with 0
    bool initialized;		///< true when \ref init() was called
    bool matSolverPackage;	///< true when using a direct solver
164

165 166
    /// Preconditioner object, that sets parameters to \ref pc
    PetscPreconditionerInterface<MatrixType, VectorType>* preconditioner;
167 168
  };

169

170 171
  /**
   * \ingroup Solver
172
   *
173 174 175 176
   * \brief
   * Wrapper for the external PETSc solver:
   *   http://www.mcs.anl.gov/petsc/
   *
177
   * This is a suite of data structures and routines for the
178 179
   * scalable (parallel) solution of scientific applications.
   */
180 181
  typedef LinearSolver< PetscMatrix, PetscVector, PetscRunner<PetscMatrix, PetscVector> > PetscSolver;
  typedef LinearSolver< PetscMatrixNested, PetscVectorNested, PetscRunner<PetscMatrixNested, PetscVectorNested> > PetscSolverNested;
182

183 184
} // end namespace AMDiS

185
#endif // HAVE_SEQ_PETSC
186

187 188
#include "solver/PetscSolver.hh"

189
#endif // AMDIS_PETSC_SOLVER_SEQ_H