Chombo + EB  3.0
PetscSolverI.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 _PETSCSOLVERI_H_
12 #define _PETSCSOLVERI_H_
13 
14 #include "LevelData.H"
15 #include "FluxBox.H"
16 #include "BoxIterator.H"
17 #include "CH_Timer.H"
18 #include "memusage.H"
19 #include "NamespaceHeader.H"
20 
21 #include <private/kspimpl.h> /*I "petscksp.h" I*/
22 #include <private/pcimpl.h>
23 
24 // *******************************************************
25 // PetscSolver Implementation
26 // *******************************************************
27 template <class T>
29  :m_homogeneous(false),
30  m_mat(0), // m_xx, m_rr, m_bb;
31  m_ksp(0),
32  m_snes(0),
33  m_defined(0),
34  m_function(0),
35  m_jacobian(0),
36  m_null(false),
37  m_nz_init_guess(false),
38  m_gid0(0)
39 {
40  m_dx = 0.;
41 }
42 // *******************************************************
43 template <class T>
45 {
46  if ( m_defined )
47  {
48 #ifdef CH_USE_PETSC
49  MatDestroy(&m_mat);
50  VecDestroy(&m_bb);
51  VecDestroy(&m_xx);
52  VecDestroy(&m_rr);
53 #endif
54  m_defined = 0;
55  }
56 #ifdef CH_USE_PETSC
57  if ( m_ksp )
58  {
59  KSPDestroy( &m_ksp );
60  m_ksp = 0;
61  }
62  if ( m_snes )
63  {
64  SNESDestroy( &m_snes );
65  m_snes = 0;
66  }
67 #endif
68 }
69 // *******************************************************
70 template <class T>
72 {
73  destroy();
74 }
75 
76 // *******************************************************
77 template <class T>
79  bool a_homogeneous)
80 {
81  m_homogeneous = a_homogeneous; // not used!!!
82  m_dx = a_operator->dx();
83  CH_assert(m_dx!=0.);
84 }
85 
86 // *******************************************************
87 template <class T>
88 void PetscSolver<T>::setNull( bool n /*= true*/ )
89 {
90  m_null = n; CH_assert(0);
91 }
92 
93 // *******************************************************
94 template <class T>
95 void PetscSolver<T>::solve( T & a_phi,
96  const T & a_rhs )
97 {
98  T& phi = a_phi;
99  const T& rhs = a_rhs;
100  solve_private( phi, rhs );
101 }
102 
103 
104 // *******************************************************
105 // create_mat_vec
106 // Create 'm_mat', 'm_xx', .... Constructs 'm_gids'.
107 //
108 template <class T>
109 int PetscSolver<T>::create_mat_vec( const T& a_phi )
110 {
111  CH_TIME("PetscSolver::create_mat_vec");
112 #ifdef CH_USE_PETSC
113  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
114  DataIterator dit( dbl );
115  int ierr;
116  const int nc = a_phi.nComp();
117 #ifdef CH_MPI
118  MPI_Comm wcomm = Chombo_MPI::comm;
119 #else
120  MPI_Comm wcomm = PETSC_COMM_SELF;
121 #endif
122 
123  if ( !m_mat )
124  {
125  print_memory_line("Before AMG set up");
126  m_defined = 2;
127  // global ids with one ghost cells
128  defineData( m_gids, a_phi );
129  //m_gids.define( dbl, 1, idghosts );
130  // get first (zero based) id on this processor
131  // defineData(m_bccode,a_phi);
132  IntVect idghosts = a_phi.ghostVect();
133  if (idghosts == IntVect::Zero)
134  {
135  MayDay::Error("PetscSolver<T>::create_mat_vec: No ghost cells in input LevelData<>.");
136  }
137  m_bccode.define(dbl, 1, idghosts);
138 
139  int data = 0;
140  for ( dit = a_phi.dataIterator() ; dit.ok() ; ++dit )
141  {
142  const Box &box = dbl.get(dit());
143  if (CH_SPACEDIM==3) data += box.size(0)*box.size(1)*box.size(2);
144  else data += box.size(0)*box.size(1);
145  BaseFab<Real> &gidsFab = getRegFab(m_gids,dit);
146  //FArrayBox& gidsFab = m_gids[dit()];
147  gidsFab.setVal(-1.0);
148  }
149  const int NN = nc*(int)data;
150 #ifdef CH_MPI
151  int result;
152  ierr = MPI_Scan( &data, &result, 1, MPI_INT, MPI_SUM, wcomm );
153  m_gid0 = result - data;
154 #else
155  m_gid0 = 0;
156 #endif
157  Real gid = (Real)m_gid0;
158  for ( dit = a_phi.dataIterator() ; dit.ok() ; ++dit )
159  {
160  BaseFab<Real> &gidsFab = getRegFab(m_gids,dit);
161  const Box& box = dbl.get(dit());
162  BoxIterator bit(box);
163  for (bit.begin(); bit.ok(); bit.next(), gid++ )
164  {
165  IntVect iv = bit();
166  gidsFab(iv,0) = (Real)gid;
167  }
168  }
169  m_gids.exchange();
170  // data
171  int nnzrow = getNNZPerRow();
172  ierr = MatCreateMPIAIJ( wcomm, NN, NN, PETSC_DECIDE, PETSC_DECIDE,
173  nnzrow, PETSC_NULL, nnzrow/2, PETSC_NULL, &m_mat ); CHKERRQ(ierr);
174  ierr = MatSetFromOptions( m_mat ); CHKERRQ(ierr);
175 
176  // create vectors
177  ierr = VecCreate( wcomm, &m_bb ); CHKERRQ(ierr);
178  ierr = VecSetFromOptions( m_bb ); CHKERRQ(ierr);
179  ierr = VecSetSizes( m_bb, NN, PETSC_DECIDE ); CHKERRQ(ierr);
180  ierr = VecDuplicate( m_bb, &m_rr ); CHKERRQ(ierr);
181  ierr = VecDuplicate( m_bb, &m_xx ); CHKERRQ(ierr);
182  }
183 #endif
184 
185  return 0;
186 }
187 
188 // *******************************************************
189 // setup_solver
190 // - creates solver if needed. forms matrix, sets up KSP
191 //
192 template <class T>
193 int PetscSolver<T>::setup_solver( const T& a_phi )
194 {
195  CH_TIMERS("PetscSolver::setup_solver");
196  CH_TIMER("solve-setup-1st", t1);
197  CH_TIMER("solve-setup-rest", t2);
198 #ifdef CH_USE_PETSC
199  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
200  DataIterator dit( dbl );
201  int ierr;
202  KSP ksp;
203 #ifdef CH_MPI
204  MPI_Comm wcomm = Chombo_MPI::comm;
205 #else
206  MPI_Comm wcomm = PETSC_COMM_SELF;
207 #endif
208  if ( m_defined == 0 )
209  {
210  m_defined = 2;
211  print_memory_line("Before AMG set up");
212 
213  ierr = create_mat_vec( a_phi ); CHKERRQ(ierr);
214 
215  CH_START(t1);
216 
217  // Add values to A
218  ierr = formMatrix( m_mat, a_phi ); CHKERRQ(ierr);
219 
220  // create solvers
221  if ( m_function && !m_snes )
222  {
223  ierr = SNESCreate( wcomm, &m_snes ); CHKERRQ(ierr);
224  ierr = SNESSetFunction( m_snes, m_rr, m_function, (void*)&a_phi ); CHKERRQ(ierr);
225  ierr = SNESSetJacobian( m_snes, m_mat, m_mat, m_jacobian, (void*)&a_phi ); CHKERRQ(ierr);
226  ierr = SNESSetFromOptions( m_snes ); CHKERRQ(ierr);
227  ierr = SNESGetKSP( m_snes, &ksp ); CHKERRQ(ierr);
228  ierr = SNESSetApplicationContext( m_snes, (void*)this ); CHKERRQ(ierr);
229  }
230  else if ( !m_function && !m_ksp )
231  {
232  // create the KSP so that we can set KSP parameters
233  KSPCreate( wcomm, &m_ksp );
234  KSPSetFromOptions( m_ksp );
235  ierr = KSPSetOperators(m_ksp,m_mat,m_mat,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
236  ierr = KSPSetInitialGuessNonzero(m_ksp, m_nz_init_guess ? PETSC_TRUE : PETSC_FALSE ); CHKERRQ(ierr);
237  ksp = m_ksp;
238  }
239  else CH_assert(0);
240 
241  { // coordinates
242  PC pc; const PCType type; PetscInt sz,ind,bs,n,m;
243  ierr = KSPGetPC( ksp, &pc ); CHKERRQ(ierr);
244  ierr = PCGetType( pc, &type ); CHKERRQ(ierr);
245  ierr = MatGetBlockSize( m_mat, &bs ); CHKERRQ( ierr );
246  if ( strcmp(type,PCGAMG) == 0 && bs > 1 )
247  {
248  PetscReal *coords;
249  DataIterator dit(a_phi.disjointBoxLayout());
250 
251  ierr = MatGetLocalSize( m_mat, &m, &n ); CHKERRQ(ierr);
252  sz = CH_SPACEDIM*(m/bs);
253  ierr = PetscMalloc( (sz+1)*sizeof(PetscReal), &coords ); CHKERRQ(ierr);
254  for ( dit = a_phi.dataIterator(), ind = 0 ; dit.ok() ; ++dit )
255  {
256  const Box &box = a_phi.getBoxes()[dit];
257  BoxIterator bit(box);
258  for (bit.begin(); bit.ok(); bit.next())
259  {
260  IntVect iv = bit(); // coordinate in any scaled, shifted, rotated frame.
261  for (int k=0; k<CH_SPACEDIM; k++) coords[ind++] = (PetscReal)iv[k];
262  }
263  }
264  CH_assert(ind==sz);
265  ierr = PCSetCoordinates( pc, CH_SPACEDIM, coords ); CHKERRQ(ierr);
266  ierr = PetscFree( coords ); CHKERRQ(ierr);
267  }
268  }
269 
270  CH_STOP(t1);
271  print_memory_line("After AMG set up");
272  }
273  else if ( m_defined == 1 )
274  {
275  m_defined = 2;
276  // form A -- m_mat
277  CH_START(t2);
278  ierr = MatZeroEntries( m_mat ); CHKERRQ(ierr);
279  ierr = formMatrix( m_mat, a_phi ); CHKERRQ(ierr);
280  if ( m_ksp )
281  {
282  ksp = m_ksp;
283  }
284  else
285  {
286  ierr = SNESGetKSP( m_snes, &ksp ); CHKERRQ(ierr);
287  }
288  ierr = KSPSetOperators(m_ksp,m_mat,m_mat,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
289  CH_STOP(t2);
290  }
291 #endif
292  return 0;
293 }
294 
295 // *******************************************************
296 template <class T>
298  const T& a_rhs )
299 {
300  CH_TIMERS("PetscSolver::solve_private");
301  CH_TIMER("formRHS", t2);
302  CH_TIMER("solve", t3);
303  CH_TIMER("output", t4);
304 #ifdef CH_USE_PETSC
305  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
306  DataIterator dit( dbl );
307  int ierr;
308  const int nc = a_rhs.nComp();
309 #ifdef CH_MPI
310  MPI_Comm wcomm = Chombo_MPI::comm;
311 #else
312  MPI_Comm wcomm = PETSC_COMM_SELF;
313 #endif
314 
315  // setup if needed
316  ierr = setup_solver( a_phi ); CHKERRQ(ierr);
317 
318  CH_START(t2);
319  // form RHS -- m_bb
320  Real idx2 = 1.e0/(this->m_dx*this->m_dx);
321  Real addV;
322 
323  // add X and B from Chombo to PETSc and add stuff for EB to B
324  ierr = VecSet( m_xx, 0.); CHKERRQ(ierr);
325  ierr = VecSet( m_bb, 0.); CHKERRQ(ierr);
326  for ( dit = a_rhs.dataIterator() ; dit.ok() ; ++dit )
327  {
328  const BaseFab<Real> &rhsFab = getRegFab(a_rhs,dit);
329  const BaseFab<Real> &xFab = getRegFab(a_phi,dit);
330  const Box& box = dbl.get(dit());
331  // const EBCellFAB& ebfab = a_rhs[dit()];
332  // const EBISBox& ebbox = ebfab.getEBISBox();
333  //const Box &box = ebbox.getRegion(); // for Helmholtz
334  const BaseFab<Real> &gidsFab = getRegFab(m_gids,dit);
335  BoxIterator bit(box);
336  for (bit.begin(); bit.ok(); bit.next())
337  {
338  IntVect iv = bit();
339  int mm, ki = nc*(int)gidsFab(iv,0);
340  for (mm=0;mm<nc;mm++,ki++)
341  {
342  PetscScalar v = rhsFab(iv,mm);
343  //if (ki==2)pout() << " *************************** b1(2)=" << v << endl;
344  // v += addBCrhsValue(iv,a_phi,dit,idx2);
345  addV = addBCrhsValue(iv,a_phi,dit,idx2);
346  if (addV == BaseFabRealSetVal)
347  {
348  v = 0;
349  }
350  else
351  {
352  v += addV;
353  }
354  ierr = VecSetValues(m_bb,1,&ki,&v,INSERT_VALUES); CHKERRQ(ierr);
355  //if (ki==2)pout() << " *************************** b2(2)=" << v << endl;
356  //if ( nonZeroInit )
357  //{
358  v = xFab(iv,mm);
359  ierr = VecSetValues(m_xx,1,&ki,&v,INSERT_VALUES); CHKERRQ(ierr);
360  //}
361  }
362  }
363  }
364  ierr = VecAssemblyBegin( m_bb ); CHKERRQ(ierr);
365  ierr = VecAssemblyEnd( m_bb ); CHKERRQ(ierr);
366  //if (nonZeroInit)
367  //{
368  ierr = VecAssemblyBegin( m_xx ); CHKERRQ(ierr);
369  ierr = VecAssemblyEnd( m_xx ); CHKERRQ(ierr);
370  //}
371  CH_STOP(t2);
372 
373  // null space for periodic BCs
374  if ( m_null )
375  {
376  MatNullSpace nullsp;
377  ierr = MatNullSpaceCreate(wcomm,PETSC_TRUE,0,PETSC_NULL,&nullsp);CHKERRQ(ierr);
378  CH_assert(m_ksp); // not used yet, needs fix for nonlinear
379  ierr = KSPSetNullSpace( m_ksp, nullsp ); CHKERRQ(ierr);
380  }
381 
382  // solve
383  CH_START(t3);
384  if ( m_snes )
385  {
386  ierr = SNESSolve( m_snes, m_bb, m_xx ); CHKERRQ(ierr);
387 
388  }
389  else
390  {
391  ierr = KSPSolve( m_ksp, m_bb, m_xx ); CHKERRQ(ierr);
392  }
393  CH_STOP(t3);
394 
395  // put solution into output
396  CH_START(t4);
397  ierr = putPetscInChombo( a_phi, m_xx ); CHKERRQ(ierr);
398  a_phi.exchange();
399  CH_STOP(t4);
400 #endif
401  return 0;
402 }
403 
404 #ifdef CH_USE_PETSC
405 // *******************************************************
406 template <class T>
407 int PetscSolver<T>::putPetscInChombo( T& a_phi, const Vec xx )
408 {
409  PetscErrorCode ierr;
410  PetscScalar *arr;
411  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
412  DataIterator dit( dbl );
413  const int nc = a_phi.nComp();
414 
415  ierr = VecGetArray(xx,&arr); CHKERRQ(ierr);
416  for ( dit = a_phi.dataIterator() ; dit.ok() ; ++dit )
417  {
418  BaseFab<Real> &phiFab = getRegFab(a_phi,dit);
419  const Box& box = dbl.get(dit());
420  const BaseFab<Real> &gidsFab = getRegFab(m_gids,dit);
421  BoxIterator bit(box);
422  for (bit.begin(); bit.ok(); bit.next())
423  {
424  IntVect iv = bit();
425  int mm, ki = nc*(int)gidsFab(iv,0);
426  for (mm=0;mm<nc;mm++,ki++)
427  {
428  phiFab(iv,mm) = arr[ki - nc*m_gid0];
429  }
430  }
431  }
432  ierr = VecRestoreArray(xx,&arr); CHKERRQ(ierr);
433  return 0;
434 }
435 
436 // *******************************************************
437 template <class T>
438 int PetscSolver<T>::putChomboInPetsc( Vec out, const T& a_phi )
439 {
440  PetscErrorCode ierr;
441  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
442  DataIterator dit( dbl );
443  const int nc = a_phi.nComp();
444  Real idx2 = 1.e0/(this->m_dx*this->m_dx);
445 
446  // add BC stuff to RHS (EB)
447  ierr = VecSet( out, 0.); CHKERRQ(ierr);
448  for ( dit = a_phi.dataIterator() ; dit.ok() ; ++dit )
449  {
450  const BaseFab<Real> &phiFab = getRegFab(a_phi,dit);
451  const Box& box = dbl.get(dit());
452  const BaseFab<Real> &gidsFab = getRegFab(m_gids,dit);
453  BoxIterator bit(box);
454  for (bit.begin(); bit.ok(); bit.next())
455  {
456  IntVect iv = bit();
457  int mm, ki = nc*(int)gidsFab(iv,0);
458  for (mm=0;mm<nc;mm++,ki++)
459  {
460  PetscScalar v = phiFab(iv,mm);
461  v += addBCrhsValue(iv,a_phi,dit,idx2);
462  if (isnan(v)) v = 0;
463  ierr = VecSetValues(out,1,&ki,&v,INSERT_VALUES);
464  }
465  }
466  }
467  ierr = VecAssemblyBegin( out ); CHKERRQ(ierr);
468  ierr = VecAssemblyEnd( out ); CHKERRQ(ierr);
469 
470  return 0;
471 }
472 #endif
473 
474 // *******************************************************
475 // computes m_bb - A m_xx, utility method (not used)
476 //
477 template <class T>
479 {
480  Vec tempVec;
481  LevelData<FArrayBox > Residual;
482  PetscScalar alpha = -1;
483  PetscErrorCode ierr;
484 
485  ierr = VecDuplicate(m_xx, &tempVec);
486  ierr = MatMult(m_mat, m_xx, tempVec);
487  ierr = VecAXPY(tempVec,alpha,m_bb);
488 
489  IntVect idghosts = m_gids.ghostVect();
490  const DisjointBoxLayout &dbl = m_gids.disjointBoxLayout();
491  Residual.define( dbl, 0, idghosts );
492 
493  ierr = putPetscInChombo( Residual, tempVec ); CHKERRQ(ierr);
494  VecDestroy( &tempVec );
495  Residual.exchange();
496 
497  //viewBFR(&Residual );
498 
499  return 0;
500 }
501 
502 // *******************************************************
503 // PetscSolver<T>::applyOp
504 // apply the matrix - use for debugging
505 //
506 template <class T>
507 int PetscSolver<T>::applyOp( T & a_out, const T & a_phi )
508 {
509  const int nc = a_phi.nComp();
510  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
511  DataIterator dit( dbl );
512  Real idx2 = 1.e0/(this->m_dx*this->m_dx);
513  int ierr;
514 
515  // setup if needed
516  ierr = create_mat_vec( a_phi );
517 
518  // I do not know if this is setup yet
519  ierr = MatZeroEntries( m_mat ); CHKERRQ(ierr);
520  ierr = formMatrix( m_mat, a_phi ); CHKERRQ(ierr);
521 
522  // add BC stuff to RHS (EB)
523  ierr = putChomboInPetsc( m_bb, a_phi ); CHKERRQ(ierr);
524 
525  // apply op
526  ierr = MatMult( m_mat, m_bb, m_xx ); CHKERRQ(ierr);
527 
528  // get data back to Chombo
529  ierr = putPetscInChombo( a_out, m_xx ); CHKERRQ(ierr);
530 
531  return ierr;
532 }
533 
534 
535 #include "NamespaceFooter.H"
536 
537 #endif
bool ok()
Definition: BoxIterator.H:215
#define CH_TIMERS(name)
Definition: CH_Timer.H:70
virtual void exchange(const Interval &comps)
Definition: LevelDataI.H:302
Definition: LinearSolver.H:28
#define CH_SPACEDIM
Definition: SPACE.H:52
#define CH_assert(cond)
Definition: CHArray.H:37
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
void setVal(T a_x, const Box &a_bx, int a_nstart, int a_numcomp)
{ data modification functions}
Definition: BaseFabImplem.H:326
#define CH_START(tpointer)
Definition: CH_Timer.H:78
IntVect size() const
size functions
Definition: Box.H:1814
virtual void define(LinearOp< T > *a_operator, bool a_homogeneous=false)
Definition: PetscSolverI.H:78
Definition: DataIterator.H:140
iterates through the IntVects of a Box
Definition: BoxIterator.H:37
void print_memory_line(const char *)
virtual Real dx() const
Definition: LinearSolver.H:130
#define CH_TIMER(name, tpointer)
Definition: CH_Timer.H:55
Definition: EBInterface.H:45
int applyOp(T &a_phi, const T &a_rhs)
Definition: PetscSolverI.H:507
virtual ~PetscSolver()
Definition: PetscSolverI.H:71
Real m_dx
Definition: PetscSolver.H:107
Real computeResidual()
Definition: PetscSolverI.H:478
void setNull(bool n=true)
Definition: PetscSolverI.H:88
#define CH_TIME(name)
Definition: CH_Timer.H:59
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
A BoxLayout that has a concept of disjointedness.
Definition: DisjointBoxLayout.H:31
static void Error(const char *const a_msg=m_nullString, int m_exitCode=CH_DEFAULT_ERROR_CODE)
Print out message to cerr and exit with the specified exit code.
PetscSolver()
Definition: PetscSolverI.H:28
static const IntVect Zero
Definition: IntVect.H:627
Real BaseFabRealSetVal
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:465
void begin()
Definition: BoxIterator.H:146
#define CH_STOP(tpointer)
Definition: CH_Timer.H:80
int setup_solver(const T &a_phi)
Definition: PetscSolverI.H:193
Definition: PetscSolver.H:38
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
void destroy()
Definition: PetscSolverI.H:44
void next()
Definition: BoxIterator.H:164
virtual bool ok() const
return true if this iterator is still in its Layout
Definition: LayoutIterator.H:110
int create_mat_vec(const T &a_phi)
Definition: PetscSolverI.H:109
Box get(const LayoutIndex &it) const
Definition: BoxLayout.H:681