Chombo + EB + MF  3.2
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 #include "CH_OpenMP.H"
21 
22 // *******************************************************
23 // PetscSolver Implementation
24 // *******************************************************
25 template <class T>
27  :m_homogeneous(false),
28  m_mat(0), // m_xx, m_rr, m_bb;
29  m_snes(0),
30  m_ksp(0),
31  m_defined(0),
32  m_function(0),
33  m_jacobian(0),
34  m_null(false),
35  m_nz_init_guess(false),
36  m_gid0(0)
37 {
38  m_dx = 1.;
39  m_prestring[0] = '\0';
40 }
41 // *******************************************************
42 template <class T>
44 {
45  if ( m_defined )
46  {
47 #ifdef CH_USE_PETSC
48  if (m_mat)
49  {
50  MatDestroy(&m_mat);
51  m_mat = 0;
52  }
53  VecDestroy(&m_bb);
54  VecDestroy(&m_xx);
55  VecDestroy(&m_rr);
56 #endif
57  m_defined = 0;
58  }
59 #ifdef CH_USE_PETSC
60  if ( m_ksp )
61  {
62  KSPDestroy( &m_ksp );
63  m_ksp = 0;
64  }
65  if ( m_snes )
66  {
67  SNESDestroy( &m_snes );
68  m_snes = 0;
69  }
70 #endif
71 }
72 // *******************************************************
73 template <class T>
75 {
76  destroy();
77 }
78 
79 // *******************************************************
80 template <class T>
82  bool a_homogeneous )
83 {
84  m_homogeneous = a_homogeneous; // not used!!!
85  m_dx = a_dx;
86  CH_assert(m_dx!=0.);
87 }
88 
89 // *******************************************************
90 template <class T>
92  bool a_homogeneous )
93 {
94  define( a_op->dx(), a_homogeneous);
95 }
96 
97 // *******************************************************
98 template <class T>
99 void PetscSolver<T>::setNull( bool n /*= true*/ )
100 {
101  m_null = n; CH_assert(0);
102 }
103 
104 // *******************************************************
105 template <class T>
106 void PetscSolver<T>::solve( T & a_phi,
107  const T & a_rhs )
108 {
109  T& phi = a_phi;
110  const T& rhs = a_rhs;
111  solve_private( phi, rhs );
112 }
113 
114 // *******************************************************
115 template <class T>
117  const T & a_rhs,
118  LinearOp<T> *a_op )
119 {
120  T& phi = a_phi;
121  const T& rhs = a_rhs;
122  LinearOp<T>* op = a_op;
123  solve_mfree_private( phi, rhs , op);
124 }
125 
126 // *******************************************************
127 // create_mat_vec
128 // Create 'm_mat', 'm_xx', .... Constructs 'm_gids'.
129 //
130 template <class T>
131 int PetscSolver<T>::create_mat_vec( const T& a_phi )
132 {
133  CH_TIME("PetscSolver::create_mat_vec");
134 #ifdef CH_USE_PETSC
135  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
136  DataIterator dit( dbl );
137  PetscErrorCode ierr;
138  const PetscInt nc = a_phi.nComp();
139 #ifdef CH_MPI
140  MPI_Comm wcomm = Chombo_MPI::comm;
141 #else
142  MPI_Comm wcomm = PETSC_COMM_SELF;
143 #endif
144 
145  if ( !m_mat )
146  {
147  // print_memory_line("Before AMG set up");
148  m_defined = 2;
149 
150  IntVect idghosts = a_phi.ghostVect();
151  if (idghosts == IntVect::Zero)
152  {
153  MayDay::Error("PetscSolver<T>::create_mat_vec: No ghost cells in input LevelData<>.");
154  }
155  m_bccode.define(dbl, 1, idghosts);
156 
157  // global ids with ghost cells
158  m_gids.define(dbl, 1, idghosts);
159 
160  // get first (zero based) id on this processor
161  PetscInt data = 0;
162  for ( dit = a_phi.dataIterator() ; dit.ok() ; ++dit )
163  {
164  const Box &box = dbl.get(dit()); // not ghosted
165  data += box.numPts();
166  BaseFab<PetscInt> &gidsFab = this->m_gids[dit];
167  gidsFab.setVal(-1); // flag for BC eventually
168  }
169  const PetscInt NN = nc*data;
170 #ifdef CH_MPI
171  PetscInt result;
172  MPI_Datatype mtype;
173  PetscDataTypeToMPIDataType(PETSC_INT,&mtype);
174  MPI_Scan( &data, &result, 1, mtype, MPI_SUM, wcomm );
175  m_gid0 = result - data;
176 #else
177  m_gid0 = 0;
178 #endif
179  PetscInt gid = m_gid0;
180  for ( dit = a_phi.dataIterator() ; dit.ok() ; ++dit )
181  {
182  BaseFab<PetscInt> &gidsFab = this->m_gids[dit];
183  const Box& box = dbl.get(dit());
184  BoxIterator bit(box);
185  for (bit.begin(); bit.ok(); bit.next(), gid++ )
186  {
187  IntVect iv = bit();
188  gidsFab(iv,0) = gid;
189  }
190  }
191  m_gids.exchange();
192 
193  // create matrix
194  PetscInt nnzrow = getNNZPerRow();
195  PetscInt *d_nnz=PETSC_NULL, *o_nnz=PETSC_NULL;
196 
197  if ( supportNNZExact() )
198  {
199  PetscInt nglobal;
200 #ifdef CH_MPI
201  PetscInt nn = NN;
202  MPI_Comm wcomm = Chombo_MPI::comm;
203  MPI_Datatype mtype;
204  PetscDataTypeToMPIDataType(PETSC_INT,&mtype);
205  MPI_Allreduce(&nn,&nglobal,1,mtype,MPI_SUM,wcomm);
206 #else
207  nglobal = NN;
208 #endif
209  ierr = PetscMalloc( (NN+1)*sizeof(PetscInt), &d_nnz ); CHKERRQ(ierr);
210  ierr = PetscMalloc( (NN+1)*sizeof(PetscInt), &o_nnz ); CHKERRQ(ierr);
211  for (PetscInt kk=0;kk<NN;kk++) d_nnz[kk] = o_nnz[kk] = 0;
212 
213  ierr = formMatrix( m_mat, &a_phi, nc*m_gid0, NN, d_nnz, o_nnz ); CHKERRQ(ierr);
214  CH_assert(!m_mat);
215  // fix bounds
216  for (PetscInt kk=0;kk<NN;kk++)
217  {
218  if (d_nnz[kk] > NN) d_nnz[kk] = NN;
219  if (o_nnz[kk] > nglobal-NN) o_nnz[kk] = nglobal-NN;
220  }
221  nnzrow = 0;
222  }
223 
224  ierr = MatCreate(wcomm,&m_mat);CHKERRQ(ierr);
225  ierr = MatSetOptionsPrefix(m_mat,m_prestring);CHKERRQ(ierr);
226  //ierr = MatSetOptionsPrefix(m_mat,"");CHKERRQ(ierr);
227 
228  pout() << "debug::matsetsize NN = " << NN << endl;
229  ierr = MatSetSizes(m_mat,NN,NN,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
230  pout() << "debug::matblocksize nc = " << nc << endl;
231  ierr = MatSetBlockSize(m_mat,nc);CHKERRQ(ierr);
232  ierr = MatSetType(m_mat,MATAIJ);CHKERRQ(ierr);
233  // ierr = MatSetOption(m_mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE) ;CHKERRQ(ierr);
234  ierr = MatSetFromOptions( m_mat ); CHKERRQ(ierr);
235  pout() << "debug::matsetpreallocation nnzrow = " << nnzrow << endl;
236  ierr = MatSeqAIJSetPreallocation(m_mat,nnzrow, d_nnz);CHKERRQ(ierr);
237  ierr = MatMPIAIJSetPreallocation(m_mat,nnzrow, d_nnz, nnzrow/2, o_nnz);CHKERRQ(ierr);
238  ierr = MatSetOption(m_mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE) ;CHKERRQ(ierr);
239 
240 #if defined(PETSC_HAVE_HYPRE)
241  ierr = MatHYPRESetPreallocation(m_mat,nnzrow, d_nnz, nnzrow/2, o_nnz);CHKERRQ(ierr);
242 #endif
243 
244  if ( d_nnz )
245  {
246  ierr = PetscFree( d_nnz ); CHKERRQ(ierr);
247  ierr = PetscFree( o_nnz ); CHKERRQ(ierr);
248  }
249 
250  // create vectors
251  ierr = MatCreateVecs(m_mat,&m_bb,&m_xx);CHKERRQ(ierr);
252  //ierr = VecCreate( wcomm, &m_bb ); CHKERRQ(ierr);
253  //ierr = VecSetFromOptions( m_bb ); CHKERRQ(ierr);
254  //ierr = VecSetSizes( m_bb, NN, PETSC_DECIDE ); CHKERRQ(ierr);
255  ierr = VecDuplicate( m_bb, &m_rr ); CHKERRQ(ierr);
256  //ierr = VecDuplicate( m_bb, &m_xx ); CHKERRQ(ierr);
257  }
258 #endif
259  return 0;
260 }
261 
262 // *******************************************************
263 // setup_solver
264 // - creates solver if needed. forms matrix, sets up KSP
265 //
266 template <class T>
267 int PetscSolver<T>::setup_solver( const T& a_phi )
268 {
269  CH_TIMERS("PetscSolver::setup_solver");
270  CH_TIMER("solve-setup-1st", t1);
271  CH_TIMER("solve-setup-rest", t2);
272 #ifdef CH_USE_PETSC
273  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
274  DataIterator dit( dbl );
275  PetscErrorCode ierr;
276  KSP ksp;
277 #ifdef CH_MPI
278  MPI_Comm wcomm = Chombo_MPI::comm;
279 #else
280  MPI_Comm wcomm = PETSC_COMM_SELF;
281 #endif
282 
283  if ( m_defined == 0 )
284  {
285  m_defined = 2;
286  // print_memory_line("Before AMG set up");
287 
288  ierr = create_mat_vec( a_phi ); CHKERRQ(ierr);
289 
290  CH_START(t1);
291 
292  // Add values to A
293  ierr = formMatrix( m_mat, &a_phi ); CHKERRQ(ierr);
294  {
295  char str[256];
296  strcpy (str,"-");
297  strcat (str,m_prestring);
298 #if PETSC_VERSION_GE(3,6,0)
299  strcat (str,"pc_gamg_square_graph 20");
300 #else
301  strcat (str,"pc_gamg_square_graph true");
302 #endif
303 #if PETSC_VERSION_GE(3,7,0)
304  ierr = PetscOptionsInsertString(PETSC_NULL,str);CHKERRQ(ierr);
305 #else
306  ierr = PetscOptionsInsertString(str);CHKERRQ(ierr);
307 #endif
308  }
309 
310  // create solvers
311  PetscBool ism = PETSC_FALSE;
312 #if PETSC_VERSION_GE(3,7,0)
313  PetscOptionsGetBool(PETSC_NULL,m_prestring,"-ksp_monitor",&ism,PETSC_NULL);
314 #else
315  PetscOptionsGetBool(m_prestring,"-ksp_monitor",&ism,PETSC_NULL);
316 #endif
317  if ( m_function && !m_snes )
318  {
319  ierr = SNESCreate( wcomm, &m_snes ); CHKERRQ(ierr);
320  // these casts get rid of the 'const'. These are used in SNES as temp vector!
321  ierr = SNESSetFunction( m_snes, m_rr, m_function, (void*)this ); CHKERRQ(ierr);
322  ierr = SNESSetJacobian( m_snes, m_mat, m_mat, m_jacobian, (void*)this ); CHKERRQ(ierr);
323  //ierr = SNESSetApplicationContext( m_snes, (void*)this ); CHKERRQ(ierr);
324  ierr = SNESSetFromOptions( m_snes ); CHKERRQ(ierr);
325  ierr = SNESGetKSP( m_snes, &ksp ); CHKERRQ(ierr);
326  if (ism)
327  {
328  ierr = KSPMonitorSet(ksp,ksp_monitor_pout,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
329  }
330  }
331  else if ( !m_function && !m_ksp )
332  {
333  // create the KSP so that we can set KSP parameters
334  KSPCreate( wcomm, &m_ksp );
335  if ( strlen(m_prestring) > 0 )
336  {
337  ierr = KSPSetOptionsPrefix( m_ksp, m_prestring ); CHKERRQ(ierr);
338  }
339  ierr = KSPSetFromOptions(m_ksp);CHKERRQ(ierr);
340  if (ism)
341  {
342  ierr = KSPMonitorSet(m_ksp,ksp_monitor_pout,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
343  }
344 #if PETSC_VERSION_GE(3,5,0)
345  ierr = KSPSetOperators(m_ksp,m_mat,m_mat);CHKERRQ(ierr);
346 #else
347  ierr = KSPSetOperators(m_ksp,m_mat,m_mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
348 #endif
349  ierr = KSPSetInitialGuessNonzero(m_ksp, m_nz_init_guess ? PETSC_TRUE : PETSC_FALSE );CHKERRQ(ierr);
350  ksp = m_ksp;
351  }
352  else CH_assert(0);
353 
354  { // coordinates
355  PC pc;
356  PetscInt sz,ind,bs,n,m;
357 #if PETSC_VERSION_LT(3,4,0) & PETSC_VERSION_RELEASE
358  const PCType type;
359 #else
360  PCType type;
361 #endif
362  ierr = KSPGetPC( ksp, &pc ); CHKERRQ(ierr);
363  ierr = PCGetType( pc, &type ); CHKERRQ(ierr);
364  ierr = MatGetBlockSize( m_mat, &bs ); CHKERRQ( ierr );
365  if ( strcmp(type,PCGAMG) == 0 && bs > 1 )
366  {
367  PetscReal *coords;
368  DataIterator dit(a_phi.disjointBoxLayout());
369 
370  ierr = MatGetLocalSize( m_mat, &m, &n ); CHKERRQ(ierr);
371  sz = CH_SPACEDIM*(m/bs);
372  ierr = PetscMalloc( (sz+1)*sizeof(PetscReal), &coords ); CHKERRQ(ierr);
373  for ( dit = a_phi.dataIterator(), ind = 0 ; dit.ok() ; ++dit )
374  {
375  const Box &box = a_phi.getBoxes()[dit];
376  BoxIterator bit(box);
377  for (bit.begin(); bit.ok(); bit.next())
378  {
379  IntVect iv = bit(); // coordinate in any scaled, shifted, rotated frame.
380  for (PetscInt k=0; k<CH_SPACEDIM; k++) coords[ind++] = (PetscReal)iv[k];
381  }
382  }
383  CH_assert(ind==sz);
384  ierr = PCSetCoordinates( pc, CH_SPACEDIM, sz/CH_SPACEDIM, coords ); CHKERRQ(ierr);
385  ierr = PetscFree( coords ); CHKERRQ(ierr);
386  }
387  }
388 
389  CH_STOP(t1);
390  // print_memory_line("After AMG set up");
391  }
392  else if ( m_defined == 1 )
393  {
394  m_defined = 2;
395  // form A -- m_mat
396  CH_START(t2);
397  ierr = MatZeroEntries( m_mat ); CHKERRQ(ierr);
398  ierr = formMatrix( m_mat, &a_phi ); CHKERRQ(ierr);
399  if ( m_ksp )
400  {
401  ksp = m_ksp;
402  }
403  else
404  {
405  ierr = SNESGetKSP( m_snes, &ksp ); CHKERRQ(ierr);
406  }
407 #if PETSC_VERSION_GE(3,5,0)
408  ierr = KSPSetOperators(ksp,m_mat,m_mat); CHKERRQ(ierr);
409 #else
410  ierr = KSPSetOperators(ksp,m_mat,m_mat,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
411 #endif
412  CH_STOP(t2);
413  }
414 #endif
415  return 0;
416 }
417 
418 // *******************************************************
419 template <class T>
421  const T& a_rhs )
422 {
423  CH_TIMERS("PetscSolver::solve_private");
424  CH_TIMER("formRHS", t2);
425  CH_TIMER("solve", t3);
426  CH_TIMER("output", t4);
427 #ifdef CH_USE_PETSC
428  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
429  PetscErrorCode ierr;
430  const PetscInt nc = a_rhs.nComp();
431 #ifdef CH_MPI
432  MPI_Comm wcomm = Chombo_MPI::comm;
433 #else
434  MPI_Comm wcomm = PETSC_COMM_SELF;
435 #endif
436  // setup if needed
437  ierr = setup_solver( a_phi );CHKERRQ(ierr);
438 #ifdef CH_MPI
439  MPI_Barrier(Chombo_MPI::comm);
440 #endif
441  CH_START(t2);
442  // form RHS -- m_bb
443  Real idx2 = 1.e0/(this->m_dx*this->m_dx);
444  Real addV;
445 
446  //rhsOperation(a_rhs);
447  // this is an interface in case some operations are needed on the rhs
448  // the default does nothing
449  ierr = VecSetOption(m_bb,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
450  ierr = VecSetOption(m_xx,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
451  // add X and B from Chombo to PETSc and add stuff for EB to B
452  ierr = VecSet( m_xx, 0.);CHKERRQ(ierr);
453  ierr = VecSet( m_bb, 0.);CHKERRQ(ierr);
454 
455  DataIterator dit = a_rhs.dataIterator();
456  int nbox=dit.size();
457 #pragma omp parallel for private(addV,ierr)
458  for (int mybox=0;mybox<nbox; mybox++)
459  {
460  const DataIndex& datInd = dit[mybox];
461  const BaseFab<Real> &rhsFab = getRegFab(a_rhs,datInd);
462  const BaseFab<Real> &xFab = getRegFab(a_phi,datInd);
463  const Box& box = dbl.get(datInd);
464  const BaseFab<PetscInt> &gidsFab = this->m_gids[datInd];
465  BoxIterator bit(box);
466  for (bit.begin(); bit.ok(); bit.next())
467  {
468  IntVect iv = bit();
469  PetscInt mm, ki = nc*gidsFab(iv,0);
470  for (mm=0;mm<nc;mm++,ki++)
471  {
472  PetscScalar v = rhsFab(iv,mm);
473  addV = addBCrhsValue(iv,a_phi,datInd,idx2);
474  if (addV == BaseFabRealSetVal)
475  {
476  v = 0;
477  }
478  else
479  {
480  v += addV;
481  }
482  ierr = VecSetValues(m_bb,1,&ki,&v,INSERT_VALUES);////CHKERRQ(ierr);
483  v = xFab(iv,mm);
484  ierr = VecSetValues(m_xx,1,&ki,&v,INSERT_VALUES);////CHKERRQ(ierr);
485  }
486  }
487  }//dit
488 
489  ierr = VecAssemblyBegin( m_bb );CHKERRQ(ierr);
490  ierr = VecAssemblyEnd( m_bb );CHKERRQ(ierr);
491  //if (nonZeroInit)
492  //{
493  ierr = VecAssemblyBegin( m_xx );CHKERRQ(ierr);
494  ierr = VecAssemblyEnd( m_xx );CHKERRQ(ierr);
495  //}
496  CH_STOP(t2);
497  // null space for periodic BCs
498  if ( m_null )
499  {
500  MatNullSpace nullsp;
501  ierr = MatNullSpaceCreate(wcomm,PETSC_TRUE,0,PETSC_NULL,&nullsp);CHKERRQ(ierr);
502  CH_assert(m_ksp); // not used yet, needs fix for nonlinear
503  ierr = MatSetNullSpace( m_mat, nullsp );CHKERRQ(ierr);
504  }
505  // solve
506  CH_START(t3);
507 #ifdef CH_MPI
508  MPI_Barrier(Chombo_MPI::comm);
509 #endif
510 // VecView(m_bb,PETSC_VIEWER_STDOUT_WORLD);
511  if ( m_snes )
512  {
513  ierr = SNESSolve( m_snes, m_bb, m_xx );CHKERRQ(ierr);
514 
515  }
516  else
517  {
518  ierr = KSPSolve( m_ksp, m_bb, m_xx );CHKERRQ(ierr);
519  }
520  CH_STOP(t3);
521  // put solution into output
522  CH_START(t4);
523 // VecView(m_xx,PETSC_VIEWER_STDOUT_WORLD);
524  ierr = putPetscInChombo( a_phi, m_xx );CHKERRQ(ierr);
525  a_phi.exchange();
526  CH_STOP(t4);
527 #endif
528  return 0;
529 }
530 #ifdef CH_USE_PETSC
531 template <class T>
532 PetscErrorCode PetscSolver<T>::apply_mfree(Mat A, Vec x, Vec f)
533 {
534  CH_TIME("PetscSolver::apply_mfree");
535  //Whenever I write any PETSc code, I look forward to pulling classes back from the void.
536  PetscFunctionBegin;
537  void *ctx;
538  MatShellGetContext(A, &ctx);
539  PetscSolver<T> *tthis = (PetscSolver<T> *)ctx;
540  tthis->putPetscInChombo( tthis->m_phi_mfree, x );
541  tthis->m_op_mfree->applyOp(tthis->m_Lphi_mfree,tthis->m_phi_mfree,tthis->m_homogeneous);
542  tthis->putChomboInPetsc(f, tthis->m_Lphi_mfree );
543  PetscFunctionReturn(0);
544 }
545 #endif
546 template <class T>
548  const T& a_rhs,
549  LinearOp<T> *a_op )
550 {
551  CH_TIME("PetscSolver::solve_mfree_private");
552 #ifdef CH_USE_PETSC
553 #ifdef CH_MPI
554  MPI_Comm wcomm = Chombo_MPI::comm;
555 #else
556  MPI_Comm wcomm = PETSC_COMM_SELF;
557 #endif
558  PetscErrorCode ierr;
559  // setup if needed
560  ierr = setup_solver( a_phi );CHKERRQ(ierr);
561 
562  //create an operator matrix shell with same dimensions as m_mat
563  PetscInt m, n, M, N;
564  ierr = MatGetSize(m_mat, &M, &N);CHKERRQ(ierr); CH_assert(M == N);
565  ierr = MatGetLocalSize(m_mat, &m, &n);CHKERRQ(ierr);
566  Mat L;
567  ierr = MatCreateShell(wcomm,m,n,N,N,(void *)this,&L);CHKERRQ(ierr);
568  ierr = MatShellSetOperation(L,MATOP_MULT,(void(*)(void))apply_mfree);
569  m_op_mfree = a_op;
570  //allocate space for a vector and a matrix-vector product in Chombo-land
571  a_op->create( m_phi_mfree , a_phi);
572  a_op->create( m_Lphi_mfree , a_rhs);
573 
574  // form RHS -- m_bb
575  Real idx2 = 1.e0/(this->m_dx*this->m_dx);
576  Real addV;
577 
578  // add X and B from Chombo to PETSc and add stuff for EB to B. This is just copy-pasted...
579  ierr = VecSetOption(m_bb,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
580  ierr = VecSetOption(m_xx,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
581  ierr = VecSet( m_xx, 0.);CHKERRQ(ierr);
582  ierr = VecSet( m_bb, 0.);CHKERRQ(ierr);
583  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
584  const PetscInt nc = a_rhs.nComp();
585 
586  DataIterator dit = a_rhs.dataIterator();
587  int nbox=dit.size();
588  for (int mybox=0;mybox<nbox; mybox++)
589  {
590  const DataIndex& datInd = dit[mybox];
591  const BaseFab<Real> &rhsFab = getRegFab(a_rhs,datInd);
592  const BaseFab<Real> &xFab = getRegFab(a_phi,datInd);
593  const Box& box = dbl.get(datInd);
594  const BaseFab<PetscInt> &gidsFab = this->m_gids[datInd];
595  BoxIterator bit(box);
596  for (bit.begin(); bit.ok(); bit.next())
597  {
598  IntVect iv = bit();
599  PetscInt mm, ki = nc*gidsFab(iv,0);
600  for (mm=0;mm<nc;mm++,ki++)
601  {
602  PetscScalar v = rhsFab(iv,mm);
603  addV = addBCrhsValue(iv,a_phi,datInd,idx2);
604  if (addV == BaseFabRealSetVal)
605  {
606  v = 0;
607  }
608  else
609  {
610  v += addV;
611  }
612  ierr = VecSetValues(m_bb,1,&ki,&v,INSERT_VALUES); CHKERRQ(ierr);
613  v = xFab(iv,mm);
614  ierr = VecSetValues(m_xx,1,&ki,&v,INSERT_VALUES); CHKERRQ(ierr);
615  }
616  }//bit
617  }//dit
618 
619  ierr = VecAssemblyBegin( m_bb );CHKERRQ(ierr);
620  ierr = VecAssemblyEnd( m_bb );CHKERRQ(ierr);
621  ierr = VecAssemblyBegin( m_xx );CHKERRQ(ierr);
622  ierr = VecAssemblyEnd( m_xx );CHKERRQ(ierr);
623  //the ksp needed here is a bit different from m_ksp, becasue the preconditioner and main operator are different
624  KSP ksp;
625  KSPCreate( wcomm, &ksp );
626  if ( strlen(m_prestring) > 0 )
627  {
628  ierr = KSPSetOptionsPrefix( m_ksp, m_prestring );CHKERRQ(ierr);
629  }
630  ierr = KSPSetFromOptions( ksp );CHKERRQ(ierr);
631  PetscBool ism;
632 #if PETSC_VERSION_GE(3,7,0)
633  PetscOptionsGetBool(PETSC_NULL,m_prestring,"-ksp_monitor",&ism,PETSC_NULL);
634 #else
635  PetscOptionsGetBool(m_prestring,"-ksp_monitor",&ism,PETSC_NULL);
636 #endif
637 
638  if (ism)
639  {
640  ierr = KSPMonitorSet(ksp,ksp_monitor_pout,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
641  }
642 #if PETSC_VERSION_GE(3,5,0)
643  ierr = KSPSetOperators(ksp,L,m_mat); CHKERRQ(ierr);
644 #else
645  ierr = KSPSetOperators(ksp,L,m_mat,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
646 #endif
647 
648  ierr = KSPSetInitialGuessNonzero(ksp, m_nz_init_guess ? PETSC_TRUE : PETSC_FALSE ); CHKERRQ(ierr);
649 
650  if (0){ // coordinates
651  PC pc; PetscInt sz,ind,bs,n,m;
652 #if PETSC_VERSION_LT(3,4,0) & PETSC_VERSION_RELEASE
653  const PCType type;
654 #else
655  PCType type;
656 #endif
657  ierr = KSPGetPC( ksp, &pc ); CHKERRQ(ierr);
658  ierr = PCGetType( pc, &type ); CHKERRQ(ierr);
659  ierr = MatGetBlockSize( m_mat, &bs ); CHKERRQ( ierr );
660  if ( strcmp(type,PCGAMG) == 0 && bs > 1 )
661  {
662  PetscReal *coords;
663  DataIterator dit(a_phi.disjointBoxLayout());
664 
665  ierr = MatGetLocalSize( m_mat, &m, &n ); CHKERRQ(ierr);
666  sz = CH_SPACEDIM*(m/bs);
667  ierr = PetscMalloc( (sz+1)*sizeof(PetscReal), &coords ); CHKERRQ(ierr);
668  for ( dit = a_phi.dataIterator(), ind = 0 ; dit.ok() ; ++dit )
669  {
670  const Box &box = a_phi.getBoxes()[dit];
671  BoxIterator bit(box);
672  for (bit.begin(); bit.ok(); bit.next())
673  {
674  IntVect iv = bit(); // coordinate in any scaled, shifted, rotated frame.
675  for (PetscInt k=0; k<CH_SPACEDIM; k++) coords[ind++] = (PetscReal)iv[k];
676  }
677  }
678  CH_assert(ind==sz);
679  ierr = PCSetCoordinates( pc, CH_SPACEDIM, sz/CH_SPACEDIM, coords ); CHKERRQ(ierr);
680  ierr = PetscFree( coords ); CHKERRQ(ierr);
681  }
682  }
683 
684  //carry out the solve
685  ierr = KSPSolve( ksp, m_bb, m_xx );CHKERRQ(ierr);
686 
687  //put solution into output
688  ierr = putPetscInChombo( a_phi, m_xx );CHKERRQ(ierr);
689  a_phi.exchange();
690 
691  //clean up
692  a_op->clear( m_phi_mfree);
693  a_op->clear( m_Lphi_mfree);
694  KSPDestroy( &ksp );
695 #endif
696  return 0;
697 }
698 
699 #ifdef CH_USE_PETSC
700 // *******************************************************
701 template <class T>
702 PetscErrorCode PetscSolver<T>::putPetscInChombo( T& a_phi, const Vec xx )
703 {
704  PetscErrorCode ierr;
705  const PetscScalar *arr;
706  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
707  const PetscInt nc = a_phi.nComp();
708 
709  ierr = VecGetArrayRead(xx,&arr); CHKERRQ(ierr);
710 
711  DataIterator dit = a_phi.dataIterator();
712  int nbox=dit.size();
713  for (int mybox=0;mybox<nbox; mybox++)
714  {
715  const DataIndex& datInd = dit[mybox];
716  BaseFab<Real> &phiFab = getRegFab(a_phi,datInd);
717  const Box& box = dbl.get(datInd);
718  const BaseFab<PetscInt> &gidsFab = this->m_gids[datInd];
719  BoxIterator bit(box);
720  for (bit.begin(); bit.ok(); bit.next())
721  {
722  IntVect iv = bit();
723  PetscInt mm, ki = nc*gidsFab(iv,0);
724  for (mm=0;mm<nc;mm++,ki++)
725  {
726  phiFab(iv,mm) = arr[ki - nc*m_gid0];
727  }
728  }
729  }
730  ierr = VecRestoreArrayRead(xx,&arr); CHKERRQ(ierr);
731  return 0;
732 }
733 
734 // *******************************************************
735 template <class T>
736 PetscErrorCode PetscSolver<T>::putChomboInPetsc( Vec out, const T& a_phi )
737 {
738  CH_TIME("PetscSolver::putChomboInPetsc");
739 
740  PetscErrorCode ierr;
741  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
742  const PetscInt nc = a_phi.nComp();
743  Real idx2 = 1.e0/(this->m_dx*this->m_dx);
744 
745  // add BC stuff to RHS (EB)
746  ierr = VecSet( out, 0.); CHKERRQ(ierr);
747 
748  DataIterator dit = a_phi.dataIterator();
749  int nbox=dit.size();
750  for (int mybox=0;mybox<nbox; mybox++)
751  {
752  const DataIndex& datInd = dit[mybox];
753  const BaseFab<Real> &phiFab = getRegFab(a_phi,datInd);
754  const Box& box = dbl.get(datInd);
755  const BaseFab<PetscInt> &gidsFab = this->m_gids[datInd];
756  BoxIterator bit(box);
757  for (bit.begin(); bit.ok(); bit.next())
758  {
759  IntVect iv = bit();
760  PetscInt mm, ki = nc*gidsFab(iv,0);
761  for (mm=0;mm<nc;mm++,ki++)
762  {
763  PetscScalar v = phiFab(iv,mm);
764  v += addBCrhsValue(iv,a_phi,datInd,idx2);
765  if (std::isnan(v)) v = 0;
766  ierr = VecSetValues(out,1,&ki,&v,INSERT_VALUES);
767  }
768  }//bit
769  }//dit
770 
771  ierr = VecAssemblyBegin( out ); CHKERRQ(ierr);
772  ierr = VecAssemblyEnd( out ); CHKERRQ(ierr);
773 
774  return 0;
775 }
776 #endif
777 
778 // *******************************************************
779 // computes m_bb - A m_xx, utility method (not used)
780 //
781 template <class T>
783 {
784  Vec tempVec;
785  LevelData<FArrayBox > Residual;
786  PetscScalar alpha = -1;
787  PetscErrorCode ierr;
788 
789  ierr = VecDuplicate(m_xx, &tempVec);CHKERRQ(ierr);
790  ierr = MatMult(m_mat, m_xx, tempVec);CHKERRQ(ierr);
791  ierr = VecAXPY(tempVec,alpha,m_bb);CHKERRQ(ierr);
792 
793  IntVect idghosts = m_gids.ghostVect();
795  Residual.define( dbl, 0, idghosts );
796 
797  ierr = putPetscInChombo( Residual, tempVec );CHKERRQ(ierr);
798  VecDestroy( &tempVec );
799  Residual.exchange();
800 
801  //viewBFR(&Residual );
802 
803  return 0;
804 }
805 
806 // *******************************************************
807 // computes |phi|_inf
808 //
809 template <class T>
810 Real PetscSolver<T>::normInfinity( const T& a_phi ) const
811 {
812  PetscErrorCode ierr;
813  PetscReal norm;
814 
815  ierr = putPetscInChombo( a_phi, m_rr ); //CHKERRQ(ierr);
816  ierr = VecNorm(m_rr,NORM_INFINITY,&norm);CHKERRQ(ierr);
817  //viewBFR(&Residual );
818  return norm;
819 }
820 
821 // *******************************************************
822 // PetscSolver<T>::applyOp
823 // apply the matrix - use for debugging
824 //
825 template <class T>
826 int PetscSolver<T>::applyOp( T & a_out, const T & a_phi )
827 {
828  //const int nc = a_phi.nComp();
829  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
830  DataIterator dit( dbl );
831  PetscErrorCode ierr;
832 
833  // setup if needed
834  ierr = create_mat_vec( a_phi );
835 
836  // add BC stuff to RHS (EB)
837  ierr = putChomboInPetsc( m_bb, a_phi ); CHKERRQ(ierr);
838 
839  // apply op
840  ierr = MatMult( m_mat, m_bb, m_xx ); CHKERRQ(ierr);
841 
842  // get data back to Chombo
843  ierr = putPetscInChombo( a_out, m_xx ); CHKERRQ(ierr);
844 
845  return ierr;
846 }
847 #include "NamespaceFooter.H"
848 #endif
std::ostream & pout()
Use this in place of std::cout for program output.
bool ok()
Definition: BoxIterator.H:281
#define CH_TIMERS(name)
Definition: CH_Timer.H:133
Definition: LinearSolver.H:28
virtual void define(Real a_dx, bool a_homogeneous=true)
Definition: PetscSolverI.H:81
#define CH_SPACEDIM
Definition: SPACE.H:51
#define CH_assert(cond)
Definition: CHArray.H:37
KSP m_ksp
Definition: PetscSolver.H:243
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
PetscInt m_gid0
Definition: PetscSolver.H:263
void setVal(T a_x, const Box &a_bx, int a_nstart, int a_numcomp)
{ data modification functions}
Definition: BaseFabImplem.H:399
Real norm(const BoxLayoutData< FArrayBox > &A, const Interval &interval, const int &p=2)
#define CH_START(tpointer)
Definition: CH_Timer.H:145
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 bool ok() const
return true if this iterator is still in its Layout
Definition: LayoutIterator.H:117
Definition: DataIterator.H:190
virtual void solve_mfree(T &a_phi, const T &a_rhs, LinearOp< T > *a_op)
Definition: PetscSolverI.H:116
iterates through the IntVects of a Box
Definition: BoxIterator.H:37
virtual Real addBCrhsValue(const IntVect &a_iv, const T &a_phi, const DataIndex &a_datInd, const Real &coeff=1)
Definition: PetscSolver.H:230
virtual void exchange(void)
Simplest case – do all components.
Definition: LevelDataI.H:451
#define CH_TIMER(name, tpointer)
Definition: CH_Timer.H:63
Definition: EBInterface.H:45
virtual int applyOp(T &a_phi, const T &a_rhs)
Definition: PetscSolverI.H:826
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
void const char const int const int const int const Real const Real * A
Definition: Lapack.H:83
int size() const
Definition: DataIterator.H:218
Real computeResidual()
Definition: PetscSolverI.H:782
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:257
LinearOp< T > * m_op_mfree
Definition: PetscSolver.H:110
#define CH_TIME(name)
Definition: CH_Timer.H:82
int solve_mfree_private(T &a_phi, const T &a_rhs, LinearOp< T > *a_op)
Definition: PetscSolverI.H:547
Mat m_mat
Definition: PetscSolver.H:238
Vec m_bb
Definition: PetscSolver.H:241
PetscErrorCode(* m_function)(SNES, Vec, Vec, void *)
Definition: PetscSolver.H:246
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:80
LevelData< BaseFab< PetscInt > > m_gids
Definition: PetscSolver.H:262
A BoxLayout that has a concept of disjointedness.
Definition: DisjointBoxLayout.H:30
Box get(const LayoutIndex &it) const
Definition: BoxLayout.H:717
const IntVect & ghostVect() const
Definition: LevelData.H:170
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.
virtual int getNNZPerRow() const
Definition: PetscSolver.H:166
Definition: BaseFab.H:81
Real normInfinity(const T &a_phi) const
Definition: PetscSolverI.H:810
PetscSolver()
Definition: PetscSolverI.H:26
static const IntVect Zero
Definition: IntVect.H:658
virtual void clear(T &a_lhs)
Definition: LinearSolver.H:70
T m_phi_mfree
Definition: PetscSolver.H:111
bool m_nz_init_guess
Definition: PetscSolver.H:256
char m_prestring[32]
Definition: PetscSolver.H:254
Real BaseFabRealSetVal
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:469
static PetscErrorCode apply_mfree(Mat A, Vec x, Vec f)
Definition: PetscSolverI.H:532
virtual void create(T &a_lhs, const T &a_rhs)=0
Vec m_xx
Definition: PetscSolver.H:241
Definition: DataIndex.H:114
T m_Lphi_mfree
Definition: PetscSolver.H:112
const DisjointBoxLayout & disjointBoxLayout() const
Definition: LevelData.H:209
void begin()
Definition: BoxIterator.H:150
#define CH_STOP(tpointer)
Definition: CH_Timer.H:150
int setup_solver(const T &a_phi)
Definition: PetscSolverI.H:267
PetscErrorCode putChomboInPetsc(Vec out, const T &a_phi)
Definition: PetscSolverI.H:736
Definition: PetscSolver.H:36
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
size_t numPts() const
void destroy()
Definition: PetscSolverI.H:43
PetscErrorCode(* m_jacobian)(SNES, Vec, Mat *, Mat *, MatStructure *, void *)
Definition: PetscSolver.H:250
void const char const int * M
Definition: Lapack.H:83
void next()
Definition: BoxIterator.H:168
void const char const int const int * N
Definition: Lapack.H:83
SNES m_snes
Definition: PetscSolver.H:242
int create_mat_vec(const T &a_phi)
Definition: PetscSolverI.H:131
PetscInt m_defined
Definition: PetscSolver.H:245
virtual Real dx() const
Definition: LinearSolver.H:128
Vec m_rr
Definition: PetscSolver.H:241
bool m_homogeneous
Definition: PetscSolver.H:188