Chombo + EB + MF  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 
74 
75  virtual void define( LinearOp<T> *, bool a_homogeneous = true ); // pure virtual in base class
76 
77 
78 
79  ///
80 
81 #ifdef CH_USE_PETSC
82 
83  /**
84  Solve
85  */
86  virtual void solve( T& a_phi, const T& a_rhs );
87 #endif
88 
89 
90 
91  int solve_private(T& a_phi, const T& a_rhs);
92  ///
93  /**
94  Solve a_op*a_phi = a_rhs using the PETSC matrix free functions
95  The preconditioner (for which a matrix is formed) need not be
96  the same as the actual operator.
97  */
98 
99 #ifdef CH_USE_PETSC
100 
101  virtual void solve_mfree( T& a_phi,
102  const T& a_rhs,
103  LinearOp<T> *a_op );
104 #endif
105 
106  int solve_mfree_private(T& a_phi,
107  const T& a_rhs,
108  LinearOp<T> *a_op );
109 
114 
115 #ifdef CH_USE_PETSC
116  static PetscErrorCode apply_mfree(Mat A, Vec x, Vec f);
117 #endif
118  ///
119 
120  /**
121  viewResidual
122  */
124 
125 #ifdef CH_USE_PETSC
126 
127  /**
128  apply
129  */
130  virtual int applyOp(T & a_phi, const T & a_rhs );
131 #endif
132 
133 
134 
135  /**
136  set initial guess non-zero
137  */
138  void setInitialGuessNonzero( bool b = true )
139  {
140  m_nz_init_guess = b;
141  }
142 
143  /**
144  set initial guess non-zero
145  */
146  void setOptionsPrefix( const char prefix[] )
147  {
148 #ifdef CH_USE_PETSC
149  if (m_ksp)
150  {
151  KSPSetOptionsPrefix(m_ksp,prefix);
152  }
153  else
154  {
155  strcpy(m_prestring, prefix); // should check len
156  }
157 #endif
158  }
159 #ifdef CH_USE_PETSC
160  KSP getKSP() { return m_ksp; }
161 #endif
162 
163  /**
164  get an estimate of the number of nnz/row for matrix allocation
165  */
166  virtual int getNNZPerRow() const
167  {
168  return 9;
169  }
170  /**
171  do I support having formMatrix precompute exact NNZ/row
172  */
173  virtual bool supportNNZExact() const
174  {
175  return false;
176  }
177 
178  /*
179  an interface to apply operations to the rhs
180  */
181  virtual void rhsOperation(const T& a_rhs)
182  {
183  }
184  ///
185  /**
186  public member data: whether or not to use inhomogeneous boundary conditions.
187  */
189  ///
190  /**
191  member data: grid spacing
192  */
193 public:
195 
196  ///
197  /**
198  handling boundary conditions, turn it into a term to be added to the diag term
199  this function coordinates with addBCrhsValue for Dirichlet BC
200  for Neumann BC no RHS modification is required
201  */
202  virtual Real addBCdiagValue(const IntVect& a_iv, const IntVect& a_jv, const T& a_rhs, const DataIndex& a_datInd, const Real coeff = 1)
203  {
204  return 0;
205  }
206 public:
208  {
209 #ifdef CH_USE_PETSC
210  if (m_defined == 2)
211  {
212  m_defined = 1;
213  }
214 #endif
215  return 0;
216  }
217 
218  ///
219  /**
220  Infinity norm
221  */
222  Real normInfinity( const T& a_phi ) const;
223 
224 protected:
225  ///
226  /**
227  handling boundary conditions, turn it into a term to be added to the RHS
228  this function coordinates with addBCdiagValue for Dirichlet BC, should return a term that is to be added to RHS
229  */
230  virtual Real addBCrhsValue(const IntVect& a_iv, const T& a_phi, const DataIndex& a_datInd, const Real& coeff = 1)
231  {
232  return 0.0;
233  }
234 
235  /* petsc data */
236 #ifdef CH_USE_PETSC
237 public:
238  Mat m_mat;
239  void *m_ctx; // pointer for nonlnear solver call backs
240 
241  Vec m_xx, m_rr, m_bb;
242  SNES m_snes;
243  KSP m_ksp;
244 protected:
245  PetscInt m_defined;
246  PetscErrorCode (*m_function)(SNES,Vec,Vec,void*);
247 #if PETSC_VERSION_GE(3,5,0)
248  PetscErrorCode (*m_jacobian)(SNES,Vec,Mat,Mat,void*);
249 #else
250  PetscErrorCode (*m_jacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
251 #endif
252 #endif
253 protected:
254  char m_prestring[32];
255  bool m_null;
258 public:
259  int setup_solver( const T& a_phi );
260  int create_mat_vec( const T& a_phi );
261 #ifdef CH_USE_PETSC
263  PetscInt m_gid0;
264  PetscErrorCode putPetscInChombo( T& a_phi, const Vec xx );
265  PetscErrorCode putChomboInPetsc( Vec out, const T& a_phi );
266  virtual int formMatrix( Mat, const T* =0, PetscInt my0=0, PetscInt nloc=0, PetscInt *d=0, PetscInt *o=0) = 0;
267 #endif
268  virtual BaseFab<Real>& getRegFab(T& a_fab, const DataIndex& a_datInd) = 0;
269  const virtual BaseFab<Real>& getRegFab(const T& a_fab, const DataIndex& a_datInd) const = 0;
270  const virtual BaseFab<Real>& getRegFab(const T& a_fab, const DataIndex& a_datInd, Box& a_box) const = 0;
271  void setNull( bool n = true );
272 };
273 
274 ////////////////////////////////////////////////////////////////////////
275 // an instantiation of PetscSolver with FArrayBox
276 ////////////////////////////////////////////////////////////////////////
277 template <class T>
278 class PetscSolverFAB : public PetscSolver<T>
279 {
280 public:
282  {
283  }
284 
286  {
287  BaseFab<Real>& fabref = static_cast<BaseFab<Real>& >(a_fab[a_datInd]);
288  return fabref;
289  }
290 
291  const BaseFab<Real>& getRegFab(const LevelData<FArrayBox>& a_fab, const DataIndex& a_datInd) const
292  {
293  const BaseFab<Real>& fabref = static_cast<const BaseFab<Real>& >(a_fab[a_datInd]);
294  return fabref;
295  }
296 
297  const BaseFab<Real>& getRegFab(const LevelData<FArrayBox>& a_fab, const DataIndex& a_datInd, Box& a_box) const
298  {
299  const BaseFab<Real>& fabref = static_cast<const BaseFab<Real>& >(a_fab[a_datInd]);
300  a_box = fabref.box();
301  return fabref;
302  }
303 };
304 
305 ////////////////////////////////////////////////////////////////////////
306 // an instantiation of PetscSolver with Poisson
307 // solve: alpha u - div( beta grad u) = f
308 // - alpha and beta are constants
309 ////////////////////////////////////////////////////////////////////////
310 template <class T>
312 {
313 public:
315  virtual void define( Real a_dx, bool a_homogeneous = true );
316  virtual void define( LinearOp<T> *a_op, bool a_homogeneous = true )
317  {
318  PetscSolver<T>::define(a_op,a_homogeneous);
319  }
320 #ifdef CH_USE_PETSC
321  virtual int formMatrix(Mat, const LevelData<FArrayBox>* =0,PetscInt my0=0, PetscInt nloc=0, PetscInt *d=0, PetscInt *o=0 );
322 #endif
325 };
326 
327 #ifdef CH_USE_PETSC
328 ////////////////////////////////////////////////////////////////////////
329 // an instantiation of PetscSolver with viscous tensor Poisson
330 // solve: alpha a(x) + beta div( eta(x)grad(u) + I * lambda(x) div(u) ) = f
331 // - alpha and beta are constants, a, eta and lambda are fields
332 ////////////////////////////////////////////////////////////////////////
333 template <class T>
334 class PetscSolverViscousTensor : public PetscSolverFAB<LevelData<FArrayBox> >
335 {
336 public:
338  virtual void define( LinearOp<T>* a_operator, bool a_homogeneous = true );
339  /**
340  get an estimate of the number of nnz/row for matrix allocation
341  */
342  virtual int getNNZPerRow() const
343  {
344  return 13;
345  }
346 #ifdef CH_USE_PETSC
347  virtual int formMatrix(Mat, const LevelData<FArrayBox>* =0, PetscInt my0=0, PetscInt nloc=0, PetscInt *d=0, PetscInt *o=0 );
348 #endif
349 protected:
355 public:
357  {
358  m_alpha = alpha;
359  m_beta = beta;
360  m_a = a;
361  m_eta = eta;
362  m_lamb = lam;
363  }
364 };
365 #endif
366 
367 #include "NamespaceFooter.H"
368 
369 #ifdef CH_USE_PETSC
370 #ifndef CH_EXPLICIT_TEMPLATES
371 #include "PetscSolverI.H"
372 #endif // CH_EXPLICIT_TEMPLATES
373 #endif
374 
375 #endif /*_PETSCSOLVER_H_*/
std::ostream & pout()
Use this in place of std::cout for program output.
virtual void define(Real a_dx, bool a_homogeneous=true)
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:202
const BaseFab< Real > & getRegFab(const LevelData< FArrayBox > &a_fab, const DataIndex &a_datInd, Box &a_box) const
Definition: PetscSolver.H:297
KSP m_ksp
Definition: PetscSolver.H:243
LevelData< FluxBox > * m_eta
Definition: PetscSolver.H:353
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:412
Real m_beta
Definition: PetscSolver.H:350
PetscInt m_gid0
Definition: PetscSolver.H:263
PetscErrorCode(* m_jacobian)(SNES, Vec, Mat *, Mat *, MatStructure *, void *)
Definition: PetscSolver.H:250
virtual void rhsOperation(const T &a_rhs)
Definition: PetscSolver.H:181
Definition: PetscSolver.H:278
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 formMatrix(Mat, const LevelData< FArrayBox > *=0, PetscInt my0=0, PetscInt nloc=0, PetscInt *d=0, PetscInt *o=0)
virtual Real addBCrhsValue(const IntVect &a_iv, const T &a_phi, const DataIndex &a_datInd, const Real &coeff=1)
Definition: PetscSolver.H:230
virtual int getNNZPerRow() const
Definition: PetscSolver.H:166
void define(Real alpha, Real beta, LevelData< FArrayBox > *a, LevelData< FluxBox > *eta, LevelData< FluxBox > *lam)
Definition: PetscSolver.H:356
virtual int applyOp(T &a_phi, const T &a_rhs)
Definition: PetscSolverI.H:818
int resetOperator()
Definition: PetscSolver.H:207
virtual ~PetscSolver()
Definition: PetscSolverI.H:74
virtual BaseFab< Real > & getRegFab(T &a_fab, const DataIndex &a_datInd)=0
Real m_dx
Definition: PetscSolver.H:194
virtual bool supportNNZExact() const
Definition: PetscSolver.H:173
bool m_null
Definition: PetscSolver.H:255
virtual int getNNZPerRow() const
Definition: PetscSolver.H:342
Real m_alpha
Definition: PetscSolver.H:350
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:113
Real computeResidual()
Definition: PetscSolverI.H:774
PetscSolverFAB()
Definition: PetscSolver.H:281
PetscErrorCode putPetscInChombo(T &a_phi, const Vec xx)
Definition: PetscSolverI.H:694
void setNull(bool n=true)
Definition: PetscSolverI.H:99
LevelData< BaseFab< bool > > m_bccode
Definition: PetscSolver.H:257
Real m_alpha
Definition: PetscSolver.H:323
LinearOp< T > * m_op_mfree
Definition: PetscSolver.H:110
void * m_ctx
Definition: PetscSolver.H:239
int solve_mfree_private(T &a_phi, const T &a_rhs, LinearOp< T > *a_op)
Definition: PetscSolverI.H:539
new code
Definition: BoxLayoutData.H:170
Mat m_mat
Definition: PetscSolver.H:238
void setOptionsPrefix(const char prefix[])
Definition: PetscSolver.H:146
Vec m_bb
Definition: PetscSolver.H:241
LevelData< FluxBox > * m_lamb
Definition: PetscSolver.H:354
double Real
Definition: REAL.H:33
virtual void define(LinearOp< T > *a_operator, bool a_homogeneous=true)
LevelData< BaseFab< PetscInt > > m_gids
Definition: PetscSolver.H:262
void setInitialGuessNonzero(bool b=true)
Definition: PetscSolver.H:138
Real m_beta
Definition: PetscSolver.H:324
Real normInfinity(const T &a_phi) const
Definition: PetscSolverI.H:802
PetscSolver()
Definition: PetscSolverI.H:26
T m_phi_mfree
Definition: PetscSolver.H:111
virtual void define(LinearOp< T > *a_op, bool a_homogeneous=true)
Definition: PetscSolver.H:316
bool m_nz_init_guess
Definition: PetscSolver.H:256
char m_prestring[32]
Definition: PetscSolver.H:254
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:469
static PetscErrorCode apply_mfree(Mat A, Vec x, Vec f)
Definition: PetscSolverI.H:524
Real m_dxCrse
Definition: PetscSolver.H:351
Definition: LinearSolver.H:156
Vec m_xx
Definition: PetscSolver.H:241
const BaseFab< Real > & getRegFab(const LevelData< FArrayBox > &a_fab, const DataIndex &a_datInd) const
Definition: PetscSolver.H:291
Definition: DataIndex.H:114
T m_Lphi_mfree
Definition: PetscSolver.H:112
PetscErrorCode(* m_function)(SNES, Vec, Vec, void *)
Definition: PetscSolver.H:246
int setup_solver(const T &a_phi)
Definition: PetscSolverI.H:259
PetscErrorCode putChomboInPetsc(Vec out, const T &a_phi)
Definition: PetscSolverI.H:728
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
Definition: PetscSolver.H:36
void destroy()
Definition: PetscSolverI.H:43
Definition: FArrayBox.H:45
KSP getKSP()
Definition: PetscSolver.H:160
SNES m_snes
Definition: PetscSolver.H:242
Definition: PetscSolver.H:311
BaseFab< Real > & getRegFab(LevelData< FArrayBox > &a_fab, const DataIndex &a_datInd)
Definition: PetscSolver.H:285
virtual int formMatrix(Mat, const LevelData< FArrayBox > *=0, PetscInt my0=0, PetscInt nloc=0, PetscInt *d=0, PetscInt *o=0)
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:245
Definition: PetscSolver.H:334
Vec m_rr
Definition: PetscSolver.H:241
LevelData< FArrayBox > * m_a
Definition: PetscSolver.H:352
bool m_homogeneous
Definition: PetscSolver.H:188
const Box & box() const
Definition: BaseFabImplem.H:271