Chombo + EB  3.0
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 // BVS, MA, Nov 18, 2009
12 
13 #ifndef _PETSCSOLVER_H_
14 #define _PETSCSOLVER_H_
15 
16 #include "LinearSolver.H"
17 #include "parstream.H"
18 #include "CH_Timer.H"
19 
20 
21 #ifdef CH_USE_PETSC
22 #include "petsc.h"
23 #include "petscmat.h"
24 #include "petscksp.h"
25 #include "petscviewer.h"
26 #endif
27 
28 #include "NamespaceHeader.H"
29 
30 ///
31 /**
32  Interface to PETSC linear solvers
33  BC implementation hooks: addBCdiagValue (formMatrix) and addBCrhsValue (solveprivate)
34  m_bccode should be encoded true for cells that are at the domain boundaries at formMatrix time.
35  */
36 ///
37 template <class T>
38 class PetscSolver : public LinearSolver<T>
39 {
40 public:
41  PetscSolver();
42  virtual ~PetscSolver();
43  void destroy();
44  ///
45  /**
46  Set whether the solver uses only homogeneous boundary conditions
47  */
48  virtual void setHomogeneous(bool a_homogeneous)
49  {
50  m_homogeneous = a_homogeneous;
51  }
52  /**
53  Set Function F(u) = 0 and Jacobian dF(u)/du for nonlinear solves
54  */
55 #ifdef CH_USE_PETSC
56  virtual void setFunctionAndJacobian( PetscErrorCode (*f)(SNES,Vec,Vec,void*),
57  PetscErrorCode (*j)(SNES,Vec,Mat*,Mat*,MatStructure*,void*)
58  )
59  {
60  m_function = f;
61  m_jacobian = j;
62  }
63 #endif
64  //
65  virtual void define(LinearOp<T>* a_operator, bool a_homogeneous = false);
66  ///
67  /**
68  Solve
69  */
70  virtual void solve( T& a_phi, const T& a_rhs );
71  int solve_private(T& a_phi, const T& a_rhs);
72  ///
73  /**
74  viewResidual
75  */
77  /**
78  apply
79  */
80  int applyOp(T & a_phi, const T & a_rhs );
81  /**
82  set initial guess non-zero
83  */
84  void setInitialGuessNonzero( bool b = true )
85  {
86  m_nz_init_guess = b;
87  }
88 
89  /**
90  get an estimate of the number of nnz/row for matrix allocation
91  */
92  virtual int getNNZPerRow() const
93  {
94  pout()<<"getNNZPerRow PetscSolver"<<endl;
95  return 9;
96  }
97  ///
98  /**
99  public member data: whether or not to use inhomogeneous boundary conditions.
100  */
102  ///
103  /**
104  member data: grid spacing
105  */
106 public:
108 protected:
109  ///
110  /**
111  handling boundary conditions, turn it into a term to be added to the diag term
112  this function coordinates with addBCrhsValue for Dirichlet BC
113  for Neumann BC no RHS modification is required
114  */
115  virtual Real addBCdiagValue(const IntVect& a_iv, const IntVect& a_jv, const T& a_rhs, DataIterator dit, const Real coeff = 1)
116  {
117  return 0;
118  }
119 
121  {
122 #ifdef CH_USE_PETSC
123  if (m_defined == 2)
124  {
125  m_defined = 1;
126  }
127 #endif
128  return 0;
129  }
130 
131  ///
132  /**
133  handling boundary conditions, turn it into a term to be added to the RHS
134  this function coordinates with addBCdiagValue for Dirichlet BC, should return a term that is to be added to RHS
135  */
136  virtual Real addBCrhsValue(const IntVect& a_iv, const T& a_phi, DataIterator dit, const Real& coeff = 1)
137  {
138  return 0.0;
139  }
140 
141  /* petsc data */
142 #ifdef CH_USE_PETSC
143  Mat m_mat;
144  void *m_ctx; // pointer for nonlnear solver call backs
145 protected:
146  Vec m_xx, m_rr, m_bb;
147  KSP m_ksp;
148  SNES m_snes;
149  PetscInt m_defined;
150  PetscErrorCode (*m_function)(SNES,Vec,Vec,void*);
151  PetscErrorCode (*m_jacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
152 #endif
153 protected:
154  bool m_null;
158  int m_gid0;
159 
160  int create_mat_vec( const T& a_phi );
161  int setup_solver( const T& a_phi );
162 public:
163 #ifdef CH_USE_PETSC
164  int putPetscInChombo( T& a_phi, const Vec xx );
165  int putChomboInPetsc( Vec out, const T& a_phi );
166  virtual int formMatrix( Mat, const T& a_rhs ) = 0;
167 #endif
168  virtual BaseFab<Real>& getRegFab(T& a_fab, DataIterator& dit) = 0;
169  const virtual BaseFab<Real>& getRegFab(const T& a_fab, DataIterator& dit) const = 0;
170  const virtual BaseFab<Real>& getRegFab(const T& a_fab, DataIterator& dit, Box& a_box) const = 0;
171  virtual void defineData(T& a_data, const T& a_template) = 0;
172  void setNull( bool n = true );
173 };
174 
175 ////////////////////////////////////////////////////////////////////////
176 // an instantialtion of PetscSolver with FArrayBox
177 // this looks like terrible C++, no reason for a template here -- mfa
178 ////////////////////////////////////////////////////////////////////////
179 template <class T>
180 class PetscSolverFAB : public PetscSolver<T>
181 {
182 public:
184  {
185  }
186 
188  {
189  BaseFab<Real>& fabref = static_cast<BaseFab<Real>& >(a_fab[dit()]);
190  return fabref;
191  }
192 
193  const BaseFab<Real>& getRegFab(const LevelData<FArrayBox>& a_fab, DataIterator& dit) const
194  {
195  const BaseFab<Real>& fabref = static_cast<const BaseFab<Real>& >(a_fab[dit()]);
196  return fabref;
197  }
198 
199  const BaseFab<Real>& getRegFab(const LevelData<FArrayBox>& a_fab, DataIterator& a_dit, Box& a_box) const
200  {
201  const BaseFab<Real>& fabref = static_cast<const BaseFab<Real>& >(a_fab[a_dit()]);
202  a_box = fabref.box();
203  return fabref;
204  }
205 
207  {
208  IntVect idghosts = a_phi.ghostVect();
209  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
210  a_fab.define( dbl, 1, idghosts );
211  }
212 };
213 
214 ////////////////////////////////////////////////////////////////////////
215 // an instantialtion of PetscSolver with Poisson
216 // solve: alpha u - div( beta grad u) = f
217 // - alpha and beta are constants
218 ////////////////////////////////////////////////////////////////////////
219 template <class T>
221 {
222 public:
224  virtual void define(LinearOp<T>* a_operator, bool a_homogeneous = false);
225 #ifdef CH_USE_PETSC
226  virtual int formMatrix( Mat, const LevelData<FArrayBox>& a_rhs );
227 #endif
230 };
231 
232 ////////////////////////////////////////////////////////////////////////
233 // an instantialtion of PetscSolver with viscouse tensor Poisson
234 // solve: alpha a(x) + beta div( eta(x)grad(u) + I * lambda(x) div(u) ) = f
235 // - alpha and beta are constants, a, eta and lambda are fields
236 ////////////////////////////////////////////////////////////////////////
237 template <class T>
238 class PetscSolverViscousTensor : public PetscSolverFAB<LevelData<FArrayBox> >
239 {
240 public:
242  virtual void define(LinearOp<LevelData<FArrayBox> >* a_operator, bool a_homogeneous = false);
243  /**
244  get an estimate of the number of nnz/row for matrix allocation
245  */
246  virtual int getNNZPerRow() const
247  {
248  return 13;
249  }
250 #ifdef CH_USE_PETSC
251  virtual int formMatrix( Mat, const LevelData<FArrayBox>& a_rhs );
252 #endif
253 protected:
259 public:
261  {
262  m_alpha = alpha;
263  m_beta = beta;
264  m_a = a;
265  m_eta = eta;
266  m_lamb = lam;
267  }
268 };
269 
270 #include "NamespaceFooter.H"
271 
272 #ifdef CH_USE_PETSC
273 #ifndef CH_EXPLICIT_TEMPLATES
274 #include "PetscSolverI.H"
275 #endif // CH_EXPLICIT_TEMPLATES
276 #endif
277 
278 #endif /*_PETSCSOLVER_H_*/
std::ostream & pout()
Use this in place of std::cout for program output.
Definition: LinearSolver.H:28
LevelData< FluxBox > * m_eta
Definition: PetscSolver.H:257
virtual void solve(T &a_phi, const T &a_rhs)
Definition: PetscSolverI.H:95
int solve_private(T &a_phi, const T &a_rhs)
Definition: PetscSolverI.H:297
virtual Real addBCdiagValue(const IntVect &a_iv, const IntVect &a_jv, const T &a_rhs, DataIterator dit, const Real coeff=1)
Definition: PetscSolver.H:115
Real m_beta
Definition: PetscSolver.H:254
const BaseFab< Real > & getRegFab(const LevelData< FArrayBox > &a_fab, DataIterator &dit) const
Definition: PetscSolver.H:193
Definition: PetscSolver.H:180
int m_gid0
Definition: PetscSolver.H:158
virtual void define(LinearOp< T > *a_operator, bool a_homogeneous=false)
Definition: PetscSolverI.H:78
Definition: DataIterator.H:140
BaseFab< Real > & getRegFab(LevelData< FArrayBox > &a_fab, DataIterator &dit)
Definition: PetscSolver.H:187
virtual void setHomogeneous(bool a_homogeneous)
Definition: PetscSolver.H:48
virtual int getNNZPerRow() const
Definition: PetscSolver.H:92
int applyOp(T &a_phi, const T &a_rhs)
Definition: PetscSolverI.H:507
int resetOperator()
Definition: PetscSolver.H:120
virtual ~PetscSolver()
Definition: PetscSolverI.H:71
Real m_dx
Definition: PetscSolver.H:107
bool m_null
Definition: PetscSolver.H:154
virtual int getNNZPerRow() const
Definition: PetscSolver.H:246
Real m_alpha
Definition: PetscSolver.H:254
Real computeResidual()
Definition: PetscSolverI.H:478
PetscSolverFAB()
Definition: PetscSolver.H:183
void setNull(bool n=true)
Definition: PetscSolverI.H:88
LevelData< BaseFab< bool > > m_bccode
Definition: PetscSolver.H:157
Real m_alpha
Definition: PetscSolver.H:228
Definition: BoxLayoutData.H:136
virtual void define(LinearOp< T > *a_operator, bool a_homogeneous=false)
virtual void defineData(T &a_data, const T &a_template)=0
LevelData< FluxBox > * m_lamb
Definition: PetscSolver.H:258
double Real
Definition: REAL.H:33
virtual void define(const DisjointBoxLayout &dp, int comps, const IntVect &ghost=IntVect::Zero, const DataFactory< T > &a_factory=DefaultDataFactory< T >())
Definition: LevelDataI.H:70
virtual BaseFab< Real > & getRegFab(T &a_fab, DataIterator &dit)=0
const DisjointBoxLayout & disjointBoxLayout() const
Definition: LevelData.H:196
void setInitialGuessNonzero(bool b=true)
Definition: PetscSolver.H:84
A BoxLayout that has a concept of disjointedness.
Definition: DisjointBoxLayout.H:31
Real m_beta
Definition: PetscSolver.H:229
void defineData(LevelData< FArrayBox > &a_fab, const LevelData< FArrayBox > &a_phi)
Definition: PetscSolver.H:206
const BaseFab< Real > & getRegFab(const LevelData< FArrayBox > &a_fab, DataIterator &a_dit, Box &a_box) const
Definition: PetscSolver.H:199
PetscSolver()
Definition: PetscSolverI.H:28
virtual Real addBCrhsValue(const IntVect &a_iv, const T &a_phi, DataIterator dit, const Real &coeff=1)
Definition: PetscSolver.H:136
bool m_nz_init_guess
Definition: PetscSolver.H:155
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:465
Real m_dxCrse
Definition: PetscSolver.H:255
Definition: LinearSolver.H:158
int setup_solver(const T &a_phi)
Definition: PetscSolverI.H:193
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
Definition: PetscSolver.H:38
void destroy()
Definition: PetscSolverI.H:44
Definition: FArrayBox.H:44
virtual void define(LinearOp< LevelData< FArrayBox > > *a_operator, bool a_homogeneous=false)
Definition: PetscSolver.H:220
int create_mat_vec(const T &a_phi)
Definition: PetscSolverI.H:109
Definition: PetscSolver.H:238
LevelData< FArrayBox > * m_a
Definition: PetscSolver.H:256
bool m_homogeneous
Definition: PetscSolver.H:101
void setVTParams(Real alpha, Real beta, LevelData< FArrayBox > *a, LevelData< FluxBox > *eta, LevelData< FluxBox > *lam)
Definition: PetscSolver.H:260
T m_gids
Definition: PetscSolver.H:156
const IntVect & ghostVect() const
Definition: LevelData.H:157
const Box & box() const
Definition: BaseFabImplem.H:202