Chombo + EB: EBAMRPoissonOp.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 _EBAMRPOISSONOP_H_
00012 #define _EBAMRPOISSONOP_H_
00013 
00014 #include "REAL.H"
00015 #include "Box.H"
00016 #include "FArrayBox.H"
00017 #include "Vector.H"
00018 #include <map>
00019 #include "RefCountedPtr.H"
00020 
00021 #include "AMRMultiGrid.H"
00022 
00023 #include "EBIndexSpace.H"
00024 #include "EBCellFAB.H"
00025 #include "EBCellFactory.H"
00026 ///only use the non-aggregated for testing
00027 #define USE_NONAGG 0
00028 
00029 #if USE_NONAGG==1
00030 #define EBSTENCIL_T NonAggregatedEBStencil
00031 #else
00032 #define EBSTENCIL_T EBStencil
00033 #endif
00034 
00035 #include "EBStencil.H"
00036 #include "NonAggregatedEBStencil.H"
00037 
00038 #include "EBLevelDataOps.H"
00039 #include "BaseEBBC.H"
00040 #include "BaseDomainBC.H"
00041 #include "CFIVS.H"
00042 #include "EBFastFR.H"
00043 #include "EBMGAverage.H"
00044 #include "EBMGInterp.H"
00045 #include "PolyGeom.H"
00046 #include "EBQuadCFInterp.H"
00047 #include "EBLevelGrid.H"
00048 #include "AMRTGA.H"
00049 #include "EBAMRIO.H"
00050 #include "AMRPoissonOp.H"
00051 #include "CFRegion.H"
00052 #include "NamespaceHeader.H"
00053 
00054 #ifdef CH_USE_PETSC
00055 #include "petsc.h"
00056 #include "petscmat.h"
00057 #include "petscksp.h"
00058 #include "petscviewer.h"
00059 #endif
00060 
00061 #if CH_SPACEDIM==2
00062 #define EBAMRPO_NUMSTEN 4
00063 #elif CH_SPACEDIM==3
00064 #define EBAMRPO_NUMSTEN 8
00065 #else
00066 void THIS_IS_AN_ERROR_MESSAGE(void)
00067 {
00068   THIS_WILL_ONLY_COMPILE_WHEN_CH_SPACEDIM_IS_2_OR_3;
00069 }
00070 #endif
00071 
00072 ///
00073 /**
00074    Operator to solve (alpha + beta lapl)phi = rhs.   This follows the AMRLevelOp interface.
00075 */
00076 class EBAMRPoissonOp: public LevelTGAHelmOp<LevelData<EBCellFAB>, EBFluxFAB >
00077 {
00078 public:
00079 
00080 #ifdef CH_USE_HDF5
00081   ///
00082   virtual void dumpAMR(Vector<LevelData<EBCellFAB>*>& a_data, string name)
00083   {
00084     writeEBAMRname(&a_data, name.c_str());
00085   }
00086 
00087   virtual void dumpLevel(LevelData<EBCellFAB>& a_data, string name)
00088   {
00089     writeEBLevelname(&a_data, name.c_str());
00090   }
00091 #endif
00092 
00093 #ifdef CH_USE_PETSC
00094   ///Fill a_petsc_mat with a matrix that describes the operator stencil
00095   int getPetscMatrix(Mat& a_petsc_mat);
00096 
00097   ///put the real components into comp 0, the imaginary ones into comp 1
00098   int getLevelDataFromPetscVector(LevelData<EBCellFAB>& a_data, const Vec& a_petsc_vec_real, const Vec& a_petsc_vec_imag );
00099 
00100   ///
00101   /**
00102      get the matrix indicies associated with every point on the level
00103      Not really going to work if level does not cover the domain
00104   */
00105   static int getMatrixIndexingLD(LevelData<BaseEBCellFAB<int> >&  a_gids, int & a_data,
00106                                  const EBLevelGrid&               a_eblg,
00107                                  const IntVect    &               a_ghostCellsPhi,
00108                                  const bool       &               a_hasCoar);
00109 #endif
00110   static void
00111   getAggregatedLayout(DisjointBoxLayout           & a_dblCoar,
00112                       const ProblemDomain         & a_domainCoar,
00113                       const EBIndexSpace * const    a_ebisPtr,
00114                       const int                   & a_maxBoxSize);
00115 
00116   ///
00117   /**
00118      version that does not fill ebislcoar
00119    */
00120   static bool getCoarserLayouts(DisjointBoxLayout&       a_dblCoar,
00121                                 ProblemDomain&           a_domainCoar,
00122                                 const DisjointBoxLayout& a_dblFine,
00123                                 const EBISLayout&        a_ebislFine,
00124                                 const ProblemDomain&     a_domainFine,
00125                                 int                      a_refToCoar,
00126                                 const EBIndexSpace*      a_ebisPtr,
00127                                 int                      a_maxBoxSize,
00128                                 bool&                    a_layoutChanged,
00129                                 int                      a_testRef = 2);
00130 
00131   static Real staticMaxNorm(const LevelData<EBCellFAB>& a_rhs, const EBLevelGrid& a_eblg);
00132 
00133   //for tga to reset stuff
00134   virtual void setAlphaAndBeta(const Real& a_alpha,
00135                                const Real& a_beta);
00136 
00137   //another tgaism
00138   virtual void diagonalScale(LevelData<EBCellFAB>& a_rhs,
00139                              bool                  a_kappaWeighted = true);
00140   virtual void divideByIdentityCoef(LevelData<EBCellFAB> & a_rhs)
00141   {
00142     //no acoef here.
00143   }
00144 
00145   virtual void kappaScale(LevelData<EBCellFAB> & a_rhs)
00146   {
00147     //since this is a constant coefficient operator with a=1,
00148     //this means the same thing as diagonal scale
00149     diagonalScale(a_rhs);
00150 
00151   }
00152 
00153   ///returns m_dx, such function is required by some LinearSolvers
00154   Real dx() const
00155   {
00156     return m_dx[0];
00157   }
00158 
00159   ///a leveltgaism
00160   virtual void fillGrad(const LevelData<EBCellFAB>& a_phi)
00161   {;}
00162 
00163   ///
00164   /**
00165      dump stencil as matrix to stdout
00166    */
00167   class StencilIndex
00168   {
00169   public:
00170     StencilIndex(VolIndex& a_vof,
00171                  int&      a_phase)
00172       {
00173         m_vof   = a_vof;
00174         m_phase = a_phase;
00175       }
00176 
00177     StencilIndex()
00178     {
00179       MayDay::Error("No weak construction of StencilIndex class.");
00180     }
00181 
00182     StencilIndex& operator= (const StencilIndex& a_sin)
00183     {
00184       m_vof   =a_sin.m_vof;
00185       m_phase = a_sin.m_phase;
00186       return *this;
00187     }
00188 
00189     VolIndex vof() const
00190     {
00191       return m_vof;
00192     }
00193 
00194     int phase() const
00195     {
00196       return m_phase;
00197     }
00198 
00199     bool operator!=(const StencilIndex& a_sin) const
00200     {
00201       return m_vof!= a_sin.m_vof || m_phase!=a_sin.m_phase;
00202     }
00203 
00204   protected:
00205     VolIndex m_vof;
00206     int      m_phase;
00207   };
00208 
00209   class StencilIndexComparator
00210   {
00211   public:
00212     bool operator() (const StencilIndex& a_s1, const StencilIndex& a_s2) const
00213     {
00214       /*
00215       int p1 = a_s1.phase();
00216       int p2 = a_s2.phase();
00217       if (p1 == p2)
00218         {
00219           return (a_s1.vof()<a_s2.vof());
00220         }
00221       else
00222         {
00223           const IntVect& iv1 = a_s1.vof().gridIndex();
00224           const IntVect& iv2 = a_s2.vof().gridIndex();
00225           if (iv1 == iv2)
00226             {
00227               return (p1<p2);
00228             }
00229           else
00230             {
00231               return (iv1<iv2);
00232             }
00233         }
00234       */
00235       const IntVect& iv1 = a_s1.vof().gridIndex();
00236       const IntVect& iv2 = a_s2.vof().gridIndex();
00237       for (int idir=0; idir<SpaceDim; ++idir)
00238         {
00239           if (iv1[idir] != iv2[idir])
00240             {
00241               return (iv1[idir]<iv2[idir]);
00242             }
00243         }
00244       int p1 = a_s1.phase();
00245       int p2 = a_s2.phase();
00246       if (p1 == p2)
00247         {
00248           int ci1 = a_s1.vof().cellIndex();
00249           int ci2 = a_s2.vof().cellIndex();
00250           return (ci1<ci2);
00251         }
00252       else
00253         {
00254           return (p1<p2);
00255         }
00256     }
00257   };
00258 
00259   /*
00260     Calculate stencil weights for each VoF. A mapping of each VolIndex
00261     an integer is also performed. The first output is this mapping, as
00262     a matrix that can be read into Matlab/Octave. The second output is
00263     the stencil weight for each VoF pair, in matrix form.
00264     The volume fraction is an additional output of the mapping matrix.
00265     The matrices are output in the format:
00266 
00267     M_(RES)(imap, :) = [ gridIndex[0] gridIndex[1] (gridIndex[2]) alphaWeight kappa];
00268     where 'RES' is the resolution at the EBPoissonOp's level, 'imap' is
00269     the integer the VoF has been mapper to, 'alphaWeight' is the
00270     weighting used to get good conditioning of the matrix, and kappa is
00271     the VoF's volume fraction.
00272 
00273     L_(RES)(imap, jmap) = value;
00274     where 'RES' is as above, 'imap' is the integer that the VoF whose
00275     stencil is being computed has been mapped to, 'jmap' is the integer
00276     for a VoF in that stencil. 'value' is the stencil weight.
00277    */
00278   void dumpStencilMatrix();
00279 
00280   /*
00281     Calculates the stencil for VoFs at the domain boundary.
00282    */
00283   void getDomainFluxStencil(      VoFStencil& a_stencil,
00284                             const VolIndex&   a_vof,
00285                             const int         a_comp,
00286                             const DataIndex&  a_dit);
00287   /*
00288     A method for testing EBAMRPoissonOp::dumpStencilMatrix. Since it
00289     uses calls to 'applyOp' to get the stencil weights, it is slow
00290     but correct.
00291    */
00292   void dumpReferenceStencilMatrix();
00293 
00294   virtual void getFlux(EBFluxFAB&                    a_flux,
00295                        const LevelData<EBCellFAB>&   a_data,
00296                        const Box&                    a_grid,
00297                        const DataIndex&              a_dit,
00298                        Real                          a_scale)
00299   {
00300     for (int idir = 0; idir < SpaceDim; idir++)
00301       {
00302         Box ghostedBox = a_grid;
00303         ghostedBox.grow(1);
00304         ghostedBox.grow(idir,-1);
00305         ghostedBox &= m_eblg.getDomain();
00306 
00307         getFlux(a_flux[idir], a_data[a_dit], ghostedBox, a_grid,
00308                 m_eblg.getDomain(),
00309                 m_eblg.getEBISL()[a_dit], m_dx, idir);
00310       }
00311   }
00312   EBLevelGrid getEBLG()
00313   {
00314     return m_eblg;
00315   }
00316   EBLevelGrid getEBLGCoarMG()
00317   {
00318     return m_eblgCoarMG;
00319   }
00320   static void setOperatorTime(Real a_time)
00321   {
00322     s_time = a_time;
00323   }
00324 
00325   //conforms to AMRTGA
00326   virtual void setTime(Real a_time)
00327   {
00328     setOperatorTime(a_time);
00329   }
00330   ///
00331   virtual ~EBAMRPoissonOp();
00332 
00333   ///
00334   EBAMRPoissonOp();
00335 
00336   ///
00337   /** a_residual = a_rhs - L(a_phiFine, a_phi)   no coaser AMR level*/
00338   void AMRResidualNC(LevelData<EBCellFAB>&       a_residual,
00339                      const LevelData<EBCellFAB>& a_phiFine,
00340                      const LevelData<EBCellFAB>& a_phi,
00341                      const LevelData<EBCellFAB>& a_rhs,
00342                      bool a_homogeneousBC,
00343                      AMRLevelOp<LevelData<EBCellFAB> >* a_finerOp);
00344 
00345 
00346   ///
00347   /** apply AMR operator   no coaser AMR level*/
00348   void AMROperatorNC(LevelData<EBCellFAB>&       a_LofPhi,
00349                      const LevelData<EBCellFAB>& a_phiFine,
00350                      const LevelData<EBCellFAB>& a_phi,
00351                      bool a_homogeneousBC,
00352                      AMRLevelOp<LevelData<EBCellFAB> >* a_finerOp);
00353 
00354   ///
00355   /**
00356      If you are approaching this operator from this interface, consider backing away and using
00357      EBAMRPoissonOpFactory to generate these objects.   Really.
00358      a_eblgFine,    : grid at finer  level \\
00359      a_eblg,        : grid at this  level \\
00360      a_eblgCoar,    : grid at coarser level \\
00361      a_eblgCoarMG,  : grid at intermediate multigrid level \\
00362      a_domainBC,    : domain boundary conditions at this level   \\
00363      a_ebBC:          eb boundary conditions at this level   \\
00364      a_dx:            grid spacing at this level \\
00365      a_origin:        offset to lowest corner of the domain \\
00366      a_refToFine:     refinement ratio to finer level \\
00367      a_refToCoar:     refinement ratio to coarser level \\
00368      a_hasFiner:      true if there is a finer AMR level, false otherwise. \\
00369      a_hasCoarser:    true  if there is a coarser AMR level. \\
00370      a_hasCoarserMG:    true  if there is a coarser MultiGrid level. \\
00371      a_preCondIters:  number of iterations to do for pre-conditioning \\
00372      a_relaxType:     0 means point Jacobi, 1 is Gauss-Seidel. \\
00373      a_alpha:         coefficent of identity \\
00374      a_beta:          coefficient of laplacian.\\
00375      a_ghostCellsPhi:  Number of ghost cells in phi, correction\\
00376      a_ghostCellsRhs:  Number of ghost cells in RHS, residual, lphi\\
00377      Ghost cell arguments are there for caching reasons.  Once you set them, an error is thrown if
00378      you send in data that does not match.
00379   */
00380 EBAMRPoissonOp(const EBLevelGrid &                  a_eblgFine,
00381                const EBLevelGrid &                  a_eblg,
00382                const EBLevelGrid &                  a_eblgCoar,
00383                const EBLevelGrid &                  a_eblgCoarMG,
00384                const RefCountedPtr<EBQuadCFInterp>& a_quadCFI,
00385                const RefCountedPtr<BaseDomainBC>&   a_domainBC,
00386                const RefCountedPtr<BaseEBBC>&       a_ebBC,
00387                const RealVect&                      a_dx,
00388                const RealVect&                      a_dxCoar,
00389                const RealVect&                      a_origin,
00390                const int&                           a_refToFine,
00391                const int&                           a_refToCoar,
00392                const bool&                          a_hasFine,
00393                const bool&                          a_hasCoar,
00394                const bool&                          a_hasMGObjects,
00395                const bool&                          a_layoutChanged,
00396                const int&                           a_numPreCondIters,
00397                const int&                           a_relaxType,
00398                const Real&                          a_alpha,
00399                const Real&                          a_beta,
00400                const IntVect&                       a_ghostCellsPhi,
00401                const IntVect&                       a_ghostCellsRHS,
00402                int                                  a_testRef = 2);
00403 
00404   //MGOp operations.  no finer or coarser
00405 
00406   ///
00407   /**
00408    */
00409   virtual void residual(LevelData<EBCellFAB>&       a_residual,
00410                         const LevelData<EBCellFAB>& a_phi,
00411                         const LevelData<EBCellFAB>& a_rhs,
00412                         bool                        a_homogeneousPhysBC=false);
00413 
00414   ///
00415   /**
00416      For debugging purpose, get the matrix of the op (could be change from time to time due to EBBC update) by calculating 0 - L(phi(i)) for i=0,1,...,n, each time phi(i)=1 for i-th cell and 0 for all other cells
00417    */
00418 
00419   virtual void getOpMatrix(const LevelData<EBCellFAB>& a_phi,
00420                            const LevelData<EBCellFAB>& a_rhs);
00421 
00422   ///
00423   /**
00424    */
00425   virtual void preCond(LevelData<EBCellFAB>&       a_opPhi,
00426                        const LevelData<EBCellFAB>& a_phi);
00427 
00428   ///
00429   /**
00430      This function assumes that coarse-fine boundary condtions have
00431      been dealt with.
00432   */
00433   virtual void applyOp(LevelData<EBCellFAB>&             a_opPhi,
00434                        const LevelData<EBCellFAB>&       a_phi,
00435                        const LevelData<EBCellFAB>* const a_phiCoarse,
00436                        const bool&                       a_homogeneousPhysBC,
00437                        const bool&                       a_homogeneousCFBC);
00438 
00439   /// virtual function called by LevelTGA
00440   virtual void applyOpNoBoundary(LevelData<EBCellFAB>&        a_opPhi,
00441                                  const LevelData<EBCellFAB>&  a_phi);
00442 
00443   ///
00444   /**
00445      get EB flux from leveldata if we have one.  otherwise use m_ebBC
00446      ignore input data in the case of homogeneous Phys BC
00447   */
00448   virtual void
00449   applyOp(LevelData<EBCellFAB>&                    a_opPhi,
00450           const LevelData<EBCellFAB>&              a_phi,
00451           const LevelData<EBCellFAB>* const        a_phiCoar,
00452           const bool&                              a_homogeneousPhysBC,
00453           const bool&                              a_homogeneousCFBC,
00454           const LevelData<BaseIVFAB<Real> >* const a_ebFluxBCLD);
00455 
00456   virtual void
00457   applyOp(LevelData<EBCellFAB>&                    a_opPhi,
00458           const LevelData<EBCellFAB>&              a_phi,
00459           const LevelData<EBCellFAB>* const        a_phiCoar,
00460           DataIterator&                            a_dit,
00461           const bool&                              a_homogeneousPhysBC,
00462           const bool&                              a_homogeneousCFBC,
00463           const LevelData<BaseIVFAB<Real> >* const a_ebFluxBCLD);
00464 
00465   virtual void
00466   GSColorAllRegular(BaseFab<Real>&                    a_phi,
00467                     const BaseFab<Real>&              a_rhs,
00468                     const int&                        a_icolor,
00469                     const Real&                       a_weight,
00470                     const bool&                       a_homogeneousPhysBC,
00471                     const DataIndex&                  a_dit);
00472 
00473   virtual void
00474   GSColorAllIrregular(EBCellFAB&                        a_phi,
00475                       const EBCellFAB&                  a_rhs,
00476                       const int&                        a_icolor,
00477                       const bool&                       a_homogeneousPhysBC,
00478                       const DataIndex&                  a_dit);
00479 
00480   virtual void
00481   GSColorAllRegularClone(LevelData<EBCellFAB>&       a_phi,
00482                          const LevelData<EBCellFAB>& a_phiOld,
00483                          const LevelData<EBCellFAB>& a_rhs,
00484                          const int&                  a_icolor,
00485                          const Real&                 a_weight,
00486                          const bool&                 a_homogeneousPhysBC);
00487 
00488   virtual void
00489   GSColorAllIrregularClone(LevelData<EBCellFAB>&       a_phi,
00490                            const LevelData<EBCellFAB>& a_phiOld,
00491                            const LevelData<EBCellFAB>& a_rhs,
00492                            const int&                  a_icolor,
00493                            const bool&                 a_homogeneousPhysBC);
00494 
00495   ///apply domainflux in multivariable mode
00496   void
00497   mvApplyDomainFlux(BaseFab<Real> & a_phiFAB,
00498                     const Box&       a_grid,
00499                     const DataIndex& a_dit);
00500   /***/
00501   ///for taking the multi-variable laplacian (think source term in viscous flow)
00502   /// evaluates beta*laplacian (with inhomogeneous bcs).
00503   /***/
00504   void 
00505   mvBetaLaplacianGrid(EBCellFAB         &              a_lph,
00506                       const EBCellFAB   &              a_phi,
00507                       const DataIndex   &              a_dit);
00508   ///
00509   /**
00510      this is the linearop function.  CFBC is set to homogeneous.  phic is null
00511   */
00512   virtual void applyOp(LevelData<EBCellFAB>&             a_opPhi,
00513                        const LevelData<EBCellFAB>&       a_phi,
00514                        bool                              a_homogeneousPhysBC);
00515 
00516   /// no exchange of cf interp
00517   void
00518   applyOpNoCFBCs(LevelData<EBCellFAB>&                    a_opPhi,
00519                  const LevelData<EBCellFAB>&              a_phi,
00520                  const LevelData<EBCellFAB>* const        a_phiCoar,
00521                  DataIterator&                            a_dit,
00522                  const bool&                              a_homogeneousPhysBC,
00523                  const bool&                              a_homogeneousCFBC,
00524                  const LevelData<BaseIVFAB<Real> >* const a_ebFluxBCLD //only non null in multifluid
00525                  );
00526 
00527   ///
00528   /**
00529    */
00530   virtual void create(LevelData<EBCellFAB>&       a_lhs,
00531                       const LevelData<EBCellFAB>& a_rhs);
00532 
00533   ///
00534   virtual void createCoarsened(LevelData<EBCellFAB>&       a_lhs,
00535                                const LevelData<EBCellFAB>& a_rhs,
00536                                const int&                  a_refRat);
00537 
00538   virtual void buildCopier(Copier                    & a_copier, 
00539                            const LevelData<EBCellFAB>& a_lhs, 
00540                            const LevelData<EBCellFAB>& a_rhs);
00541 
00542 
00543   virtual void assignCopier(LevelData<EBCellFAB>      & a_lhs, 
00544                             const LevelData<EBCellFAB>& a_rhs, 
00545                             const Copier              & a_copier);
00546 
00547   Real
00548   AMRNorm(const LevelData<EBCellFAB>& a_coarResid,
00549           const LevelData<EBCellFAB>& a_fineResid,
00550           const int& a_refRat,
00551           const int& a_ord);
00552 
00553   ///
00554   /**
00555    */
00556   virtual void assign(LevelData<EBCellFAB>&       a_lhs,
00557                       const LevelData<EBCellFAB>& a_rhs);
00558 
00559   ///copier definition was killing us.
00560   virtual void 
00561   assignLocal(LevelData<EBCellFAB>&       a_lhs,
00562               const LevelData<EBCellFAB>& a_rhs);
00563 
00564 
00565   ///
00566   /**
00567    */
00568   virtual Real dotProduct(const LevelData<EBCellFAB>& a_1,
00569                           const LevelData<EBCellFAB>& a_2);
00570 
00571   ///
00572   /**
00573    */
00574   virtual void incr(LevelData<EBCellFAB>&       a_lhs,
00575                     const LevelData<EBCellFAB>& a_x,
00576                     Real                        a_scale);
00577 
00578   ///
00579   /**
00580    */
00581   virtual void axby(LevelData<EBCellFAB>&       a_lhs,
00582                     const LevelData<EBCellFAB>& a_x,
00583                     const LevelData<EBCellFAB>& a_y,
00584                     Real                        a_a,
00585                     Real                        a_b);
00586 
00587   ///
00588   /**
00589    */
00590   virtual void scale(LevelData<EBCellFAB>& a_lhs,
00591                      const Real&           a_scale);
00592 
00593   ///
00594   /**
00595    */
00596   virtual Real norm(const LevelData<EBCellFAB>& a_rhs,
00597                     int                         a_ord);
00598 
00599   ///
00600   /**
00601    */
00602   virtual Real localMaxNorm(const LevelData<EBCellFAB>& a_rhs);
00603 
00604   ///
00605   /**
00606    */
00607   virtual void setToZero(LevelData<EBCellFAB>& a_lhs);
00608 
00609   ///
00610   /**
00611    */
00612   virtual void setVal(LevelData<EBCellFAB>& a_lhs, const Real& a_value);
00613 
00614   ///
00615   /**
00616    */
00617   virtual void createCoarser(LevelData<EBCellFAB>&       a_coarse,
00618                              const LevelData<EBCellFAB>& a_fine,
00619                              bool                        a_ghosted);
00620 
00621   ///
00622   /**
00623    */
00624   virtual void relax(LevelData<EBCellFAB>&       a_e,
00625                      const LevelData<EBCellFAB>& a_residual,
00626                      int                         a_iterations);
00627 
00628   ///
00629   /**
00630      Calculate restricted residual:
00631      a_resCoarse[2h] = I[h->2h] (a_rhsFine[h] - L[h](a_phiFine[h]))
00632   */
00633   virtual void restrictResidual(LevelData<EBCellFAB>&       a_resCoarse,
00634                                 LevelData<EBCellFAB>&       a_phiFine,
00635                                 const LevelData<EBCellFAB>& a_rhsFine);
00636 
00637   ///
00638   /**
00639      Correct the fine solution based on coarse correction:
00640      a_phiThisLevel += I[2h->h] (a_correctCoarse)
00641   */
00642   virtual void prolongIncrement(LevelData<EBCellFAB>&       a_phiThisLevel,
00643                                 const LevelData<EBCellFAB>& a_correctCoarse);
00644 
00645   ///
00646   /** Refinement ratio between this level and coarser level.
00647       Returns 1 when there are no coarser AMRLevelOp objects */
00648   virtual int refToCoarser();
00649 
00650   ///
00651   /** Refinement ratio between this level and coarser level.
00652       Returns 1 when there are no coarser AMRLevelOp objects */
00653   virtual int refToFiner();
00654 
00655   ///
00656   /** a_residual = a_rhs - L(a_phi, a_phiFine, a_phiCoarse) */
00657   virtual void AMRResidual(LevelData<EBCellFAB>& a_residual,
00658                            const LevelData<EBCellFAB>& a_phiFine,
00659                            const LevelData<EBCellFAB>& a_phi,
00660                            const LevelData<EBCellFAB>& a_phiCoarse,
00661                            const LevelData<EBCellFAB>& a_rhs,
00662                            bool a_homogeneousBC,
00663                            AMRLevelOp<LevelData<EBCellFAB> >* a_finerOp);
00664 
00665   ///
00666   /** a_residual = a_rhs - L(a_phi, a_phiCoarse)  */
00667   virtual void AMRResidualNF(LevelData<EBCellFAB>& a_residual,
00668                              const LevelData<EBCellFAB>& a_phi,
00669                              const LevelData<EBCellFAB>& a_phiCoarse,
00670                              const LevelData<EBCellFAB>& a_rhs,
00671                              bool a_homogeneousBC);
00672 
00673 
00674   ///
00675   /** a_residual = a_rhs - L(a_phi, a_phiFine, a_phiCoarse) */
00676   virtual void AMROperator(LevelData<EBCellFAB>& a_LofPhi,
00677                            const LevelData<EBCellFAB>& a_phiFine,
00678                            const LevelData<EBCellFAB>& a_phi,
00679                            const LevelData<EBCellFAB>& a_phiCoarse,
00680                            bool a_homogeneousBC,
00681                            AMRLevelOp<LevelData<EBCellFAB> >* a_finerOp);
00682 
00683   ///
00684   /** a_residual = a_rhs - L(a_phi, a_phiCoarse)  */
00685   virtual void AMROperatorNF(LevelData<EBCellFAB>& a_LofPhi,
00686                              const LevelData<EBCellFAB>& a_phi,
00687                              const LevelData<EBCellFAB>& a_phiCoarse,
00688                              bool a_homogeneousBC);
00689 
00690   ///
00691   /** a_resCoarse = I[h-2h] (a_residual - L(a_correction, a_coarseCorrection)) */
00692   virtual void AMRRestrict(LevelData<EBCellFAB>& a_resCoarse,
00693                            const LevelData<EBCellFAB>& a_residual,
00694                            const LevelData<EBCellFAB>& a_correction,
00695                            const LevelData<EBCellFAB>& a_coarseCorrection, 
00696                            bool a_skip_res = false );
00697 
00698   ///
00699   /** a_correction += I[2h->h](a_coarseCorrection) */
00700   virtual void AMRProlong(LevelData<EBCellFAB>&       a_correction,
00701                           const LevelData<EBCellFAB>& a_coarseCorrection);
00702 
00703   ///
00704   /** a_residual = a_residual - L(a_correction, a_coarseCorrection) */
00705   virtual void AMRUpdateResidual(LevelData<EBCellFAB>&       a_residual,
00706                                  const LevelData<EBCellFAB>& a_correction,
00707                                  const LevelData<EBCellFAB>& a_coarseCorrection);
00708 
00709   ///
00710   /** a_residual = a_residual - L(a_correction, a_coarseCorrection)
00711       used in multifluid */
00712   void AMRUpdateResidual(LevelData<EBCellFAB>&       a_residual,
00713                          const LevelData<EBCellFAB>& a_correction,
00714                          const LevelData<EBCellFAB>& a_coarseCorrection,
00715                          const LevelData<BaseIVFAB<Real> >* const a_ebFluxBCLD);
00716 
00717 
00718   void reflux(LevelData<EBCellFAB>& a_residual,
00719               const LevelData<EBCellFAB>& a_phiFine,
00720               const LevelData<EBCellFAB>& a_phi,
00721               AMRLevelOp<LevelData<EBCellFAB> >* a_finerOp);
00722 
00723 
00724   void fast_reflux(LevelData<EBCellFAB>& a_residual,
00725                   const LevelData<EBCellFAB>& a_phiFine,
00726                   const LevelData<EBCellFAB>& a_phi,
00727                   AMRLevelOp<LevelData<EBCellFAB> >* a_finerOp);
00728 
00729   void setEBBC(const RefCountedPtr<BaseEBBC>&      a_ebBC);
00730 
00731   //for debugging other operators
00732   void slowGSRBColor(LevelData<EBCellFAB>&       a_phi,
00733                      const LevelData<EBCellFAB>& a_lph,
00734                      const LevelData<EBCellFAB>& a_rhs,
00735                      const IntVect&              a_color);
00736 
00737   //for debugging other operators
00738   void levelSlowRelax(LevelData<EBCellFAB>&       a_phi,
00739                       const LevelData<EBCellFAB>& a_rhs);
00740 
00741   void levelJacobi(LevelData<EBCellFAB>&       a_phi,
00742                    const LevelData<EBCellFAB>& a_rhs,
00743                    int                         a_iterations);
00744 
00745   void levelMultiColorGS(LevelData<EBCellFAB>&       a_phi,
00746                          const LevelData<EBCellFAB>& a_rhs);
00747 
00748   void levelMultiColorGS(LevelData<EBCellFAB>&       a_phi,
00749                          const LevelData<EBCellFAB>& a_resid,
00750                          const IntVect& color);
00751 
00752   void colorGS(LevelData<EBCellFAB>&       a_phi,
00753                const LevelData<EBCellFAB>& a_rhs,
00754                const int&                  a_icolor);
00755 
00756   void levelGSRB(LevelData<EBCellFAB>&       a_phi,
00757                  const LevelData<EBCellFAB>& a_rhs);
00758 
00759   void levelGSRB(LevelData<EBCellFAB>&       a_phi,
00760                  const LevelData<EBCellFAB>& a_rhs,
00761                  const int                   a_color);
00762 
00763   void levelMultiColorGSClone(LevelData<EBCellFAB>&       a_phi,
00764                               const LevelData<EBCellFAB>& a_rhs);
00765 
00766   void colorGSClone(LevelData<EBCellFAB>&       a_phi,
00767                     const LevelData<EBCellFAB>& a_phiOld,
00768                     const LevelData<EBCellFAB>& a_rhs,
00769                     const int&                  a_icolor);
00770 
00771   static void doLazyRelax(bool a_doLazyRelax);
00772   static void doEBEllipticLoadBalance(bool a_doEBEllipticLoadBalance);
00773   static void areaFracWeighted(bool a_areaFracWeighted);
00774 
00775   static void  getDivFStencil(VoFStencil&      a_vofStencil,
00776                               const VolIndex&  a_vof,
00777                               const EBISBox&   a_ebisBox,
00778                               const IntVectSet& a_cfivs,
00779                               const RealVect&  a_dx,
00780                               bool doFaceInterp = true);
00781 
00782   void  getDivFStencil(VoFStencil&      a_vofStencil,
00783                        const VolIndex&  a_vof,
00784                        const DataIndex& a_dit,
00785                        bool doFaceInterp);
00786 
00787   //CP recently added
00788   void getVoFStencil(LayoutData<BaseIVFAB<VoFStencil> > const*& a_vofStencil)
00789   {
00790     //I'm defining const* to make sure the content in m_opStencil is not modified
00791     a_vofStencil = &m_opStencil;
00792   }
00793 
00794   void getAlphaDiagWeight(LayoutData<BaseIVFAB<Real> > const*& a_alphaDiagWeight)
00795   {
00796     a_alphaDiagWeight = &m_alphaDiagWeight;
00797   }
00798 
00799   void getAlphaBeta(Real& a_alpha, Real& a_beta)
00800   {
00801     a_alpha = m_alpha;
00802     a_beta  = m_beta;
00803   }
00804   //CP recently added
00805   const RefCountedPtr<BaseDomainBC> getDomainBC()
00806   {
00807     //I'm defining const* to make sure the content in m_opStencil is not modified
00808     return m_domainBC;
00809   }
00810 
00811   void setListValue(const LevelData<EBCellFAB>& a_data, Real a_value);
00812 
00813   void setListValue(EBCellFAB& a_data, const Vector<VolIndex>& a_setList, Real a_value);
00814 
00815   void setRhsSetList(const LayoutData<Vector<VolIndex> >& a_list);
00816 
00817   LayoutData<Vector<VolIndex> >& getRhsSetList()
00818   {    
00819     return m_rhsSetList;
00820   }
00821 
00822   const LayoutData<Vector<VolIndex> >& getListValue()
00823   {    
00824     return m_rhsSetList;
00825   }
00826 
00827   static void getFluxStencil(VoFStencil&       a_fluxStencil,
00828                              const FaceIndex&  a_face,
00829                              const EBISBox&    a_ebisBox,
00830                              const IntVectSet& a_cfivs,
00831                              const RealVect&   a_dx,
00832                              bool a_doFaceInterp);
00833 
00834   static void getFaceCenteredFluxStencil(VoFStencil&      a_fluxStencil,
00835                                          const FaceIndex& a_face,
00836                                          const RealVect&  a_dx);
00837 
00838   void applyCFBCs(LevelData<EBCellFAB>&             a_phi,
00839                   const LevelData<EBCellFAB>* const a_phiCoarse,
00840                   bool a_homogeneousCFBC,
00841                   bool a_doOnlyRegularInterp = false);
00842 
00843   ///
00844   /**
00845      This function computes: a_lhs = (1/diagonal)*a_rhs
00846      It is used to initialize the preconditioner, and by
00847      MFPoissonOp::levelJacobi.
00848      Consider using one of the level Gauss-Seidel methods
00849      instead of monkeying with this.
00850   */
00851   void getInvDiagRHS(LevelData<EBCellFAB>&       a_lhs,
00852                      const LevelData<EBCellFAB>& a_rhs);
00853 
00854   static int                      s_numComps;
00855   static int                      s_whichComp;
00856 
00857 protected:
00858 
00859   LevelData<EBCellFAB> m_resThisLevel;
00860   void defineMGObjects(const EBLevelGrid& a_eblgCoarMG);
00861   void defineWithCoarser(const EBLevelGrid& a_eblgCoar, const int& a_refToCoar);
00862   void defineWithFiner(const EBLevelGrid& a_eblgFine,
00863                        const int& a_refToFine);
00864 
00865   static bool                     s_turnOffBCs;
00866   static bool                     s_doEBEllipticLoadBalance;
00867   static bool                     s_areaFracWeighted;
00868   void defineStencils();
00869   void defineEBCFStencils();
00870 
00871   int                             m_testRef;
00872 
00873   const IntVect                   m_ghostCellsPhi;
00874   const IntVect                   m_ghostCellsRHS;
00875 
00876   RefCountedPtr<EBQuadCFInterp>   m_quadCFIWithCoar;
00877 
00878   EBLevelGrid                     m_eblg;
00879   EBLevelGrid                     m_eblgFine;
00880   EBLevelGrid                     m_eblgCoar;
00881   EBLevelGrid                     m_eblgCoarMG;
00882   EBLevelGrid                     m_eblgCoarsenedFine;
00883 
00884   RefCountedPtr<BaseDomainBC>     m_domainBC;
00885   RefCountedPtr<BaseEBBC>         m_ebBC;
00886   LayoutData<Vector<VolIndex> >   m_rhsSetList;
00887 
00888   RealVect                        m_dxFine;
00889   RealVect                        m_dx;
00890   RealVect                        m_dxCoar;
00891 
00892   RealVect                        m_invDx;
00893   RealVect                        m_invDx2;
00894   Real                            m_dxScale;
00895   Real                            m_alpha;
00896   Real                            m_aCoef;
00897   Real                            m_beta;
00898   Real                            m_bCoef;
00899   static Real                     s_time;
00900   static bool                     s_doLazyRelax;
00901   static bool                     s_doInconsistentRelax;
00902   static bool                     s_doTrimEdges;
00903   static bool                     s_doSetListValueResid;  // if this variable is set to true, it will apply setListValue in residual()
00904   RealVect                        m_origin;
00905   int                             m_refToFine;
00906   int                             m_refToCoar;
00907   bool                            m_hasFine;
00908   bool                            m_hasInterpAve;
00909   bool                            m_hasCoar;
00910   int                             m_numPreCondIters;
00911   int                             m_relaxType;
00912 
00913   bool                            m_hasEBCF;
00914 
00915   Copier                          m_exchangeCopier;
00916   //restriction object
00917   EBMGAverage                    m_ebAverage;
00918   //prolongation object
00919   EBMGInterp                     m_ebInterp;
00920 
00921   //stencils for operator evaluation
00922   LayoutData<BaseIVFAB<VoFStencil> > m_opStencil;
00923   LayoutData<RefCountedPtr<EBSTENCIL_T> >  m_opEBStencil;
00924   LayoutData<RefCountedPtr<EBSTENCIL_T> >  m_opEBStencilInhomDomLo[SpaceDim];
00925   LayoutData<RefCountedPtr<EBSTENCIL_T> >  m_opEBStencilInhomDomHi[SpaceDim];
00926   LayoutData<RefCountedPtr<EBSTENCIL_T> >  m_colorEBStencil[EBAMRPO_NUMSTEN];
00927   LayoutData<RefCountedPtr<EBSTENCIL_T> >  m_colorEBStencilDomLo[EBAMRPO_NUMSTEN][SpaceDim];
00928   LayoutData<RefCountedPtr<EBSTENCIL_T> >  m_colorEBStencilDomHi[EBAMRPO_NUMSTEN][SpaceDim];
00929   LayoutData<RefCountedPtr<EBSTENCIL_T> >  m_invDiagEBStencil;
00930   //weights that get multiplied by alpha
00931   LayoutData<BaseIVFAB<Real> >       m_alphaDiagWeight;
00932   //weights that get multiplied by beta
00933   LayoutData<BaseIVFAB<Real> >       m_betaDiagWeight;
00934   //constant for using in place of weight in EBStencil->apply for inhom dom bcs
00935   LayoutData<BaseIVFAB<Real> >       m_one;
00936 
00937   //cache the irreg vofiterator
00938   LayoutData<VoFIterator >     m_vofItIrreg;
00939   LayoutData<VoFIterator >     m_vofItIrregColor[EBAMRPO_NUMSTEN];
00940 
00941   LayoutData<VoFIterator >     m_vofItIrregDomLo[SpaceDim];
00942   LayoutData<VoFIterator >     m_vofItIrregDomHi[SpaceDim];
00943 
00944   LayoutData<VoFIterator >     m_vofItIrregColorDomLo[EBAMRPO_NUMSTEN][SpaceDim];
00945   LayoutData<VoFIterator >     m_vofItIrregColorDomHi[EBAMRPO_NUMSTEN][SpaceDim];
00946   Vector<Vector<RefCountedPtr<LayoutData<EBCellFAB> > > > m_cacheInhomDomBCLo;
00947   Vector<Vector<RefCountedPtr<LayoutData<EBCellFAB> > > > m_cacheInhomDomBCHi;
00948   // LayoutData<BaseIVFAB<Real> > m_cacheIrregFluxDomLo[EBAMRPO_NUMSTEN][SpaceDim];
00949   // LayoutData<BaseIVFAB<Real> > m_cacheIrregFluxDomHi[EBAMRPO_NUMSTEN][SpaceDim];
00950 
00951   // Coarse-fine stencils for homogeneous CFInterp
00952   LayoutData<CFIVS> m_loCFIVS[SpaceDim];
00953   LayoutData<CFIVS> m_hiCFIVS[SpaceDim];
00954 
00955   //flux register with finer level
00956   EBFastFR       m_fastFR;
00957 
00958   //stuff to make EBCF go faster
00959   LayoutData< Vector<FaceIndex>  >   m_faceitCoar[2*SpaceDim];
00960   LayoutData< Vector<VoFStencil> >  m_stencilCoar[2*SpaceDim];
00961 
00962 
00963   //special mg objects for when we do not have
00964   //a coarser level or when the refinement to coarse
00965   //is greater than two
00966   //flag for when we need special MG objects
00967   bool                        m_hasMGObjects;
00968   bool                        m_layoutChanged;
00969 
00970   Vector<IntVect> m_colors;
00971 
00972 
00973   //stuff below is only defined if m_hasMGObjects==true
00974   EBMGAverage                 m_ebAverageMG;
00975   EBMGInterp                  m_ebInterpMG;
00976   DisjointBoxLayout           m_dblCoarMG;
00977   EBISLayout                  m_ebislCoarMG;
00978   ProblemDomain               m_domainCoarMG;
00979 
00980 private:
00981   //internally useful but not for general consumption
00982   //lots of hidden assumptions and all that
00983 
00984   void fast_incrementFRCoar(const LevelData<EBCellFAB>& a_phiFine,
00985                             const LevelData<EBCellFAB>& a_phi);
00986 
00987   void fast_incrementFRFine(const LevelData<EBCellFAB>& a_phiFine,
00988                             const LevelData<EBCellFAB>& a_phi,
00989                             AMRLevelOp<LevelData<EBCellFAB> >* a_finerOp);
00990 
00991   ///this one gets called by base level tga (called by the public getFlux)
00992   void getFlux(EBFaceFAB&                    a_flux,
00993                const EBCellFAB&              a_phi,
00994                const Box&                    a_ghostedBox,
00995                const Box&                    a_fabBox,
00996                const ProblemDomain&          a_domainBox,
00997                const EBISBox&                a_ebisBox,
00998                const RealVect&               a_dx,
00999                const int&                    a_idir);
01000 
01001   ///this one is internal (called by refluxing)
01002   void getFluxEBCF(EBFaceFAB&                    a_flux,
01003                    const EBCellFAB&              a_phi,
01004                    const Box&                    a_ghostedBox,
01005                    Vector<FaceIndex>&            a_faceitEBCF,
01006                    Vector<VoFStencil>&           a_ebcfsten,
01007                    const RealVect&               a_dx);
01008 
01009   ///this one is internal (called by refluxing)
01010   void getFluxRegO(EBFaceFAB&                    a_flux,
01011                   const EBCellFAB&              a_phi,
01012                   const Box&                    a_ghostedBox,
01013                   const RealVect&               a_dx);
01014 
01015   void getOpVoFStencil(VoFStencil&     a_stencil,
01016                        const EBISBox&  a_ebisbox,
01017                        const VolIndex& a_vof);
01018 
01019   void getOpVoFStencil(VoFStencil&             a_stencil,
01020                        const int&              a_dir,
01021                        const Vector<VolIndex>& a_allMonotoneVoFs,
01022                        const EBISBox&          a_ebisbox,
01023                        const VolIndex&         a_vof,
01024                        const bool&             a_lowOrder);
01025 
01026 
01027   void getOpFaceStencil(VoFStencil&             a_stencil,
01028                         const Vector<VolIndex>& a_allMonotoneVofs,
01029                         const EBISBox&          a_ebisbox,
01030                         const VolIndex&         a_vof,
01031                         int                     a_dir,
01032                         const Side::LoHiSide&   a_side,
01033                         const FaceIndex&        a_face,
01034                         const bool&             a_lowOrder);
01035 
01036   void levelJacobi(LevelData<EBCellFAB>&       a_phi,
01037                    const LevelData<EBCellFAB>& a_rhs);
01038 
01039   void applyHomogeneousCFBCs(LevelData<EBCellFAB>&   a_phi);
01040 
01041   void applyHomogeneousCFBCs(EBCellFAB&            a_phi,
01042                              const DataIndex&      a_datInd,
01043                              int                   a_idir,
01044                              Side::LoHiSide        a_hiorlo);
01045 
01046   Real getRadius(const FaceIndex& a_face, const RealVect& a_centroid);
01047 
01048   ///
01049   /** applyOp() on the regular cells for all directions
01050    opPhi comes in holding alpha*phi.  this adds on beta*lapl phi*/
01051   void applyOpRegularAllDirs(Box * a_loBox,
01052                              Box * a_hiBox,
01053                              int * a_hasLo,
01054                              int * a_hasHi,
01055                              Box & a_curOpPhiBox,
01056                              Box & a_curPhiBox,
01057                              int a_nComps,
01058                              BaseFab<Real> & a_curOpPhiFAB,
01059                              const BaseFab<Real> & a_curPhiFAB,
01060                              bool a_homogeneousPhysBC,
01061                              const DataIndex& a_dit,
01062                              const Real& a_beta);
01063 
01064   void applyDomainFlux(Box * a_loBox,
01065                        Box * a_hiBox,
01066                        int * a_hasLo,
01067                        int * a_hasHi,
01068                        Box & a_curPhiBox,
01069                        int a_nComps,
01070                        BaseFab<Real> & a_phiFAB,
01071                        bool a_homogeneousPhysBC,
01072                        const DataIndex& a_dit,
01073                        const Real& a_beta);
01074 
01075 private:
01076   //copy constructor and operator= disallowed for all the usual reasons
01077   EBAMRPoissonOp(const EBAMRPoissonOp& a_opin)
01078   {
01079     MayDay::Error("invalid operator");
01080   }
01081 
01082   void operator=(const EBAMRPoissonOp& a_opin)
01083   {
01084     MayDay::Error("invalid operator");
01085   }
01086 };
01087 
01088 
01089 #include "NamespaceFooter.H"
01090 #endif