Chombo + EB: PetscSolverI.H Source File

00001 #ifdef CH_LANG_CC
00002 /*
00003  *      _______              __
00004  *     / ___/ /  ___  __ _  / /  ___
00005  *    / /__/ _ \/ _ \/  V \/ _ \/ _ \
00006  *    \___/_//_/\___/_/_/_/_.__/\___/
00007  *    Please refer to Copyright.txt, in Chombo's root directory.
00008  */
00009 #endif
00010 
00011 #ifndef _PETSCSOLVERI_H_
00012 #define _PETSCSOLVERI_H_
00013 
00014 #include "LevelData.H"
00015 #include "FluxBox.H"
00016 #include "BoxIterator.H"
00017 #include "CH_Timer.H"
00018 #include "memusage.H"
00019 #include "NamespaceHeader.H"
00020 #include "CH_OpenMP.H"
00021 
00022 // *******************************************************
00023 // PetscSolver Implementation
00024 // *******************************************************
00025 template <class T>
00026 PetscSolver<T>::PetscSolver()
00027   :m_homogeneous(false),
00028    m_mat(0), // m_xx, m_rr, m_bb;
00029    m_snes(0),
00030    m_ksp(0),
00031    m_defined(0),
00032    m_function(0),
00033    m_jacobian(0),
00034    m_null(false),
00035    m_nz_init_guess(false),
00036    m_gid0(0)
00037 {
00038   m_dx = 1.;
00039   m_prestring[0] = '\0';
00040 }
00041 // *******************************************************
00042 template <class T>
00043 void PetscSolver<T>::destroy()
00044 {
00045   if ( m_defined )
00046   {
00047 #ifdef CH_USE_PETSC
00048     if (m_mat)
00049       {
00050         MatDestroy(&m_mat);
00051         m_mat = 0;
00052       }
00053     VecDestroy(&m_bb);
00054     VecDestroy(&m_xx);
00055     VecDestroy(&m_rr);
00056 #endif
00057     m_defined = 0;
00058   }
00059 #ifdef CH_USE_PETSC
00060   if ( m_ksp )
00061     {
00062       KSPDestroy( &m_ksp );
00063       m_ksp = 0;
00064     }
00065   if ( m_snes )
00066     {
00067       SNESDestroy( &m_snes );
00068       m_snes = 0;
00069     }
00070 #endif
00071 }
00072 // *******************************************************
00073 template <class T>
00074 PetscSolver<T>::~PetscSolver()
00075 {
00076   destroy();
00077 }
00078 
00079 // *******************************************************
00080 template <class T>
00081 void PetscSolver<T>::define( Real a_dx,
00082                              bool a_homogeneous )
00083 {
00084   m_homogeneous = a_homogeneous; // not used!!!
00085   m_dx = a_dx;
00086   CH_assert(m_dx!=0.);
00087 }
00088 
00089 // *******************************************************
00090 template <class T>
00091 void PetscSolver<T>::define( LinearOp<T> *a_op,
00092                              bool a_homogeneous )
00093 {
00094   define( a_op->dx(), a_homogeneous);
00095 }
00096 
00097 // *******************************************************
00098 template <class T>
00099 void PetscSolver<T>::setNull( bool n /*= true*/ )
00100 {
00101   m_null = n; CH_assert(0);
00102 }
00103 
00104 // *******************************************************
00105 template <class T>
00106 void PetscSolver<T>::solve( T & a_phi,
00107                             const T & a_rhs )
00108 {
00109   T& phi = a_phi;
00110   const T& rhs = a_rhs;
00111   solve_private( phi, rhs );
00112 }
00113 
00114 // *******************************************************
00115 template <class T>
00116 void PetscSolver<T>::solve_mfree( T & a_phi,
00117                                   const T & a_rhs, 
00118                                   LinearOp<T> *a_op )
00119 {
00120   T& phi = a_phi;
00121   const T& rhs = a_rhs;
00122   LinearOp<T>* op = a_op;
00123   solve_mfree_private( phi, rhs , op);
00124 }
00125 
00126 // *******************************************************
00127 //   create_mat_vec
00128 //     Create 'm_mat', 'm_xx', ....  Constructs 'm_gids'.
00129 //
00130 template <class T>
00131 int PetscSolver<T>::create_mat_vec( const T& a_phi )
00132 {
00133   CH_TIME("PetscSolver::create_mat_vec");
00134 #ifdef CH_USE_PETSC
00135   const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
00136   DataIterator dit( dbl );
00137   PetscErrorCode ierr;
00138   const PetscInt nc = a_phi.nComp();
00139 #ifdef CH_MPI
00140   MPI_Comm wcomm = Chombo_MPI::comm;
00141 #else
00142   MPI_Comm wcomm = PETSC_COMM_SELF;
00143 #endif
00144 
00145   if ( !m_mat )
00146     {
00147       // print_memory_line("Before AMG set up");
00148       m_defined = 2;
00149  
00150       IntVect idghosts = a_phi.ghostVect();
00151       if (idghosts == IntVect::Zero)
00152         {
00153           MayDay::Error("PetscSolver<T>::create_mat_vec: No ghost cells in input LevelData<>.");
00154         }
00155       m_bccode.define(dbl, 1, idghosts);
00156 
00157       // global ids with ghost cells
00158       m_gids.define(dbl, 1, idghosts);
00159 
00160       // get first (zero based) id on this processor      
00161       PetscInt data = 0;
00162       for ( dit = a_phi.dataIterator() ; dit.ok() ; ++dit )
00163         {
00164           const Box &box = dbl.get(dit()); // not ghosted
00165           data += box.numPts();
00166           BaseFab<PetscInt> &gidsFab = this->m_gids[dit]; 
00167           gidsFab.setVal(-1); // flag for BC eventually
00168         }
00169       const PetscInt NN = nc*data;
00170 #ifdef CH_MPI
00171       PetscInt result;
00172       MPI_Datatype mtype;
00173       PetscDataTypeToMPIDataType(PETSC_INT,&mtype);
00174       MPI_Scan( &data, &result, 1, mtype, MPI_SUM, wcomm );
00175       m_gid0 = result - data;
00176 #else
00177       m_gid0 = 0;
00178 #endif
00179       PetscInt gid = m_gid0;
00180       for ( dit = a_phi.dataIterator() ; dit.ok() ; ++dit )
00181         {
00182           BaseFab<PetscInt> &gidsFab = this->m_gids[dit]; 
00183           const Box& box = dbl.get(dit());
00184           BoxIterator bit(box);
00185           for (bit.begin(); bit.ok(); bit.next(), gid++ )
00186             {
00187               IntVect iv = bit();
00188               gidsFab(iv,0) = gid;
00189             }
00190         }
00191       m_gids.exchange();
00192 
00193       // create matrix
00194       PetscInt nnzrow = getNNZPerRow();
00195       PetscInt *d_nnz=PETSC_NULL, *o_nnz=PETSC_NULL;
00196       
00197       if ( supportNNZExact() )
00198         { 
00199           PetscInt nglobal;
00200 #ifdef CH_MPI
00201           PetscInt nn = NN;
00202           MPI_Comm wcomm = Chombo_MPI::comm;
00203           MPI_Datatype mtype;
00204           PetscDataTypeToMPIDataType(PETSC_INT,&mtype);
00205           MPI_Allreduce(&nn,&nglobal,1,mtype,MPI_SUM,wcomm);
00206 #else
00207           nglobal = NN;
00208 #endif    
00209           ierr = PetscMalloc( (NN+1)*sizeof(PetscInt), &d_nnz ); CHKERRQ(ierr);
00210           ierr = PetscMalloc( (NN+1)*sizeof(PetscInt), &o_nnz ); CHKERRQ(ierr);
00211           for (PetscInt kk=0;kk<NN;kk++) d_nnz[kk] = o_nnz[kk] = 0;
00212 
00213           ierr = formMatrix( m_mat, &a_phi, nc*m_gid0, NN, d_nnz, o_nnz  ); CHKERRQ(ierr);
00214           CH_assert(!m_mat);
00215           // fix bounds
00216           for (PetscInt kk=0;kk<NN;kk++) 
00217             {
00218               if (d_nnz[kk] > NN) d_nnz[kk] = NN;
00219               if (o_nnz[kk] > nglobal-NN) o_nnz[kk] = nglobal-NN;
00220             }
00221           nnzrow = 0;
00222         }
00223 
00224       ierr = MatCreate(wcomm,&m_mat);CHKERRQ(ierr);
00225       ierr = MatSetSizes(m_mat,NN,NN,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
00226       ierr = MatSetBlockSize(m_mat,nc);CHKERRQ(ierr);
00227       ierr = MatSetType(m_mat,MATAIJ);CHKERRQ(ierr);
00228       ierr = MatSeqAIJSetPreallocation(m_mat,nnzrow, d_nnz);CHKERRQ(ierr);
00229       ierr = MatMPIAIJSetPreallocation(m_mat,nnzrow, d_nnz, nnzrow/2, o_nnz);CHKERRQ(ierr);
00230       ierr = MatSetOption(m_mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE) ;CHKERRQ(ierr);
00231       ierr = MatSetFromOptions( m_mat ); CHKERRQ(ierr);
00232       
00233       if ( d_nnz )
00234         {
00235           ierr = PetscFree( d_nnz );  CHKERRQ(ierr);
00236           ierr = PetscFree( o_nnz );  CHKERRQ(ierr);
00237         }
00238 
00239       // create vectors
00240       ierr = VecCreate( wcomm, &m_bb ); CHKERRQ(ierr);
00241       ierr = VecSetFromOptions( m_bb ); CHKERRQ(ierr);
00242       ierr = VecSetSizes( m_bb, NN, PETSC_DECIDE ); CHKERRQ(ierr);
00243       ierr = VecDuplicate( m_bb, &m_rr ); CHKERRQ(ierr);
00244       ierr = VecDuplicate( m_bb, &m_xx ); CHKERRQ(ierr);
00245     }
00246 #endif
00247   return 0;
00248 }
00249 
00250 // *******************************************************
00251 //    setup_solver
00252 //      - creates solver if needed.  forms matrix, sets up KSP
00253 //
00254 template <class T>
00255 int PetscSolver<T>::setup_solver( const T& a_phi )
00256 {
00257   CH_TIMERS("PetscSolver::setup_solver");
00258   CH_TIMER("solve-setup-1st", t1);
00259   CH_TIMER("solve-setup-rest", t2);
00260 #ifdef CH_USE_PETSC
00261   const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
00262   DataIterator dit( dbl );
00263   PetscErrorCode ierr;
00264   KSP ksp;
00265 #ifdef CH_MPI
00266   MPI_Comm wcomm = Chombo_MPI::comm;
00267 #else
00268   MPI_Comm wcomm = PETSC_COMM_SELF;
00269 #endif
00270 
00271   if ( m_defined == 0 )
00272     {
00273       m_defined = 2;
00274       // print_memory_line("Before AMG set up");
00275 
00276       ierr = create_mat_vec( a_phi ); CHKERRQ(ierr);
00277 
00278       CH_START(t1);
00279 
00280       // Add values to A
00281       ierr = formMatrix( m_mat, &a_phi ); CHKERRQ(ierr);
00282       {
00283         char str[256];
00284         strcpy (str,"-");
00285         strcat (str,m_prestring);
00286 #if PETSC_VERSION_GE(3,6,0)
00287         strcat (str,"pc_gamg_square_graph 20");
00288 #else
00289         strcat (str,"pc_gamg_square_graph true");
00290 #endif
00291 #if PETSC_VERSION_GE(3,7,0)
00292         ierr = PetscOptionsInsertString(PETSC_NULL,str);CHKERRQ(ierr);
00293 #else
00294         ierr = PetscOptionsInsertString(str);CHKERRQ(ierr);
00295 #endif
00296       }
00297 
00298       // create solvers
00299       PetscBool ism = PETSC_FALSE;
00300 #if PETSC_VERSION_GE(3,7,0)
00301       PetscOptionsGetBool(PETSC_NULL,m_prestring,"-ksp_monitor",&ism,PETSC_NULL);
00302 #else 
00303       PetscOptionsGetBool(m_prestring,"-ksp_monitor",&ism,PETSC_NULL);
00304 #endif
00305       if ( m_function && !m_snes )
00306         {
00307           ierr = SNESCreate( wcomm, &m_snes );                                       CHKERRQ(ierr);
00308           // these casts get rid of the 'const'. These are used in SNES as temp vector!  
00309           ierr = SNESSetFunction( m_snes, m_rr, m_function, (void*)this );         CHKERRQ(ierr);
00310           ierr = SNESSetJacobian( m_snes, m_mat, m_mat, m_jacobian, (void*)this ); CHKERRQ(ierr);
00311           //ierr = SNESSetApplicationContext( m_snes, (void*)this );                   CHKERRQ(ierr);
00312           ierr = SNESSetFromOptions( m_snes );                                       CHKERRQ(ierr);
00313           ierr = SNESGetKSP( m_snes, &ksp );                                         CHKERRQ(ierr);
00314           if (ism)
00315             {
00316               ierr = KSPMonitorSet(ksp,ksp_monitor_pout,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
00317             }
00318         }
00319       else if ( !m_function && !m_ksp )
00320         {
00321           // create the KSP so that we can set KSP parameters
00322           KSPCreate( wcomm, &m_ksp );
00323           if ( strlen(m_prestring) > 0 )
00324             {
00325               ierr = KSPSetOptionsPrefix( m_ksp, m_prestring );    CHKERRQ(ierr);
00326             }
00327           ierr = KSPSetFromOptions(m_ksp);CHKERRQ(ierr);
00328           if (ism)
00329             {
00330               ierr = KSPMonitorSet(m_ksp,ksp_monitor_pout,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00331             }
00332 #if PETSC_VERSION_GE(3,5,0)
00333           ierr = KSPSetOperators(m_ksp,m_mat,m_mat);CHKERRQ(ierr);
00334 #else
00335           ierr = KSPSetOperators(m_ksp,m_mat,m_mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
00336 #endif
00337           ierr = KSPSetInitialGuessNonzero(m_ksp, m_nz_init_guess ? PETSC_TRUE : PETSC_FALSE );CHKERRQ(ierr);
00338           ksp = m_ksp;
00339         }
00340       else CH_assert(0);
00341 
00342       { // coordinates
00343         PC pc; 
00344         PetscInt sz,ind,bs,n,m;
00345 #if PETSC_VERSION_LT(3,4,0) & PETSC_VERSION_RELEASE
00346         const PCType type;
00347 #else
00348         PCType type;
00349 #endif  
00350         ierr = KSPGetPC( ksp, &pc );     CHKERRQ(ierr);
00351         ierr = PCGetType( pc, &type );    CHKERRQ(ierr);
00352         ierr = MatGetBlockSize( m_mat, &bs );               CHKERRQ( ierr );
00353         if ( strcmp(type,PCGAMG) == 0 && bs > 1 )
00354           {
00355             PetscReal    *coords;
00356             DataIterator dit(a_phi.disjointBoxLayout());
00357         
00358             ierr = MatGetLocalSize( m_mat, &m, &n );  CHKERRQ(ierr);
00359             sz = CH_SPACEDIM*(m/bs);
00360             ierr = PetscMalloc( (sz+1)*sizeof(PetscReal), &coords ); CHKERRQ(ierr);
00361             for ( dit = a_phi.dataIterator(), ind = 0 ; dit.ok() ; ++dit )
00362               {
00363                 const Box &box = a_phi.getBoxes()[dit];
00364                 BoxIterator bit(box);
00365                 for (bit.begin(); bit.ok(); bit.next())
00366                   {
00367                     IntVect iv = bit(); // coordinate in any scaled, shifted, rotated frame.
00368                     for (PetscInt k=0; k<CH_SPACEDIM; k++) coords[ind++] = (PetscReal)iv[k];
00369                   }
00370               }
00371             CH_assert(ind==sz);
00372             ierr = PCSetCoordinates( pc, CH_SPACEDIM, sz/CH_SPACEDIM, coords ); CHKERRQ(ierr);
00373             ierr = PetscFree( coords );  CHKERRQ(ierr);
00374           }
00375       }
00376 
00377       CH_STOP(t1);
00378       // print_memory_line("After AMG set up");
00379     }
00380   else if ( m_defined == 1 )
00381     {
00382       m_defined = 2;
00383       // form A -- m_mat
00384       CH_START(t2);
00385       ierr = MatZeroEntries( m_mat );        CHKERRQ(ierr);
00386       ierr = formMatrix( m_mat, &a_phi );     CHKERRQ(ierr);
00387       if ( m_ksp )
00388         {
00389           ksp = m_ksp;
00390         }
00391       else
00392         {
00393           ierr = SNESGetKSP( m_snes, &ksp ); CHKERRQ(ierr);
00394         }
00395 #if PETSC_VERSION_GE(3,5,0)
00396       ierr = KSPSetOperators(ksp,m_mat,m_mat); CHKERRQ(ierr);
00397 #else
00398       ierr = KSPSetOperators(ksp,m_mat,m_mat,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
00399 #endif
00400       CH_STOP(t2);
00401     }
00402 #endif
00403   return 0;
00404 }
00405 
00406 // *******************************************************
00407 template <class T>
00408 int PetscSolver<T>::solve_private( T& a_phi,
00409                                    const T& a_rhs )
00410 {
00411   CH_TIMERS("PetscSolver::solve_private");
00412   CH_TIMER("formRHS", t2);
00413   CH_TIMER("solve", t3);
00414   CH_TIMER("output", t4);
00415 #ifdef CH_USE_PETSC
00416   const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
00417   PetscErrorCode ierr;
00418   const PetscInt nc = a_rhs.nComp();
00419 #ifdef CH_MPI
00420   MPI_Comm wcomm = Chombo_MPI::comm;
00421 #else
00422   MPI_Comm wcomm = PETSC_COMM_SELF;
00423 #endif
00424   // setup if needed
00425   ierr = setup_solver( a_phi );CHKERRQ(ierr);
00426 #ifdef CH_MPI
00427   MPI_Barrier(Chombo_MPI::comm);
00428 #endif
00429   CH_START(t2);
00430   // form RHS -- m_bb
00431   Real idx2 = 1.e0/(this->m_dx*this->m_dx);
00432   Real addV;
00433 
00434   //rhsOperation(a_rhs); 
00435   // this is an interface in case some operations are needed on the rhs
00436   // the default does nothing
00437   ierr = VecSetOption(m_bb,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
00438   ierr = VecSetOption(m_xx,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
00439   // add X and B from Chombo to PETSc and add stuff for EB to B
00440   ierr = VecSet( m_xx, 0.);CHKERRQ(ierr);
00441   ierr = VecSet( m_bb, 0.);CHKERRQ(ierr);
00442 
00443   DataIterator dit = a_rhs.dataIterator();
00444   int nbox=dit.size();
00445 #pragma omp parallel for private(addV,ierr)
00446   for (int mybox=0;mybox<nbox; mybox++)
00447     {
00448       const DataIndex& datInd = dit[mybox];
00449       const BaseFab<Real> &rhsFab = getRegFab(a_rhs,datInd);
00450       const BaseFab<Real> &xFab = getRegFab(a_phi,datInd);
00451       const Box& box = dbl.get(datInd);
00452       const BaseFab<PetscInt> &gidsFab = this->m_gids[datInd]; 
00453       BoxIterator bit(box);
00454       for (bit.begin(); bit.ok(); bit.next())
00455         {
00456           IntVect iv = bit();
00457           PetscInt mm, ki = nc*gidsFab(iv,0);
00458           for (mm=0;mm<nc;mm++,ki++)
00459             {
00460               PetscScalar v = rhsFab(iv,mm);
00461               addV = addBCrhsValue(iv,a_phi,datInd,idx2);
00462               if (addV == BaseFabRealSetVal)
00463                 {
00464                   v = 0;
00465                 }
00466               else
00467                 {
00468                   v += addV;
00469                 }
00470               ierr = VecSetValues(m_bb,1,&ki,&v,INSERT_VALUES);////CHKERRQ(ierr);
00471               v = xFab(iv,mm);
00472               ierr = VecSetValues(m_xx,1,&ki,&v,INSERT_VALUES);////CHKERRQ(ierr);
00473             }
00474         }
00475     }//dit
00476 
00477   ierr = VecAssemblyBegin( m_bb );CHKERRQ(ierr);
00478   ierr = VecAssemblyEnd( m_bb );CHKERRQ(ierr);
00479   //if (nonZeroInit)
00480   //{
00481   ierr = VecAssemblyBegin( m_xx );CHKERRQ(ierr);
00482   ierr = VecAssemblyEnd( m_xx );CHKERRQ(ierr);
00483   //}
00484   CH_STOP(t2);
00485   // null space for periodic BCs
00486   if ( m_null )
00487     {
00488     MatNullSpace nullsp;
00489     ierr = MatNullSpaceCreate(wcomm,PETSC_TRUE,0,PETSC_NULL,&nullsp);CHKERRQ(ierr);
00490     CH_assert(m_ksp); // not used yet, needs fix for nonlinear
00491     ierr = MatSetNullSpace( m_mat, nullsp );CHKERRQ(ierr);
00492   }
00493   // solve
00494   CH_START(t3);
00495 #ifdef CH_MPI
00496   MPI_Barrier(Chombo_MPI::comm);
00497 #endif
00498   if ( m_snes )
00499     {
00500       ierr = SNESSolve( m_snes, m_bb, m_xx );CHKERRQ(ierr);
00501 
00502     }
00503   else
00504     {
00505       ierr = KSPSolve( m_ksp, m_bb, m_xx );CHKERRQ(ierr);
00506     }
00507   CH_STOP(t3);
00508   // put solution into output
00509   CH_START(t4);
00510   ierr = putPetscInChombo( a_phi, m_xx );CHKERRQ(ierr);
00511   a_phi.exchange();
00512   CH_STOP(t4);
00513 #endif
00514   return 0;
00515 }
00516 #ifdef CH_USE_PETSC
00517 template <class T>
00518 PetscErrorCode PetscSolver<T>::apply_mfree(Mat A, Vec x, Vec f)
00519 {
00520   CH_TIME("PetscSolver::apply_mfree");
00521   //Whenever I write any PETSc code, I look forward to pulling classes back from the void.
00522   PetscFunctionBegin;
00523   void *ctx;
00524   MatShellGetContext(A, &ctx);
00525   PetscSolver<T> *tthis = (PetscSolver<T> *)ctx;
00526   tthis->putPetscInChombo( tthis->m_phi_mfree, x );
00527   tthis->m_op_mfree->applyOp(tthis->m_Lphi_mfree,tthis->m_phi_mfree,tthis->m_homogeneous);
00528   tthis->putChomboInPetsc(f, tthis->m_Lphi_mfree );
00529   PetscFunctionReturn(0);
00530 }
00531 #endif
00532 template <class T>
00533 int PetscSolver<T>::solve_mfree_private( T& a_phi,
00534                                          const T& a_rhs, 
00535                                          LinearOp<T> *a_op )
00536 {
00537   CH_TIME("PetscSolver::solve_mfree_private");
00538 #ifdef CH_USE_PETSC
00539 #ifdef CH_MPI
00540   MPI_Comm wcomm = Chombo_MPI::comm;
00541 #else
00542   MPI_Comm wcomm = PETSC_COMM_SELF;
00543 #endif
00544   PetscErrorCode ierr;
00545   // setup if needed
00546   ierr = setup_solver( a_phi );CHKERRQ(ierr);
00547 
00548   //create an operator matrix shell with same dimensions as m_mat
00549   PetscInt m, n, M, N;
00550   ierr = MatGetSize(m_mat, &M, &N);CHKERRQ(ierr); CH_assert(M == N);
00551   ierr = MatGetLocalSize(m_mat, &m, &n);CHKERRQ(ierr);
00552   Mat L; 
00553   ierr = MatCreateShell(wcomm,m,n,N,N,(void *)this,&L);CHKERRQ(ierr);
00554   ierr = MatShellSetOperation(L,MATOP_MULT,(void(*)(void))apply_mfree);
00555   m_op_mfree = a_op;
00556   //allocate space for a vector and a matrix-vector product in Chombo-land
00557   a_op->create( m_phi_mfree , a_phi);
00558   a_op->create( m_Lphi_mfree , a_rhs);
00559   
00560   // form RHS -- m_bb
00561   Real idx2 = 1.e0/(this->m_dx*this->m_dx);
00562   Real addV;
00563 
00564   // add X and B from Chombo to PETSc and add stuff for EB to B. This is just copy-pasted...
00565   ierr = VecSetOption(m_bb,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
00566   ierr = VecSetOption(m_xx,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
00567   ierr = VecSet( m_xx, 0.);CHKERRQ(ierr);
00568   ierr = VecSet( m_bb, 0.);CHKERRQ(ierr);
00569   const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout(); 
00570   const PetscInt nc = a_rhs.nComp();
00571 
00572   DataIterator dit = a_rhs.dataIterator();
00573   int nbox=dit.size();
00574   for (int mybox=0;mybox<nbox; mybox++)
00575     {
00576       const DataIndex& datInd = dit[mybox];
00577       const BaseFab<Real> &rhsFab = getRegFab(a_rhs,datInd);
00578       const BaseFab<Real> &xFab = getRegFab(a_phi,datInd);
00579       const Box& box = dbl.get(datInd);
00580       const BaseFab<PetscInt> &gidsFab = this->m_gids[datInd]; 
00581       BoxIterator bit(box);
00582       for (bit.begin(); bit.ok(); bit.next())
00583         {
00584           IntVect iv = bit();
00585           PetscInt mm, ki = nc*gidsFab(iv,0);
00586           for (mm=0;mm<nc;mm++,ki++)
00587             {
00588               PetscScalar v = rhsFab(iv,mm);
00589               addV = addBCrhsValue(iv,a_phi,datInd,idx2);
00590               if (addV == BaseFabRealSetVal)
00591                 {
00592                   v = 0;
00593                 }
00594               else
00595                 {
00596                   v += addV;
00597                 }
00598               ierr = VecSetValues(m_bb,1,&ki,&v,INSERT_VALUES); CHKERRQ(ierr);
00599               v = xFab(iv,mm);
00600               ierr = VecSetValues(m_xx,1,&ki,&v,INSERT_VALUES); CHKERRQ(ierr);
00601             }
00602         }//bit
00603     }//dit
00604 
00605   ierr = VecAssemblyBegin( m_bb );CHKERRQ(ierr);
00606   ierr = VecAssemblyEnd( m_bb );CHKERRQ(ierr);
00607   ierr = VecAssemblyBegin( m_xx );CHKERRQ(ierr);
00608   ierr = VecAssemblyEnd( m_xx );CHKERRQ(ierr);
00609   //the ksp needed here is a bit different from m_ksp, becasue the preconditioner and main operator are different 
00610   KSP ksp;
00611   KSPCreate( wcomm, &ksp );
00612   if ( strlen(m_prestring) > 0 )
00613     {
00614       ierr = KSPSetOptionsPrefix( m_ksp, m_prestring );CHKERRQ(ierr);
00615     }
00616   ierr = KSPSetFromOptions( ksp );CHKERRQ(ierr);
00617   PetscBool ism;
00618 #if PETSC_VERSION_GE(3,7,0)
00619   PetscOptionsGetBool(PETSC_NULL,m_prestring,"-ksp_monitor",&ism,PETSC_NULL);
00620 #else
00621   PetscOptionsGetBool(m_prestring,"-ksp_monitor",&ism,PETSC_NULL);
00622 #endif
00623   
00624   if (ism)
00625     {
00626       ierr = KSPMonitorSet(ksp,ksp_monitor_pout,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00627     }
00628 #if PETSC_VERSION_GE(3,5,0)
00629   ierr = KSPSetOperators(ksp,L,m_mat); CHKERRQ(ierr);
00630 #else
00631   ierr = KSPSetOperators(ksp,L,m_mat,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
00632 #endif
00633 
00634   ierr = KSPSetInitialGuessNonzero(ksp, m_nz_init_guess ? PETSC_TRUE : PETSC_FALSE ); CHKERRQ(ierr);
00635 
00636   if (0){ // coordinates
00637     PC pc; PetscInt sz,ind,bs,n,m;
00638 #if PETSC_VERSION_LT(3,4,0) & PETSC_VERSION_RELEASE
00639     const PCType type;
00640 #else
00641     PCType type;
00642 #endif  
00643     ierr = KSPGetPC( ksp, &pc );     CHKERRQ(ierr);
00644     ierr = PCGetType( pc, &type );    CHKERRQ(ierr);
00645     ierr = MatGetBlockSize( m_mat, &bs );               CHKERRQ( ierr );
00646     if ( strcmp(type,PCGAMG) == 0 && bs > 1 )
00647       {
00648         PetscReal    *coords;
00649         DataIterator dit(a_phi.disjointBoxLayout());
00650         
00651         ierr = MatGetLocalSize( m_mat, &m, &n );  CHKERRQ(ierr);
00652         sz = CH_SPACEDIM*(m/bs);
00653         ierr = PetscMalloc( (sz+1)*sizeof(PetscReal), &coords ); CHKERRQ(ierr);
00654         for ( dit = a_phi.dataIterator(), ind = 0 ; dit.ok() ; ++dit )
00655           {
00656             const Box &box = a_phi.getBoxes()[dit];
00657             BoxIterator bit(box);
00658             for (bit.begin(); bit.ok(); bit.next())
00659               {
00660                 IntVect iv = bit(); // coordinate in any scaled, shifted, rotated frame.
00661                 for (PetscInt k=0; k<CH_SPACEDIM; k++) coords[ind++] = (PetscReal)iv[k];
00662               }
00663           }
00664         CH_assert(ind==sz);
00665         ierr = PCSetCoordinates( pc, CH_SPACEDIM, sz/CH_SPACEDIM, coords ); CHKERRQ(ierr);
00666         ierr = PetscFree( coords );  CHKERRQ(ierr);
00667       }
00668   }
00669 
00670   //carry out the solve
00671   ierr = KSPSolve( ksp, m_bb, m_xx );CHKERRQ(ierr);
00672 
00673   //put solution into output
00674   ierr = putPetscInChombo( a_phi, m_xx );CHKERRQ(ierr);
00675   a_phi.exchange();
00676 
00677   //clean up 
00678   a_op->clear( m_phi_mfree);
00679   a_op->clear( m_Lphi_mfree);
00680   KSPDestroy( &ksp );
00681 #endif
00682   return 0;
00683 }
00684 
00685 #ifdef CH_USE_PETSC
00686 // *******************************************************
00687 template <class T>
00688 PetscErrorCode PetscSolver<T>::putPetscInChombo( T& a_phi, const Vec xx )
00689 {
00690   PetscErrorCode ierr;
00691   const PetscScalar *arr;
00692   const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
00693   const PetscInt nc = a_phi.nComp();
00694 
00695   ierr = VecGetArrayRead(xx,&arr);  CHKERRQ(ierr);
00696 
00697   DataIterator dit = a_phi.dataIterator();
00698   int nbox=dit.size();
00699   for (int mybox=0;mybox<nbox; mybox++)
00700     {
00701       const DataIndex& datInd = dit[mybox];
00702       BaseFab<Real> &phiFab = getRegFab(a_phi,datInd);
00703       const Box& box = dbl.get(datInd);
00704       const BaseFab<PetscInt> &gidsFab = this->m_gids[datInd]; 
00705       BoxIterator bit(box);
00706       for (bit.begin(); bit.ok(); bit.next())
00707         {
00708           IntVect iv = bit();
00709           PetscInt mm, ki = nc*gidsFab(iv,0);
00710           for (mm=0;mm<nc;mm++,ki++)
00711             {
00712               phiFab(iv,mm) = arr[ki - nc*m_gid0];
00713             }
00714         }
00715     }
00716   ierr =  VecRestoreArrayRead(xx,&arr); CHKERRQ(ierr);
00717   return 0;
00718 }
00719 
00720 // *******************************************************
00721 template <class T>
00722 PetscErrorCode PetscSolver<T>::putChomboInPetsc( Vec out, const T& a_phi )
00723 {
00724   CH_TIME("PetscSolver::putChomboInPetsc");
00725 
00726   PetscErrorCode ierr;
00727   const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
00728   const PetscInt nc = a_phi.nComp();
00729   Real idx2 = 1.e0/(this->m_dx*this->m_dx);
00730 
00731   // add BC stuff to RHS (EB)
00732   ierr = VecSet( out, 0.);  CHKERRQ(ierr);
00733 
00734   DataIterator dit = a_phi.dataIterator();
00735   int nbox=dit.size();
00736   for (int mybox=0;mybox<nbox; mybox++)
00737     {
00738       const DataIndex& datInd = dit[mybox];
00739       const BaseFab<Real> &phiFab = getRegFab(a_phi,datInd);
00740       const Box& box = dbl.get(datInd);
00741       const BaseFab<PetscInt> &gidsFab = this->m_gids[datInd]; 
00742       BoxIterator bit(box);
00743       for (bit.begin(); bit.ok(); bit.next())
00744         {
00745           IntVect iv = bit();
00746           PetscInt mm, ki = nc*gidsFab(iv,0);
00747           for (mm=0;mm<nc;mm++,ki++)
00748             {
00749               PetscScalar v = phiFab(iv,mm);
00750               v += addBCrhsValue(iv,a_phi,datInd,idx2);
00751               if (std::isnan(v)) v = 0;
00752               ierr = VecSetValues(out,1,&ki,&v,INSERT_VALUES);
00753             }
00754         }//bit
00755     }//dit
00756 
00757   ierr = VecAssemblyBegin( out );  CHKERRQ(ierr);
00758   ierr = VecAssemblyEnd( out );  CHKERRQ(ierr);
00759 
00760   return 0;
00761 }
00762 #endif
00763 
00764 // *******************************************************
00765 //  computes m_bb - A m_xx, utility method (not used)
00766 //
00767 template <class T>
00768 Real PetscSolver<T>::computeResidual()
00769 {
00770   Vec tempVec;
00771   LevelData<FArrayBox > Residual;
00772   PetscScalar alpha = -1;
00773   PetscErrorCode ierr;
00774 
00775   ierr = VecDuplicate(m_xx, &tempVec);CHKERRQ(ierr);
00776   ierr = MatMult(m_mat, m_xx, tempVec);CHKERRQ(ierr);
00777   ierr = VecAXPY(tempVec,alpha,m_bb);CHKERRQ(ierr);
00778 
00779   IntVect idghosts = m_gids.ghostVect();
00780   const DisjointBoxLayout &dbl = m_gids.disjointBoxLayout();
00781   Residual.define( dbl, 0, idghosts );
00782 
00783   ierr = putPetscInChombo( Residual, tempVec );CHKERRQ(ierr);
00784   VecDestroy( &tempVec );
00785   Residual.exchange();
00786 
00787   //viewBFR(&Residual );
00788 
00789   return 0;
00790 }
00791 
00792 // *******************************************************
00793 //  computes |phi|_inf
00794 //
00795 template <class T>
00796 Real PetscSolver<T>::normInfinity( const T& a_phi ) const
00797 {
00798   PetscErrorCode ierr;
00799   PetscReal      norm;
00800   
00801   ierr = putPetscInChombo( a_phi, m_rr );     //CHKERRQ(ierr);
00802   ierr = VecNorm(m_rr,NORM_INFINITY,&norm);CHKERRQ(ierr);
00803   //viewBFR(&Residual );
00804   return norm;
00805 }
00806 
00807 // *******************************************************
00808 // PetscSolver<T>::applyOp
00809 //  apply the matrix - use for debugging
00810 //
00811 template <class T>
00812 int PetscSolver<T>::applyOp( T & a_out, const T & a_phi )
00813 {
00814   //const int nc = a_phi.nComp();
00815   const DisjointBoxLayout &dbl = a_phi.disjointBoxLayout();
00816   DataIterator dit( dbl );
00817   PetscErrorCode ierr;
00818 
00819   // setup if needed
00820   ierr = create_mat_vec( a_phi );
00821 
00822   // add BC stuff to RHS (EB)
00823   ierr = putChomboInPetsc( m_bb, a_phi );  CHKERRQ(ierr);
00824 
00825   // apply op
00826   ierr = MatMult( m_mat, m_bb, m_xx );  CHKERRQ(ierr);
00827 
00828   // get data back to Chombo
00829   ierr = putPetscInChombo( a_out, m_xx );     CHKERRQ(ierr);
00830 
00831   return ierr;
00832 }
00833 #include "NamespaceFooter.H"
00834 #endif