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