PetscSolver.h 4.79 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
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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90

#include "solver/LinearSolver.h"
#include "solver/MatrixStreams.h"
#include "solver/PetscTypes.h"
#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 {
  
  class PetscRunner : public OEMRunner
  {
  public:
    /// Constructor of standard PETSc runner. Reads ksp and pc parameters from initfile.
    PetscRunner(LinearSolver* oem_);

    typedef SolverMatrix<Matrix<DOFMatrix*> > BlockMatrix;
    
    /// initialize the solver \ref ksp and preconditioner \ref pc
    void init(const BlockMatrix& A, const Mat& fullMatrix);
      
    /// solve the linear equation \f$ A\cdot x = b \f$ by applying the PETSc solver \ref ksp
    int solve(const Mat& A, Vec& x, const Vec& b);
    
    /// destroy solver \ref ksp and preconditioner \ref pc
    virtual void exit()
    {
      KSPDestroy(&ksp);
    }

    /// get the PETSc solver \ref ksp
    KSP getSolver() 
    { 
      return ksp; 
    }

    /// get the PETSc preconditioner \ref pc
    PC getPc() 
    { 
      return pc; 
    }

    ~PetscRunner() { }
    
  protected:    
    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);

    LinearSolver& oem;
    
91
  //private:    
92
93
94
95
96
97
98
    /// PETSc solver object
    KSP ksp;

    /// PETSc preconditioner object
    PC pc;

    /// KSP database prefix
99
    std::string kspPrefix;
100
101
    
    bool zeroStartVector;
102
    bool initialized;
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
  };

  
  /**
   * \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.
   */
  template< typename Runner = PetscRunner >
  class PetscSolver : public LinearSolver 
  {
  public:
    /// Creator class used in the LinearSolverMap.
    class Creator : public LinearSolverCreator
    {
    public:
      virtual ~Creator() {}
      
      /// Returns a new PetscSolver object.
      LinearSolver* create() 
      { 
	return new PetscSolver<Runner>(this->name); 
      }
    };
    
    PetscSolver(std::string n)
    : LinearSolver(n),
      runner(this)
    {}
    
    ~PetscSolver() {}
    
    virtual OEMRunner* getRunner()
    {
      return &runner;
    }
    
145
    void setNestedVectors(bool nested = true)
146
    {
147
148
      vecSol.isNested = nested;
      vecRhs.isNested = nested;
149
150
    }
 
151
    void setNestedMatrix(bool nested = true)
152
    {
153
      petscMat.isNested=nested;
154
155
156
157
    } 

    void setNested(bool n)
    {
158
159
      setNestedVectors(n); 
      setNestedMatrix(n); 
160
    }
161
    
162
163
164
165
166
167
168
169
170
171
172
173
174
175
  protected:    
    Runner runner;

    int solveLinearSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
			  SystemVector& x, 
			  SystemVector& b,
			  bool createMatrixData,
			  bool storeMatrixData) 
    { FUNCNAME("PetscSolver::solveLinearSystem()");

      Timer t;
      // transfer matrix and rhs-vector to PETSc data-structures
      if (createMatrixData) {
	petscMat << A;
176
	runner.init(A, petscMat.matrix);
177
178
179
180
181
182
183
184
      }

      vecSol << x;
      vecRhs << b;
      INFO(info, 8)("fill PETSc matrix needed %.5f seconds\n", t.elapsed());

      // solve the linear system using PETSc solvers
      t.reset();
185
      error = runner.solve(petscMat.matrix, vecSol.vector, vecRhs.vector);
186
187
188
189
190
191

      // transfer solution from PETSc vector to SystemVector
      vecSol.vector >> x;

      vecSol.destroy();
      vecRhs.destroy();
192
193
      
      if (!storeMatrixData) {
194
	petscMat.destroy();
195
196
	runner.exit();
      }
197
198
199
200
201
202
203
204
205
206
207
208
209
      return error;
    }
    
  private:
    /// PETSc System-Matrix
    PetscMatrix petscMat;

    /// Solution and RHS vectors.
    PetscVector vecSol, vecRhs;
  };
  
} // end namespace AMDiS

210
#endif // HAVE_SEQ_PETSC
211
212

#endif // AMDIS_PETSC_SOLVER_SEQ_H