Chombo + EB  3.2
PetscSolver.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 #ifndef _PETSCSOLVER_H_
12 #define _PETSCSOLVER_H_
13 
14 #include "LinearSolver.H"
15 #include "parstream.H"
16 #include "CH_Timer.H"
17 #include "LevelData.H"
18 
19 #ifdef CH_USE_PETSC
20 #include "petsc.h"
21 #include "petscmat.h"
22 #include "petscksp.h"
23 #include "petscviewer.h"
24 #endif
25 
26 #include "NamespaceHeader.H"
27 
28 ///
29 /**
30  Interface to PETSC linear solvers
31  BC implementation hooks: addBCdiagValue (formMatrix) and addBCrhsValue (solveprivate)
32  m_bccode should be encoded true for cells that are at the domain boundaries at formMatrix time.
33  */
34 ///
35 template <class T>
36 class PetscSolver : public LinearSolver<T>
37 {
38 public:
39  PetscSolver();
40  virtual ~PetscSolver();
41  void destroy();
42  ///
43  /**
44  Set whether the solver uses only homogeneous boundary conditions
45  */
46  virtual void setHomogeneous(bool a_homogeneous)
47  {
48  m_homogeneous = a_homogeneous;
49  }
50  /**
51  Set Function F(u) = 0 and Jacobian dF(u)/du for nonlinear solves
52  */
53 #ifdef CH_USE_PETSC
54  virtual void setFunctionAndJacobian( PetscErrorCode (*f)(SNES,Vec,Vec,void*),
55 #if PETSC_VERSION_GE(3,5,0)
56  PetscErrorCode (*j)(SNES,Vec,Mat,Mat,void*)
57 #else
58  PetscErrorCode (*j)(SNES,Vec,Mat*,Mat*,MatStructure*,void*)
59 #endif
60  )
61  {
62  m_function = f;
63  m_jacobian = j;
64  }
65  static PetscErrorCode ksp_monitor_pout(KSP ksp, PetscInt it, PetscReal rnorm ,void *ctx)
66  {
67  pout() << " KSP:: iteration = " << it << " residual norm = " << rnorm << std::endl;
68  return 0;
69  }
70 #endif
71  //
72  virtual void define( Real a_dx, bool a_homogeneous = true ); // convenience define
73  virtual void define( LinearOp<T> *, bool a_homogeneous = true ); // pure virtual in base class
74  ///
75  /**
76  Solve
77  */
78  virtual void solve( T& a_phi, const T& a_rhs );
79  int solve_private(T& a_phi, const T& a_rhs);
80  ///
81  /**
82  Solve a_op*a_phi = a_rhs using the PETSC matrix free functions
83  The preconditioner (for which a matrix is formed) need not be
84  the same as the actual operator.
85  */
86 
87  virtual void solve_mfree( T& a_phi,
88  const T& a_rhs,
89  LinearOp<T> *a_op );
90 
91  int solve_mfree_private(T& a_phi,
92  const T& a_rhs,
93  LinearOp<T> *a_op );
94 
99 
100 #ifdef CH_USE_PETSC
101  static PetscErrorCode apply_mfree(Mat A, Vec x, Vec f);
102 #endif
103  ///
104 
105  /**
106  viewResidual
107  */
109  /**
110  apply
111  */
112  virtual int applyOp(T & a_phi, const T & a_rhs );
113  /**
114  set initial guess non-zero
115  */
116  void setInitialGuessNonzero( bool b = true )
117  {
118  m_nz_init_guess = b;
119  }
120 
121  /**
122  set initial guess non-zero
123  */
124  void setOptionsPrefix( const char prefix[] )
125  {
126 #ifdef CH_USE_PETSC
127  if (m_ksp)
128  {
129  KSPSetOptionsPrefix(m_ksp,prefix);
130  }
131  else
132  {
133  strcpy(m_prestring, prefix); // should check len
134  }
135 #endif
136  }
137 #ifdef CH_USE_PETSC
138  KSP getKSP() { return m_ksp; }
139 #endif
140 
141  /**
142  get an estimate of the number of nnz/row for matrix allocation
143  */
144  virtual int getNNZPerRow() const
145  {
146  return 9;
147  }
148  /**
149  do I support having formMatrix precompute exact NNZ/row
150  */
151  virtual bool supportNNZExact() const
152  {
153  return false;
154  }
155 
156  /*
157  an interface to apply operations to the rhs
158  */
159  virtual void rhsOperation(const T& a_rhs)
160  {
161  }
162  ///
163  /**
164  public member data: whether or not to use inhomogeneous boundary conditions.
165  */
167  ///
168  /**
169  member data: grid spacing
170  */
171 public:
173 
174  ///
175  /**
176  handling boundary conditions, turn it into a term to be added to the diag term
177  this function coordinates with addBCrhsValue for Dirichlet BC
178  for Neumann BC no RHS modification is required
179  */
180  virtual Real addBCdiagValue(const IntVect& a_iv, const IntVect& a_jv, const T& a_rhs, const DataIndex& a_datInd, const Real coeff = 1)
181  {
182  return 0;
183  }
184 public:
186  {
187 #ifdef CH_USE_PETSC
188  if (m_defined == 2)
189  {
190  m_defined = 1;
191  }
192 #endif
193  return 0;
194  }
195 
196  ///
197  /**
198  Infinity norm
199  */
200  Real normInfinity( const T& a_phi ) const;
201 
202 protected:
203  ///
204  /**
205  handling boundary conditions, turn it into a term to be added to the RHS
206  this function coordinates with addBCdiagValue for Dirichlet BC, should return a term that is to be added to RHS
207  */
208  virtual Real addBCrhsValue(const IntVect& a_iv, const T& a_phi, const DataIndex& a_datInd, const Real& coeff = 1)
209  {
210  return 0.0;
211  }
212 
213  /* petsc data */
214 #ifdef CH_USE_PETSC
215 public:
216  Mat m_mat;
217  void *m_ctx; // pointer for nonlnear solver call backs
218 
219  Vec m_xx, m_rr, m_bb;
220  SNES m_snes;
221  KSP m_ksp;
222 protected:
223  PetscInt m_defined;
224  PetscErrorCode (*m_function)(SNES,Vec,Vec,void*);
225 #if PETSC_VERSION_GE(3,5,0)
226  PetscErrorCode (*m_jacobian)(SNES,Vec,Mat,Mat,void*);
227 #else
228  PetscErrorCode (*m_jacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
229 #endif
230 #endif
231 protected:
232  char m_prestring[32];
233  bool m_null;
236 public:
237  int setup_solver( const T& a_phi );
238  int create_mat_vec( const T& a_phi );
239 #ifdef CH_USE_PETSC
241  PetscInt m_gid0;
242  PetscErrorCode putPetscInChombo( T& a_phi, const Vec xx );
243  PetscErrorCode putChomboInPetsc( Vec out, const T& a_phi );
244  virtual int formMatrix( Mat, const T* =0, PetscInt my0=0, PetscInt nloc=0, PetscInt *d=0, PetscInt *o=0) = 0;
245 #endif
246  virtual BaseFab<Real>& getRegFab(T& a_fab, const DataIndex& a_datInd) = 0;
247  const virtual BaseFab<Real>& getRegFab(const T& a_fab, const DataIndex& a_datInd) const = 0;
248  const virtual BaseFab<Real>& getRegFab(const T& a_fab, const DataIndex& a_datInd, Box& a_box) const = 0;
249  void setNull( bool n = true );
250 };
251 
252 ////////////////////////////////////////////////////////////////////////
253 // an instantiation of PetscSolver with FArrayBox
254 ////////////////////////////////////////////////////////////////////////
255 template <class T>
256 class PetscSolverFAB : public PetscSolver<T>
257 {
258 public:
260  {
261  }
262 
264  {
265  BaseFab<Real>& fabref = static_cast<BaseFab<Real>& >(a_fab[a_datInd]);
266  return fabref;
267  }
268 
269  const BaseFab<Real>& getRegFab(const LevelData<FArrayBox>& a_fab, const DataIndex& a_datInd) const
270  {
271  const BaseFab<Real>& fabref = static_cast<const BaseFab<Real>& >(a_fab[a_datInd]);
272  return fabref;
273  }
274 
275  const BaseFab<Real>& getRegFab(const LevelData<FArrayBox>& a_fab, const DataIndex& a_datInd, Box& a_box) const
276  {
277  const BaseFab<Real>& fabref = static_cast<const BaseFab<Real>& >(a_fab[a_datInd]);
278  a_box = fabref.box();
279  return fabref;
280  }
281 };
282 
283 ////////////////////////////////////////////////////////////////////////
284 // an instantiation of PetscSolver with Poisson
285 // solve: alpha u - div( beta grad u) = f
286 // - alpha and beta are constants
287 ////////////////////////////////////////////////////////////////////////
288 template <class T>
290 {
291 public:
293  virtual void define( Real a_dx, bool a_homogeneous = true );
294  virtual void define( LinearOp<T> *a_op, bool a_homogeneous = true )
295  {
296  PetscSolver<T>::define(a_op,a_homogeneous);
297  }
298 #ifdef CH_USE_PETSC
299  virtual int formMatrix(Mat, const LevelData<FArrayBox>* =0,PetscInt my0=0, PetscInt nloc=0, PetscInt *d=0, PetscInt *o=0 );
300 #endif
303 };
304 
305 ////////////////////////////////////////////////////////////////////////
306 // an instantiation of PetscSolver with viscous tensor Poisson
307 // solve: alpha a(x) + beta div( eta(x)grad(u) + I * lambda(x) div(u) ) = f
308 // - alpha and beta are constants, a, eta and lambda are fields
309 ////////////////////////////////////////////////////////////////////////
310 template <class T>
311 class PetscSolverViscousTensor : public PetscSolverFAB<LevelData<FArrayBox> >
312 {
313 public:
315  virtual void define( LinearOp<T>* a_operator, bool a_homogeneous = true );
316  /**
317  get an estimate of the number of nnz/row for matrix allocation
318  */
319  virtual int getNNZPerRow() const
320  {
321  return 13;
322  }
323 #ifdef CH_USE_PETSC
324  virtual int formMatrix(Mat, const LevelData<FArrayBox>* =0, PetscInt my0=0, PetscInt nloc=0, PetscInt *d=0, PetscInt *o=0 );
325 #endif
326 protected:
327  Real m_alpha, m_beta;
332 public:
334  {
335  m_alpha = alpha;
336  m_beta = beta;
337  m_a = a;
338  m_eta = eta;
339  m_lamb = lam;
340  }
341 };
342 
343 #include "NamespaceFooter.H"
344 
345 #ifdef CH_USE_PETSC
346 #ifndef CH_EXPLICIT_TEMPLATES
347 #include "PetscSolverI.H"
348 #endif // CH_EXPLICIT_TEMPLATES
349 #endif
350 
351 #endif /*_PETSCSOLVER_H_*/
std::ostream & pout()
Use this in place of std::cout for program output.
Definition: LinearSolver.H:28
virtual void define(Real a_dx, bool a_homogeneous=true)
Definition: PetscSolverI.H:81
virtual Real addBCdiagValue(const IntVect &a_iv, const IntVect &a_jv, const T &a_rhs, const DataIndex &a_datInd, const Real coeff=1)
Definition: PetscSolver.H:180
KSP m_ksp
Definition: PetscSolver.H:221
LevelData< FluxBox > * m_eta
Definition: PetscSolver.H:330
virtual void solve(T &a_phi, const T &a_rhs)
Definition: PetscSolverI.H:106
int solve_private(T &a_phi, const T &a_rhs)
Definition: PetscSolverI.H:420
Real m_beta
Definition: PetscSolver.H:327
PetscInt m_gid0
Definition: PetscSolver.H:241
virtual void rhsOperation(const T &a_rhs)
Definition: PetscSolver.H:159
Definition: PetscSolver.H:256
static PetscErrorCode ksp_monitor_pout(KSP ksp, PetscInt it, PetscReal rnorm, void *ctx)
Definition: PetscSolver.H:65
virtual int formMatrix(Mat, const T *=0, PetscInt my0=0, PetscInt nloc=0, PetscInt *d=0, PetscInt *o=0)=0
virtual void solve_mfree(T &a_phi, const T &a_rhs, LinearOp< T > *a_op)
Definition: PetscSolverI.H:116
virtual void setHomogeneous(bool a_homogeneous)
Definition: PetscSolver.H:46
virtual int getNNZPerRow() const
Definition: PetscSolver.H:319
const Box & box() const
Definition: BaseFabImplem.H:271
virtual Real addBCrhsValue(const IntVect &a_iv, const T &a_phi, const DataIndex &a_datInd, const Real &coeff=1)
Definition: PetscSolver.H:208
const BaseFab< Real > & getRegFab(const LevelData< FArrayBox > &a_fab, const DataIndex &a_datInd, Box &a_box) const
Definition: PetscSolver.H:275
void define(Real alpha, Real beta, LevelData< FArrayBox > *a, LevelData< FluxBox > *eta, LevelData< FluxBox > *lam)
Definition: PetscSolver.H:333
virtual int applyOp(T &a_phi, const T &a_rhs)
Definition: PetscSolverI.H:826
int resetOperator()
Definition: PetscSolver.H:185
virtual ~PetscSolver()
Definition: PetscSolverI.H:74
virtual BaseFab< Real > & getRegFab(T &a_fab, const DataIndex &a_datInd)=0
Real m_dx
Definition: PetscSolver.H:172
virtual bool supportNNZExact() const
Definition: PetscSolver.H:151
bool m_null
Definition: PetscSolver.H:233
void const char const int const int const int const Real const Real * A
Definition: Lapack.H:83
bool m_mfree_homogeneous
Definition: PetscSolver.H:98
Real computeResidual()
Definition: PetscSolverI.H:782
PetscSolverFAB()
Definition: PetscSolver.H:259
PetscErrorCode putPetscInChombo(T &a_phi, const Vec xx)
Definition: PetscSolverI.H:702
void setNull(bool n=true)
Definition: PetscSolverI.H:99
LevelData< BaseFab< bool > > m_bccode
Definition: PetscSolver.H:235
Real m_alpha
Definition: PetscSolver.H:301
LinearOp< T > * m_op_mfree
Definition: PetscSolver.H:95
void * m_ctx
Definition: PetscSolver.H:217
int solve_mfree_private(T &a_phi, const T &a_rhs, LinearOp< T > *a_op)
Definition: PetscSolverI.H:547
new code
Definition: BoxLayoutData.H:170
Mat m_mat
Definition: PetscSolver.H:216
void setOptionsPrefix(const char prefix[])
Definition: PetscSolver.H:124
Vec m_bb
Definition: PetscSolver.H:219
LevelData< FluxBox > * m_lamb
Definition: PetscSolver.H:331
PetscErrorCode(* m_function)(SNES, Vec, Vec, void *)
Definition: PetscSolver.H:224
double Real
Definition: REAL.H:33
LevelData< BaseFab< PetscInt > > m_gids
Definition: PetscSolver.H:240
void setInitialGuessNonzero(bool b=true)
Definition: PetscSolver.H:116
Real m_beta
Definition: PetscSolver.H:302
virtual int getNNZPerRow() const
Definition: PetscSolver.H:144
Real normInfinity(const T &a_phi) const
Definition: PetscSolverI.H:810
PetscSolver()
Definition: PetscSolverI.H:26
T m_phi_mfree
Definition: PetscSolver.H:96
virtual void define(LinearOp< T > *a_op, bool a_homogeneous=true)
Definition: PetscSolver.H:294
bool m_nz_init_guess
Definition: PetscSolver.H:234
char m_prestring[32]
Definition: PetscSolver.H:232
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:465
static PetscErrorCode apply_mfree(Mat A, Vec x, Vec f)
Definition: PetscSolverI.H:532
Real m_dxCrse
Definition: PetscSolver.H:328
Definition: LinearSolver.H:156
Vec m_xx
Definition: PetscSolver.H:219
Definition: DataIndex.H:112
T m_Lphi_mfree
Definition: PetscSolver.H:97
int setup_solver(const T &a_phi)
Definition: PetscSolverI.H:267
PetscErrorCode putChomboInPetsc(Vec out, const T &a_phi)
Definition: PetscSolverI.H:736
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
Definition: PetscSolver.H:36
void destroy()
Definition: PetscSolverI.H:43
PetscErrorCode(* m_jacobian)(SNES, Vec, Mat *, Mat *, MatStructure *, void *)
Definition: PetscSolver.H:228
Definition: FArrayBox.H:45
KSP getKSP()
Definition: PetscSolver.H:138
SNES m_snes
Definition: PetscSolver.H:220
Definition: PetscSolver.H:289
BaseFab< Real > & getRegFab(LevelData< FArrayBox > &a_fab, const DataIndex &a_datInd)
Definition: PetscSolver.H:263
virtual void setFunctionAndJacobian(PetscErrorCode(*f)(SNES, Vec, Vec, void *), PetscErrorCode(*j)(SNES, Vec, Mat *, Mat *, MatStructure *, void *))
Definition: PetscSolver.H:54
int create_mat_vec(const T &a_phi)
Definition: PetscSolverI.H:131
PetscInt m_defined
Definition: PetscSolver.H:223
Definition: PetscSolver.H:311
const BaseFab< Real > & getRegFab(const LevelData< FArrayBox > &a_fab, const DataIndex &a_datInd) const
Definition: PetscSolver.H:269
Vec m_rr
Definition: PetscSolver.H:219
LevelData< FArrayBox > * m_a
Definition: PetscSolver.H:329
bool m_homogeneous
Definition: PetscSolver.H:166