Chombo + EB: EBPoissonOp.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 _EBPOISSONOP_H__
00012 #define _EBPOISSONOP_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 
00027 #include "EBStencil.H"
00028 
00029 #include "EBLevelDataOps.H"
00030 #include "BaseEBBC.H"
00031 #include "BaseDomainBC.H"
00032 #include "CFIVS.H"
00033 #include "EBFluxRegister.H"
00034 #include "EBMGAverage.H"
00035 #include "EBMGInterp.H"
00036 #include "EBCoarsen.H"
00037 #include "PolyGeom.H"
00038 #include "EBQuadCFInterp.H"
00039 #include "EBLevelGrid.H"
00040 #include "NamespaceHeader.H"
00041 
00042 
00043 #if CH_SPACEDIM==2
00044 #define EBPO_NUMSTEN 4
00045 #elif CH_SPACEDIM==3
00046 #define EBPO_NUMSTEN 8
00047 #else
00048 void THIS_IS_AN_ERROR_MESSAGE(void)
00049 {
00050   THIS_WILL_ONLY_COMPILE_WHEN_CH_SPACEDIM_IS_2_OR_3;
00051 }
00052 #endif
00054 
00057 class EBPoissonOp: public MGLevelOp<LevelData<EBCellFAB> >
00058 {
00059 public:
00061   ~EBPoissonOp();
00062 
00064   EBPoissonOp();
00065 
00067 
00087   EBPoissonOp(   const EBLevelGrid &                  a_eblg,
00088                  const EBLevelGrid &                  a_eblgCoarMG,
00089                  const RefCountedPtr<BaseDomainBC>&   a_domainBC,
00090                  const RefCountedPtr<BaseEBBC>&       a_ebBC,
00091                  const RealVect&                      a_dx,
00092                  const RealVect&                      a_origin,
00093                  const bool&                          a_hasMGObjects,
00094                  const int&                           a_numPreCondIters,
00095                  const int&                           a_relaxType,
00096                  const int&                           a_orderEB,
00097                  const Real&                          a_alpha,
00098                  const Real&                          a_beta,
00099                  const IntVect&                       a_ghostCellsPhi,
00100                  const IntVect&                       a_ghostCellsRHS);
00101 
00102   //MGOp operations.  no finer or coarser
00103 
00105 
00107   virtual void residual(LevelData<EBCellFAB>&       a_residual,
00108                         const LevelData<EBCellFAB>& a_phi,
00109                         const LevelData<EBCellFAB>& a_rhs,
00110                         bool                        a_homogeneousPhysBC=false);
00111 
00113 
00115   virtual void preCond(LevelData<EBCellFAB>&       a_opPhi,
00116                        const LevelData<EBCellFAB>& a_phi);
00117 
00118 
00119   virtual void
00120   GSColorAllRegular(LevelData<EBCellFAB>&             a_phi,
00121                     const LevelData<EBCellFAB>&       a_rhs,
00122                     const IntVect&                    a_color,
00123                     const Real&                       a_weight,
00124                     const bool&                       a_homogeneousPhysBC);
00125 
00126   virtual void
00127   GSColorAllIrregular(LevelData<EBCellFAB>&             a_phi,
00128                       const LevelData<EBCellFAB>&       a_rhs,
00129                       const IntVect&                    a_color,
00130                       const bool&                       a_homogeneousPhysBC,
00131                       int icolor);
00132 
00133 
00135 
00138   virtual void applyOp(LevelData<EBCellFAB>&             a_opPhi,
00139                        const LevelData<EBCellFAB>&       a_phi,
00140                        bool                              a_homogeneousPhysBC,
00141                        DataIterator& a_dit,
00142                        bool do_exchange = true);
00143 
00144   virtual void applyOp(LevelData<EBCellFAB>&             a_opPhi,
00145                        const LevelData<EBCellFAB>&       a_phi,
00146                        bool                              a_homogeneousPhysBC)
00147   {
00148     DataIterator dit = a_opPhi.dataIterator();
00149     applyOp(a_opPhi, a_phi, a_homogeneousPhysBC, dit, true);
00150   }
00151 
00153 
00155   virtual void create(LevelData<EBCellFAB>&       a_lhs,
00156                       const LevelData<EBCellFAB>& a_rhs);
00157 
00159 
00161   virtual void assign(LevelData<EBCellFAB>&       a_lhs,
00162                       const LevelData<EBCellFAB>& a_rhs);
00163 
00165 
00167   virtual Real dotProduct(const LevelData<EBCellFAB>& a_1,
00168                           const LevelData<EBCellFAB>& a_2);
00169 
00171 
00173   virtual void incr(LevelData<EBCellFAB>&       a_lhs,
00174                     const LevelData<EBCellFAB>& a_x,
00175                     Real                        a_scale);
00176 
00178 
00180   virtual void axby(LevelData<EBCellFAB>&       a_lhs,
00181                     const LevelData<EBCellFAB>& a_x,
00182                     const LevelData<EBCellFAB>& a_y,
00183                     Real                        a_a,
00184                     Real                        a_b);
00185 
00187 
00189   virtual void scale(LevelData<EBCellFAB>& a_lhs,
00190                      const Real&           a_scale);
00191 
00193 
00195   virtual Real norm(const LevelData<EBCellFAB>& a_rhs,
00196                     int                         a_ord);
00197 
00199 
00201   virtual void setToZero(LevelData<EBCellFAB>& a_lhs);
00202 
00204 
00206   virtual void setVal(LevelData<EBCellFAB>& a_lhs, const Real& a_value);
00207 
00209 
00211   virtual void createCoarser(LevelData<EBCellFAB>&       a_coarse,
00212                              const LevelData<EBCellFAB>& a_fine,
00213                              bool                        a_ghosted);
00214 
00216 
00218   virtual void relax(LevelData<EBCellFAB>&       a_e,
00219                      const LevelData<EBCellFAB>& a_residual,
00220                      int                         a_iterations);
00221 
00223 
00227   virtual void restrictResidual(LevelData<EBCellFAB>&       a_resCoarse,
00228                                 LevelData<EBCellFAB>&       a_phiFine,
00229                                 const LevelData<EBCellFAB>& a_rhsFine);
00230 
00232 
00236   virtual void prolongIncrement(LevelData<EBCellFAB>&       a_phiThisLevel,
00237                                 const LevelData<EBCellFAB>& a_correctCoarse);
00238 
00239 
00240   static bool nextColor(IntVect& color,
00241                         const IntVect& limit);
00242 
00243   static int getIndex(const IntVect& a_color);
00244 
00245   static IntVect getColor(const int& a_icolor);
00246 
00247   void levelMulticolorGS(LevelData<EBCellFAB>&       a_phi,
00248                          const LevelData<EBCellFAB>& a_rhs);
00249 
00250   void colorGS(LevelData<EBCellFAB>&       a_phi,
00251                const LevelData<EBCellFAB>& a_rhs,
00252                const IntVect& color, int icolor);
00253 
00254 
00255 
00256 protected:
00257 
00258   void defineStencils();
00259 
00260   const IntVect                   m_ghostCellsPhi;
00261   const IntVect                   m_ghostCellsRHS;
00262   Vector<IntVect> m_colors;
00263 
00264 
00265   EBLevelGrid                     m_eblg;
00266   EBLevelGrid                     m_eblgCoarMG;
00267 
00268   RefCountedPtr<BaseDomainBC>     m_domainBC;
00269   RefCountedPtr<BaseEBBC>         m_ebBC;
00270 
00271   RealVect                        m_dx;
00272 
00273   RealVect                        m_invDx;
00274   RealVect                        m_invDx2;
00275   Real                            m_dxScale;
00276   Real                            m_alpha;
00277   Real                            m_beta;
00278   Real                            m_time;
00279   RealVect                        m_origin;
00280   int                             m_orderEB;
00281   int                             m_numPreCondIters;
00282   int                             m_relaxType;
00283 
00284   Copier                          m_exchangeCopier;
00285 
00286   //stencils for operator evaluation
00287   LayoutData<RefCountedPtr<EBStencil> >  m_opEBStencil;
00288   //stencils for level Jacobi evaluation
00289   LayoutData<RefCountedPtr<EBStencil> >  m_jacobiEBStencil;
00290   LayoutData<RefCountedPtr<EBStencil> >  m_rhsJacobiEBStencil;
00291   //stencils for operator evaluation on gauss-seidel colors
00292   LayoutData<RefCountedPtr<EBStencil> >  m_colorEBStencil[EBPO_NUMSTEN];
00293   //stencils for lambda*rhs addition to opphi
00294   LayoutData<RefCountedPtr<EBStencil> >  m_rhsColorEBStencil[EBPO_NUMSTEN];
00295   //weights that get multiplied by alpha
00296   LayoutData<BaseIVFAB<Real> >       m_alphaDiagWeight;
00297   //weights that get multiplied by beta
00298   LayoutData<BaseIVFAB<Real> >       m_betaDiagWeight;
00299 
00300   //cache the irreg vofiterator
00301   LayoutData<VoFIterator > m_vofItIrreg;
00302   LayoutData<VoFIterator > m_vofItIrregColor[EBPO_NUMSTEN];
00303 
00304   LayoutData<VoFIterator > m_vofItIrregDomLo[SpaceDim];
00305   LayoutData<VoFIterator > m_vofItIrregDomHi[SpaceDim];
00306 
00307   LayoutData<VoFIterator > m_vofItIrregColorDomLo[EBPO_NUMSTEN][SpaceDim];
00308   LayoutData<VoFIterator > m_vofItIrregColorDomHi[EBPO_NUMSTEN][SpaceDim];
00309 
00310   bool                        m_hasMGObjects;
00311   //stuff below is only defined if m_hasMGObjects==true
00312   EBMGAverage                 m_ebAverageMG;
00313   EBMGInterp                  m_ebInterpMG;
00314   DisjointBoxLayout           m_dblCoarMG;
00315   EBISLayout                  m_ebislCoarMG;
00316   ProblemDomain               m_domainCoarMG;
00317 
00318 private:
00319 
00320 
00321   //internally useful but not for general consumption
00322   //lots of hidden assumptions and all that
00323 
00324   void getOpVoFStencil(VoFStencil&     a_stencil,
00325                        const EBISBox&  a_ebisbox,
00326                        const VolIndex& a_vof);
00327 
00328   void getOpVoFStencil(VoFStencil&             a_stencil,
00329                        const int&              a_dir,
00330                        const Vector<VolIndex>& a_allMonotoneVoFs,
00331                        const EBISBox&          a_ebisbox,
00332                        const VolIndex&         a_vof,
00333                        const bool&             a_lowOrder);
00334 
00335 
00336   void getOpFaceStencil(VoFStencil&             a_stencil,
00337                         const Vector<VolIndex>& a_allMonotoneVofs,
00338                         const EBISBox&          a_ebisbox,
00339                         const VolIndex&         a_vof,
00340                         int                     a_dir,
00341                         const Side::LoHiSide&   a_side,
00342                         const FaceIndex&        a_face,
00343                         const bool&             a_lowOrder);
00344 
00345 
00346   void levelJacobi(LevelData<EBCellFAB>&       a_phi,
00347                    const LevelData<EBCellFAB>& a_rhs);
00348 
00350 
00352   void applyOpRegular( int idir,
00353                        Box * a_loBox,
00354                        Box * a_hiBox,
00355                        int * a_hasLo,
00356                        int * a_hasHi,
00357                        Box & a_curOpPhiBox,
00358                        Box & a_curPhiBox,
00359                        int a_nComps,
00360                        BaseFab<Real> & a_curOpPhiFAB,
00361                        const BaseFab<Real> & a_curPhiFAB,
00362                        bool a_homogeneousPhysBC,
00363                        const DataIndex& a_dit,
00364                        const Real& a_beta);
00365 
00367 
00369   void applyOpRegularAllDirs(Box * a_loBox,
00370                              Box * a_hiBox,
00371                              int * a_hasLo,
00372                              int * a_hasHi,
00373                              Box & a_curOpPhiBox,
00374                              Box & a_curPhiBox,
00375                              int a_nComps,
00376                              BaseFab<Real> & a_curOpPhiFAB,
00377                              const BaseFab<Real> & a_curPhiFAB,
00378                              bool a_homogeneousPhysBC,
00379                              const DataIndex& a_dit,
00380                              const Real& a_beta);
00381 
00382   void applyDomainFlux(Box * a_loBox,
00383                        Box * a_hiBox,
00384                        int * a_hasLo,
00385                        int * a_hasHi,
00386                        Box & a_curPhiBox,
00387                        int a_nComps,
00388                        BaseFab<Real> & a_phiFAB,
00389                        bool a_homogeneousPhysBC,
00390                        const DataIndex& a_dit,
00391                        const Real& a_beta);
00392 
00394 
00398   void getInvDiagRHS(LevelData<EBCellFAB>&       a_lhs,
00399                      const LevelData<EBCellFAB>& a_rhs);
00400 
00401 private:
00402   //copy constructor and operator= disallowed for all the usual reasons
00403   EBPoissonOp(const EBPoissonOp& a_opin)
00404   {
00405     MayDay::Error("invalid operator");
00406   }
00407 
00408   void operator=(const EBPoissonOp& a_opin)
00409   {
00410     MayDay::Error("invalid operator");
00411   }
00412 };
00413 
00414 
00415 #include "NamespaceFooter.H"
00416 #endif