PetscSolver.h 4.88 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
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * 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.
 * 
 ******************************************************************************/


/** \file PetscSolver.h */

#ifndef AMDIS_PETSC_SOLVER_SEQ_H
#define AMDIS_PETSC_SOLVER_SEQ_H

27
#ifdef HAVE_SEQ_PETSC
28

29
#include "solver/MTL4Solver.h"
30
#include "solver/PetscTypes.h"
31
#include "solver/MatrixStreams.h"
32 33 34 35 36 37 38 39 40 41 42 43 44
#include "Timer.h"
#include <vector>
#include <iostream>
#include <boost/mpl/bool.hpp>

#include <petsc.h>
#include <petscksp.h>
#include <petscmat.h> 
#include <petscvec.h>
#include <petscsys.h>
#include <petscao.h>

namespace AMDiS {
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 71 72 73
    
  /**
   * \ingroup Solver
   * 
   * \brief Common base class for wrappers to use Petsc preconditioners in AMDiS.
   */
  template< class MatrixType, class VectorType >
  struct Petsc_BasePreconditioner : public OEMPreconditioner
  {
    Petsc_BasePreconditioner(std::string prefix, std::string name) : prefix(prefix), name(name) {}
    virtual ~Petsc_BasePreconditioner() {}
    
    virtual void init(PC pc, const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) 
    {
      PetscOptionsInsertString(("-" + prefix + "pc_type " + name).c_str());
      PCSetFromOptions(pc);
    }
    
    virtual void exit() {}
    
  protected:
    std::string prefix;
    std::string name;
  };

  
  typedef Petsc_BasePreconditioner<PetscMatrix, PetscVector> PetscPreconditioner;
  typedef Petsc_BasePreconditioner<PetscMatrixNested, PetscVectorNested> PetscPreconditionerNested;
  
74
  
75
  template<typename MatrixType, typename VectorType>
76 77 78 79
  class PetscRunner : public OEMRunner
  {
  public:
    /// Constructor of standard PETSc runner. Reads ksp and pc parameters from initfile.
80 81 82 83 84 85 86 87 88 89
    PetscRunner(LinearSolver* oem);
    
    /// Destructor, deletes preconditioner object \ref preconditioner
    ~PetscRunner() 
    {
      if (preconditioner) {
	delete preconditioner;
	preconditioner = NULL;
      }
    }
90
    
91 92
    /// Initialize the solver \ref ksp and preconditioner \ref pc
    void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix);
93
      
94 95
    /// 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);
96
    
97
    /// Destroy solver \ref ksp and preconditioner \ref pc
98 99
    virtual void exit()
    {
100
      preconditioner->exit();
101 102 103
      KSPDestroy(&ksp);
    }

104
    /// Get the PETSc solver \ref ksp
105 106 107 108 109
    KSP getSolver() 
    { 
      return ksp; 
    }

110
    /// Get the PETSc preconditioner \ref pc
111 112 113 114 115
    PC getPc() 
    { 
      return pc; 
    }

116 117 118 119 120 121 122 123 124 125 126 127
    /// Returns preconditioner object \ref preconditioner
    OEMPreconditioner* getLeftPrecon()
    {
      return preconditioner;
    }

    /// Returns preconditioner object \ref preconditioner
    OEMPreconditioner* getRightPrecon()
    {
      return preconditioner;
    }
    
128 129 130 131 132 133 134 135
    
    static void createSubSolver(KSP &ksp_, Mat m, std::string kspPrefix_);
    
    static void setSolver(KSP ksp_, std::string kspPrefix_,
			  KSPType kspType, PCType pcType, 
			  PetscReal rtol = PETSC_DEFAULT,
			  PetscReal atol = PETSC_DEFAULT,
			  PetscInt maxIt = PETSC_DEFAULT);
136 137 138 139 140 141 142
			  
    
  protected:    
    /// Creates a preconditioner object \ref preconditioner
    void setPrecon();
    
  protected:
143 144 145 146 147 148 149 150 151
    LinearSolver& oem;
    
    /// PETSc solver object
    KSP ksp;

    /// PETSc preconditioner object
    PC pc;

    /// KSP database prefix
152
    std::string kspPrefix;
153
    
154 155 156 157 158 159
    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
    
    /// Preconditioner object, that set parameters to \ref pc
    Petsc_BasePreconditioner<MatrixType, VectorType>* preconditioner;
160 161 162 163 164 165 166 167 168 169 170 171 172
  };

  
  /**
   * \ingroup Solver
   * 
   * \brief
   * Wrapper for the external PETSc solver:
   *   http://www.mcs.anl.gov/petsc/
   *
   * This is a suite of data structures and routines for the 
   * scalable (parallel) solution of scientific applications.
   */
173 174
  typedef MTL4Solver< PetscMatrix, PetscVector, PetscRunner<PetscMatrix, PetscVector> > PetscSolver;
  typedef MTL4Solver< PetscMatrixNested, PetscVectorNested, PetscRunner<PetscMatrixNested, PetscVectorNested> > PetscSolverNested;
175 176 177
  
} // end namespace AMDiS

178
#endif // HAVE_SEQ_PETSC
179

180 181
#include "solver/PetscSolver.hh"

182
#endif // AMDIS_PETSC_SOLVER_SEQ_H