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  ierr = MatSetSizes(m_mat,NN,NN,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
229  ierr = MatSetBlockSize(m_mat,nc);CHKERRQ(ierr);
230  ierr = MatSetType(m_mat,MATAIJ);CHKERRQ(ierr);
231  // ierr = MatSetOption(m_mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE) ;CHKERRQ(ierr);
232  ierr = MatSetFromOptions( m_mat ); CHKERRQ(ierr);
233  ierr = MatSeqAIJSetPreallocation(m_mat,nnzrow, d_nnz);CHKERRQ(ierr);
234  ierr = MatMPIAIJSetPreallocation(m_mat,nnzrow, d_nnz, nnzrow/2, o_nnz);CHKERRQ(ierr);
235  ierr = MatSetOption(m_mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE) ;CHKERRQ(ierr);
236 
237 #if defined(PETSC_HAVE_HYPRE)
238  ierr = MatHYPRESetPreallocation(m_mat,nnzrow, d_nnz, nnzrow/2, o_nnz);CHKERRQ(ierr);
239 #endif
240 
241  if ( d_nnz )
242  {
243  ierr = PetscFree( d_nnz ); CHKERRQ(ierr);
244  ierr = PetscFree( o_nnz ); CHKERRQ(ierr);
245  }
246 
247  // create vectors
248  ierr = MatCreateVecs(m_mat,&m_bb,&m_xx);CHKERRQ(ierr);
249  //ierr = VecCreate( wcomm, &m_bb ); CHKERRQ(ierr);
250  //ierr = VecSetFromOptions( m_bb ); CHKERRQ(ierr);
251  //ierr = VecSetSizes( m_bb, NN, PETSC_DECIDE ); CHKERRQ(ierr);
252  ierr = VecDuplicate( m_bb, &m_rr ); CHKERRQ(ierr);
253  //ierr = VecDuplicate( m_bb, &m_xx ); CHKERRQ(ierr);
254  }
255 #endif
256  return 0;
257 }
258 
259 // *******************************************************
260 // setup_solver
261 // - creates solver if needed. forms matrix, sets up KSP
262 //
263 template <class T>
264 int PetscSolver<T>::setup_solver( const T& a_phi )
265 {
266  CH_TIMERS("PetscSolver::setup_solver");
267  CH_TIMER("solve-setup-1st", t1);
268  CH_TIMER("solve-setup-rest", t2);
269 #ifdef CH_USE_PETSC
270  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
271  DataIterator dit( dbl );
272  PetscErrorCode ierr;
273  KSP ksp;
274 #ifdef CH_MPI
275  MPI_Comm wcomm = Chombo_MPI::comm;
276 #else
277  MPI_Comm wcomm = PETSC_COMM_SELF;
278 #endif
279 
280  if ( m_defined == 0 )
281  {
282  m_defined = 2;
283  // print_memory_line("Before AMG set up");
284 
285  ierr = create_mat_vec( a_phi ); CHKERRQ(ierr);
286 
287  CH_START(t1);
288 
289  // Add values to A
290  ierr = formMatrix( m_mat, &a_phi ); CHKERRQ(ierr);
291  {
292  char str[256];
293  strcpy (str,"-");
294  strcat (str,m_prestring);
295 #if PETSC_VERSION_GE(3,6,0)
296  strcat (str,"pc_gamg_square_graph 20");
297 #else
298  strcat (str,"pc_gamg_square_graph true");
299 #endif
300 #if PETSC_VERSION_GE(3,7,0)
301  ierr = PetscOptionsInsertString(PETSC_NULL,str);CHKERRQ(ierr);
302 #else
303  ierr = PetscOptionsInsertString(str);CHKERRQ(ierr);
304 #endif
305  }
306 
307  // create solvers
308  PetscBool ism = PETSC_FALSE;
309 #if PETSC_VERSION_GE(3,7,0)
310  PetscOptionsGetBool(PETSC_NULL,m_prestring,"-ksp_monitor",&ism,PETSC_NULL);
311 #else
312  PetscOptionsGetBool(m_prestring,"-ksp_monitor",&ism,PETSC_NULL);
313 #endif
314  if ( m_function && !m_snes )
315  {
316  ierr = SNESCreate( wcomm, &m_snes ); CHKERRQ(ierr);
317  // these casts get rid of the 'const'. These are used in SNES as temp vector!
318  ierr = SNESSetFunction( m_snes, m_rr, m_function, (void*)this ); CHKERRQ(ierr);
319  ierr = SNESSetJacobian( m_snes, m_mat, m_mat, m_jacobian, (void*)this ); CHKERRQ(ierr);
320  //ierr = SNESSetApplicationContext( m_snes, (void*)this ); CHKERRQ(ierr);
321  ierr = SNESSetFromOptions( m_snes ); CHKERRQ(ierr);
322  ierr = SNESGetKSP( m_snes, &ksp ); CHKERRQ(ierr);
323  if (ism)
324  {
325  ierr = KSPMonitorSet(ksp,ksp_monitor_pout,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
326  }
327  }
328  else if ( !m_function && !m_ksp )
329  {
330  // create the KSP so that we can set KSP parameters
331  KSPCreate( wcomm, &m_ksp );
332  if ( strlen(m_prestring) > 0 )
333  {
334  ierr = KSPSetOptionsPrefix( m_ksp, m_prestring ); CHKERRQ(ierr);
335  }
336  ierr = KSPSetFromOptions(m_ksp);CHKERRQ(ierr);
337  if (ism)
338  {
339  ierr = KSPMonitorSet(m_ksp,ksp_monitor_pout,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
340  }
341 #if PETSC_VERSION_GE(3,5,0)
342  ierr = KSPSetOperators(m_ksp,m_mat,m_mat);CHKERRQ(ierr);
343 #else
344  ierr = KSPSetOperators(m_ksp,m_mat,m_mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
345 #endif
346  ierr = KSPSetInitialGuessNonzero(m_ksp, m_nz_init_guess ? PETSC_TRUE : PETSC_FALSE );CHKERRQ(ierr);
347  ksp = m_ksp;
348  }
349  else CH_assert(0);
350 
351  { // coordinates
352  PC pc;
353  PetscInt sz,ind,bs,n,m;
354 #if PETSC_VERSION_LT(3,4,0) & PETSC_VERSION_RELEASE
355  const PCType type;
356 #else
357  PCType type;
358 #endif
359  ierr = KSPGetPC( ksp, &pc ); CHKERRQ(ierr);
360  ierr = PCGetType( pc, &type ); CHKERRQ(ierr);
361  ierr = MatGetBlockSize( m_mat, &bs ); CHKERRQ( ierr );
362  if ( strcmp(type,PCGAMG) == 0 && bs > 1 )
363  {
364  PetscReal *coords;
365  DataIterator dit(a_phi.disjointBoxLayout());
366 
367  ierr = MatGetLocalSize( m_mat, &m, &n ); CHKERRQ(ierr);
368  sz = CH_SPACEDIM*(m/bs);
369  ierr = PetscMalloc( (sz+1)*sizeof(PetscReal), &coords ); CHKERRQ(ierr);
370  for ( dit = a_phi.dataIterator(), ind = 0 ; dit.ok() ; ++dit )
371  {
372  const Box &box = a_phi.getBoxes()[dit];
373  BoxIterator bit(box);
374  for (bit.begin(); bit.ok(); bit.next())
375  {
376  IntVect iv = bit(); // coordinate in any scaled, shifted, rotated frame.
377  for (PetscInt k=0; k<CH_SPACEDIM; k++) coords[ind++] = (PetscReal)iv[k];
378  }
379  }
380  CH_assert(ind==sz);
381  ierr = PCSetCoordinates( pc, CH_SPACEDIM, sz/CH_SPACEDIM, coords ); CHKERRQ(ierr);
382  ierr = PetscFree( coords ); CHKERRQ(ierr);
383  }
384  }
385 
386  CH_STOP(t1);
387  // print_memory_line("After AMG set up");
388  }
389  else if ( m_defined == 1 )
390  {
391  m_defined = 2;
392  // form A -- m_mat
393  CH_START(t2);
394  ierr = MatZeroEntries( m_mat ); CHKERRQ(ierr);
395  ierr = formMatrix( m_mat, &a_phi ); CHKERRQ(ierr);
396  if ( m_ksp )
397  {
398  ksp = m_ksp;
399  }
400  else
401  {
402  ierr = SNESGetKSP( m_snes, &ksp ); CHKERRQ(ierr);
403  }
404 #if PETSC_VERSION_GE(3,5,0)
405  ierr = KSPSetOperators(ksp,m_mat,m_mat); CHKERRQ(ierr);
406 #else
407  ierr = KSPSetOperators(ksp,m_mat,m_mat,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
408 #endif
409  CH_STOP(t2);
410  }
411 #endif
412  return 0;
413 }
414 
415 // *******************************************************
416 template <class T>
418  const T& a_rhs )
419 {
420  CH_TIMERS("PetscSolver::solve_private");
421  CH_TIMER("formRHS", t2);
422  CH_TIMER("solve", t3);
423  CH_TIMER("output", t4);
424 #ifdef CH_USE_PETSC
425  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
426  PetscErrorCode ierr;
427  const PetscInt nc = a_rhs.nComp();
428 #ifdef CH_MPI
429  MPI_Comm wcomm = Chombo_MPI::comm;
430 #else
431  MPI_Comm wcomm = PETSC_COMM_SELF;
432 #endif
433  // setup if needed
434  ierr = setup_solver( a_phi );CHKERRQ(ierr);
435 #ifdef CH_MPI
436  MPI_Barrier(Chombo_MPI::comm);
437 #endif
438  CH_START(t2);
439  // form RHS -- m_bb
440  Real idx2 = 1.e0/(this->m_dx*this->m_dx);
441  Real addV;
442 
443  //rhsOperation(a_rhs);
444  // this is an interface in case some operations are needed on the rhs
445  // the default does nothing
446  ierr = VecSetOption(m_bb,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
447  ierr = VecSetOption(m_xx,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
448  // add X and B from Chombo to PETSc and add stuff for EB to B
449  ierr = VecSet( m_xx, 0.);CHKERRQ(ierr);
450  ierr = VecSet( m_bb, 0.);CHKERRQ(ierr);
451 
452  DataIterator dit = a_rhs.dataIterator();
453  int nbox=dit.size();
454 #pragma omp parallel for private(addV,ierr)
455  for (int mybox=0;mybox<nbox; mybox++)
456  {
457  const DataIndex& datInd = dit[mybox];
458  const BaseFab<Real> &rhsFab = getRegFab(a_rhs,datInd);
459  const BaseFab<Real> &xFab = getRegFab(a_phi,datInd);
460  const Box& box = dbl.get(datInd);
461  const BaseFab<PetscInt> &gidsFab = this->m_gids[datInd];
462  BoxIterator bit(box);
463  for (bit.begin(); bit.ok(); bit.next())
464  {
465  IntVect iv = bit();
466  PetscInt mm, ki = nc*gidsFab(iv,0);
467  for (mm=0;mm<nc;mm++,ki++)
468  {
469  PetscScalar v = rhsFab(iv,mm);
470  addV = addBCrhsValue(iv,a_phi,datInd,idx2);
471  if (addV == BaseFabRealSetVal)
472  {
473  v = 0;
474  }
475  else
476  {
477  v += addV;
478  }
479  ierr = VecSetValues(m_bb,1,&ki,&v,INSERT_VALUES);////CHKERRQ(ierr);
480  v = xFab(iv,mm);
481  ierr = VecSetValues(m_xx,1,&ki,&v,INSERT_VALUES);////CHKERRQ(ierr);
482  }
483  }
484  }//dit
485 
486  ierr = VecAssemblyBegin( m_bb );CHKERRQ(ierr);
487  ierr = VecAssemblyEnd( m_bb );CHKERRQ(ierr);
488  //if (nonZeroInit)
489  //{
490  ierr = VecAssemblyBegin( m_xx );CHKERRQ(ierr);
491  ierr = VecAssemblyEnd( m_xx );CHKERRQ(ierr);
492  //}
493  CH_STOP(t2);
494  // null space for periodic BCs
495  if ( m_null )
496  {
497  MatNullSpace nullsp;
498  ierr = MatNullSpaceCreate(wcomm,PETSC_TRUE,0,PETSC_NULL,&nullsp);CHKERRQ(ierr);
499  CH_assert(m_ksp); // not used yet, needs fix for nonlinear
500  ierr = MatSetNullSpace( m_mat, nullsp );CHKERRQ(ierr);
501  }
502  // solve
503  CH_START(t3);
504 #ifdef CH_MPI
505  MPI_Barrier(Chombo_MPI::comm);
506 #endif
507 // VecView(m_bb,PETSC_VIEWER_STDOUT_WORLD);
508  if ( m_snes )
509  {
510  ierr = SNESSolve( m_snes, m_bb, m_xx );CHKERRQ(ierr);
511 
512  }
513  else
514  {
515  ierr = KSPSolve( m_ksp, m_bb, m_xx );CHKERRQ(ierr);
516  }
517  CH_STOP(t3);
518  // put solution into output
519  CH_START(t4);
520 // VecView(m_xx,PETSC_VIEWER_STDOUT_WORLD);
521  ierr = putPetscInChombo( a_phi, m_xx );CHKERRQ(ierr);
522  a_phi.exchange();
523  CH_STOP(t4);
524 #endif
525  return 0;
526 }
527 #ifdef CH_USE_PETSC
528 template <class T>
529 PetscErrorCode PetscSolver<T>::apply_mfree(Mat A, Vec x, Vec f)
530 {
531  CH_TIME("PetscSolver::apply_mfree");
532  //Whenever I write any PETSc code, I look forward to pulling classes back from the void.
533  PetscFunctionBegin;
534  void *ctx;
535  MatShellGetContext(A, &ctx);
536  PetscSolver<T> *tthis = (PetscSolver<T> *)ctx;
537  tthis->putPetscInChombo( tthis->m_phi_mfree, x );
538  tthis->m_op_mfree->applyOp(tthis->m_Lphi_mfree,tthis->m_phi_mfree,tthis->m_homogeneous);
539  tthis->putChomboInPetsc(f, tthis->m_Lphi_mfree );
540  PetscFunctionReturn(0);
541 }
542 #endif
543 template <class T>
545  const T& a_rhs,
546  LinearOp<T> *a_op )
547 {
548  CH_TIME("PetscSolver::solve_mfree_private");
549 #ifdef CH_USE_PETSC
550 #ifdef CH_MPI
551  MPI_Comm wcomm = Chombo_MPI::comm;
552 #else
553  MPI_Comm wcomm = PETSC_COMM_SELF;
554 #endif
555  PetscErrorCode ierr;
556  // setup if needed
557  ierr = setup_solver( a_phi );CHKERRQ(ierr);
558 
559  //create an operator matrix shell with same dimensions as m_mat
560  PetscInt m, n, M, N;
561  ierr = MatGetSize(m_mat, &M, &N);CHKERRQ(ierr); CH_assert(M == N);
562  ierr = MatGetLocalSize(m_mat, &m, &n);CHKERRQ(ierr);
563  Mat L;
564  ierr = MatCreateShell(wcomm,m,n,N,N,(void *)this,&L);CHKERRQ(ierr);
565  ierr = MatShellSetOperation(L,MATOP_MULT,(void(*)(void))apply_mfree);
566  m_op_mfree = a_op;
567  //allocate space for a vector and a matrix-vector product in Chombo-land
568  a_op->create( m_phi_mfree , a_phi);
569  a_op->create( m_Lphi_mfree , a_rhs);
570 
571  // form RHS -- m_bb
572  Real idx2 = 1.e0/(this->m_dx*this->m_dx);
573  Real addV;
574 
575  // add X and B from Chombo to PETSc and add stuff for EB to B. This is just copy-pasted...
576  ierr = VecSetOption(m_bb,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
577  ierr = VecSetOption(m_xx,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
578  ierr = VecSet( m_xx, 0.);CHKERRQ(ierr);
579  ierr = VecSet( m_bb, 0.);CHKERRQ(ierr);
580  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
581  const PetscInt nc = a_rhs.nComp();
582 
583  DataIterator dit = a_rhs.dataIterator();
584  int nbox=dit.size();
585  for (int mybox=0;mybox<nbox; mybox++)
586  {
587  const DataIndex& datInd = dit[mybox];
588  const BaseFab<Real> &rhsFab = getRegFab(a_rhs,datInd);
589  const BaseFab<Real> &xFab = getRegFab(a_phi,datInd);
590  const Box& box = dbl.get(datInd);
591  const BaseFab<PetscInt> &gidsFab = this->m_gids[datInd];
592  BoxIterator bit(box);
593  for (bit.begin(); bit.ok(); bit.next())
594  {
595  IntVect iv = bit();
596  PetscInt mm, ki = nc*gidsFab(iv,0);
597  for (mm=0;mm<nc;mm++,ki++)
598  {
599  PetscScalar v = rhsFab(iv,mm);
600  addV = addBCrhsValue(iv,a_phi,datInd,idx2);
601  if (addV == BaseFabRealSetVal)
602  {
603  v = 0;
604  }
605  else
606  {
607  v += addV;
608  }
609  ierr = VecSetValues(m_bb,1,&ki,&v,INSERT_VALUES); CHKERRQ(ierr);
610  v = xFab(iv,mm);
611  ierr = VecSetValues(m_xx,1,&ki,&v,INSERT_VALUES); CHKERRQ(ierr);
612  }
613  }//bit
614  }//dit
615 
616  ierr = VecAssemblyBegin( m_bb );CHKERRQ(ierr);
617  ierr = VecAssemblyEnd( m_bb );CHKERRQ(ierr);
618  ierr = VecAssemblyBegin( m_xx );CHKERRQ(ierr);
619  ierr = VecAssemblyEnd( m_xx );CHKERRQ(ierr);
620  //the ksp needed here is a bit different from m_ksp, becasue the preconditioner and main operator are different
621  KSP ksp;
622  KSPCreate( wcomm, &ksp );
623  if ( strlen(m_prestring) > 0 )
624  {
625  ierr = KSPSetOptionsPrefix( m_ksp, m_prestring );CHKERRQ(ierr);
626  }
627  ierr = KSPSetFromOptions( ksp );CHKERRQ(ierr);
628  PetscBool ism;
629 #if PETSC_VERSION_GE(3,7,0)
630  PetscOptionsGetBool(PETSC_NULL,m_prestring,"-ksp_monitor",&ism,PETSC_NULL);
631 #else
632  PetscOptionsGetBool(m_prestring,"-ksp_monitor",&ism,PETSC_NULL);
633 #endif
634 
635  if (ism)
636  {
637  ierr = KSPMonitorSet(ksp,ksp_monitor_pout,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
638  }
639 #if PETSC_VERSION_GE(3,5,0)
640  ierr = KSPSetOperators(ksp,L,m_mat); CHKERRQ(ierr);
641 #else
642  ierr = KSPSetOperators(ksp,L,m_mat,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
643 #endif
644 
645  ierr = KSPSetInitialGuessNonzero(ksp, m_nz_init_guess ? PETSC_TRUE : PETSC_FALSE ); CHKERRQ(ierr);
646 
647  if (0){ // coordinates
648  PC pc; PetscInt sz,ind,bs,n,m;
649 #if PETSC_VERSION_LT(3,4,0) & PETSC_VERSION_RELEASE
650  const PCType type;
651 #else
652  PCType type;
653 #endif
654  ierr = KSPGetPC( ksp, &pc ); CHKERRQ(ierr);
655  ierr = PCGetType( pc, &type ); CHKERRQ(ierr);
656  ierr = MatGetBlockSize( m_mat, &bs ); CHKERRQ( ierr );
657  if ( strcmp(type,PCGAMG) == 0 && bs > 1 )
658  {
659  PetscReal *coords;
660  DataIterator dit(a_phi.disjointBoxLayout());
661 
662  ierr = MatGetLocalSize( m_mat, &m, &n ); CHKERRQ(ierr);
663  sz = CH_SPACEDIM*(m/bs);
664  ierr = PetscMalloc( (sz+1)*sizeof(PetscReal), &coords ); CHKERRQ(ierr);
665  for ( dit = a_phi.dataIterator(), ind = 0 ; dit.ok() ; ++dit )
666  {
667  const Box &box = a_phi.getBoxes()[dit];
668  BoxIterator bit(box);
669  for (bit.begin(); bit.ok(); bit.next())
670  {
671  IntVect iv = bit(); // coordinate in any scaled, shifted, rotated frame.
672  for (PetscInt k=0; k<CH_SPACEDIM; k++) coords[ind++] = (PetscReal)iv[k];
673  }
674  }
675  CH_assert(ind==sz);
676  ierr = PCSetCoordinates( pc, CH_SPACEDIM, sz/CH_SPACEDIM, coords ); CHKERRQ(ierr);
677  ierr = PetscFree( coords ); CHKERRQ(ierr);
678  }
679  }
680 
681  //carry out the solve
682  ierr = KSPSolve( ksp, m_bb, m_xx );CHKERRQ(ierr);
683 
684  //put solution into output
685  ierr = putPetscInChombo( a_phi, m_xx );CHKERRQ(ierr);
686  a_phi.exchange();
687 
688  //clean up
689  a_op->clear( m_phi_mfree);
690  a_op->clear( m_Lphi_mfree);
691  KSPDestroy( &ksp );
692 #endif
693  return 0;
694 }
695 
696 #ifdef CH_USE_PETSC
697 // *******************************************************
698 template <class T>
699 PetscErrorCode PetscSolver<T>::putPetscInChombo( T& a_phi, const Vec xx )
700 {
701  PetscErrorCode ierr;
702  const PetscScalar *arr;
703  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
704  const PetscInt nc = a_phi.nComp();
705 
706  ierr = VecGetArrayRead(xx,&arr); CHKERRQ(ierr);
707 
708  DataIterator dit = a_phi.dataIterator();
709  int nbox=dit.size();
710  for (int mybox=0;mybox<nbox; mybox++)
711  {
712  const DataIndex& datInd = dit[mybox];
713  BaseFab<Real> &phiFab = getRegFab(a_phi,datInd);
714  const Box& box = dbl.get(datInd);
715  const BaseFab<PetscInt> &gidsFab = this->m_gids[datInd];
716  BoxIterator bit(box);
717  for (bit.begin(); bit.ok(); bit.next())
718  {
719  IntVect iv = bit();
720  PetscInt mm, ki = nc*gidsFab(iv,0);
721  for (mm=0;mm<nc;mm++,ki++)
722  {
723  phiFab(iv,mm) = arr[ki - nc*m_gid0];
724  }
725  }
726  }
727  ierr = VecRestoreArrayRead(xx,&arr); CHKERRQ(ierr);
728  return 0;
729 }
730 
731 // *******************************************************
732 template <class T>
733 PetscErrorCode PetscSolver<T>::putChomboInPetsc( Vec out, const T& a_phi )
734 {
735  CH_TIME("PetscSolver::putChomboInPetsc");
736 
737  PetscErrorCode ierr;
738  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
739  const PetscInt nc = a_phi.nComp();
740  Real idx2 = 1.e0/(this->m_dx*this->m_dx);
741 
742  // add BC stuff to RHS (EB)
743  ierr = VecSet( out, 0.); CHKERRQ(ierr);
744 
745  DataIterator dit = a_phi.dataIterator();
746  int nbox=dit.size();
747  for (int mybox=0;mybox<nbox; mybox++)
748  {
749  const DataIndex& datInd = dit[mybox];
750  const BaseFab<Real> &phiFab = getRegFab(a_phi,datInd);
751  const Box& box = dbl.get(datInd);
752  const BaseFab<PetscInt> &gidsFab = this->m_gids[datInd];
753  BoxIterator bit(box);
754  for (bit.begin(); bit.ok(); bit.next())
755  {
756  IntVect iv = bit();
757  PetscInt mm, ki = nc*gidsFab(iv,0);
758  for (mm=0;mm<nc;mm++,ki++)
759  {
760  PetscScalar v = phiFab(iv,mm);
761  v += addBCrhsValue(iv,a_phi,datInd,idx2);
762  if (std::isnan(v)) v = 0;
763  ierr = VecSetValues(out,1,&ki,&v,INSERT_VALUES);
764  }
765  }//bit
766  }//dit
767 
768  ierr = VecAssemblyBegin( out ); CHKERRQ(ierr);
769  ierr = VecAssemblyEnd( out ); CHKERRQ(ierr);
770 
771  return 0;
772 }
773 #endif
774 
775 // *******************************************************
776 // computes m_bb - A m_xx, utility method (not used)
777 //
778 template <class T>
780 {
781  Vec tempVec;
782  LevelData<FArrayBox > Residual;
783  PetscScalar alpha = -1;
784  PetscErrorCode ierr;
785 
786  ierr = VecDuplicate(m_xx, &tempVec);CHKERRQ(ierr);
787  ierr = MatMult(m_mat, m_xx, tempVec);CHKERRQ(ierr);
788  ierr = VecAXPY(tempVec,alpha,m_bb);CHKERRQ(ierr);
789 
790  IntVect idghosts = m_gids.ghostVect();
792  Residual.define( dbl, 0, idghosts );
793 
794  ierr = putPetscInChombo( Residual, tempVec );CHKERRQ(ierr);
795  VecDestroy( &tempVec );
796  Residual.exchange();
797 
798  //viewBFR(&Residual );
799 
800  return 0;
801 }
802 
803 // *******************************************************
804 // computes |phi|_inf
805 //
806 template <class T>
807 Real PetscSolver<T>::normInfinity( const T& a_phi ) const
808 {
809  PetscErrorCode ierr;
810  PetscReal norm;
811 
812  ierr = putPetscInChombo( a_phi, m_rr ); //CHKERRQ(ierr);
813  ierr = VecNorm(m_rr,NORM_INFINITY,&norm);CHKERRQ(ierr);
814  //viewBFR(&Residual );
815  return norm;
816 }
817 
818 // *******************************************************
819 // PetscSolver<T>::applyOp
820 // apply the matrix - use for debugging
821 //
822 template <class T>
823 int PetscSolver<T>::applyOp( T & a_out, const T & a_phi )
824 {
825  //const int nc = a_phi.nComp();
826  const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
827  DataIterator dit( dbl );
828  PetscErrorCode ierr;
829 
830  // setup if needed
831  ierr = create_mat_vec( a_phi );
832 
833  // add BC stuff to RHS (EB)
834  ierr = putChomboInPetsc( m_bb, a_phi ); CHKERRQ(ierr);
835 
836  // apply op
837  ierr = MatMult( m_mat, m_bb, m_xx ); CHKERRQ(ierr);
838 
839  // get data back to Chombo
840  ierr = putPetscInChombo( a_out, m_xx ); CHKERRQ(ierr);
841 
842  return ierr;
843 }
844 #include "NamespaceFooter.H"
845 #endif
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:417
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:823
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:779
PetscErrorCode putPetscInChombo(T &a_phi, const Vec xx)
Definition: PetscSolverI.H:699
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:544
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:807
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:529
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:264
PetscErrorCode putChomboInPetsc(Vec out, const T &a_phi)
Definition: PetscSolverI.H:733
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