Chombo + EB: slowEBCO.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 _SLOWEBCO_H_
00012 #define _SLOWEBCO_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 #include "EBStencil.H"
00027 
00028 #include "EBLevelDataOps.H"
00029 #include "BaseEBBC.H"
00030 #include "ConductivityBaseDomainBC.H"
00031 #include "CFIVS.H"
00032 #include "EBFluxRegister.H"
00033 #include "EBFastFR.H"
00034 #include "EBMGAverage.H"
00035 #include "EBMGInterp.H"
00036 #include "PolyGeom.H"
00037 #include "EBQuadCFInterp.H"
00038 #include "EBLevelGrid.H"
00039 #include "AMRTGA.H"
00040 #include "AMRPoissonOp.H"
00041 #include "CFRegion.H"
00042 #include "NamespaceHeader.H"
00043 
00044 
00045 ///
00046 /**
00047    Operator to solve (alpha a + beta div (b grad) )phi = rhs.   This follows the AMRLevelOp interface.
00048 */
00049 class slowEBCO: public LevelTGAHelmOp<LevelData<EBCellFAB>, EBFluxFAB >
00050 {
00051 public:
00052 
00053   Real getSafety();
00054 
00055   //for tga to reset stuff
00056   virtual void setAlphaAndBeta(const Real& a_alpha,
00057                                const Real& a_beta);
00058 
00059   //another tgaism
00060   virtual void diagonalScale(LevelData<EBCellFAB> & a_rhs,
00061                              bool a_kappaWeighted);
00062   virtual void divideByIdentityCoef(LevelData<EBCellFAB> & a_rhs);
00063 
00064   ///a leveltgaism
00065   virtual void fillGrad(const LevelData<EBCellFAB>& a_phi)
00066   {;}
00067 
00068   ///a leveltgaism
00069   virtual void getFlux(EBFluxFAB&                    a_flux,
00070                        const LevelData<EBCellFAB>&   a_data,
00071                        const Box&                    a_grid,
00072                        const DataIndex&              a_dit,
00073                        Real                          a_scale);
00074 
00075   void getFlux(EBFaceFAB&                    a_fluxCentroid,
00076                const EBCellFAB&              a_phi,
00077                const Box&                    a_ghostedBox,
00078                const Box&                    a_fabBox,
00079                const ProblemDomain&          a_domain,
00080                const EBISBox&                a_ebisBox,
00081                const Real&                   a_dx,
00082                const DataIndex&              a_datInd,
00083                const int&                    a_idir);
00084 
00085   ///
00086   virtual ~slowEBCO()
00087   {;}
00088 
00089   ///
00090   slowEBCO()
00091   {;}
00092 
00093   ///
00094   /** a_residual = a_rhs - L(a_phiFine, a_phi)   no coaser AMR level*/
00095   void AMRResidualNC(LevelData<EBCellFAB>&       a_residual,
00096                      const LevelData<EBCellFAB>& a_phiFine,
00097                      const LevelData<EBCellFAB>& a_phi,
00098                      const LevelData<EBCellFAB>& a_rhs,
00099                      bool a_homogeneousBC,
00100                      AMRLevelOp<LevelData<EBCellFAB> >* a_finerOp);
00101 
00102 
00103   ///
00104   /** apply AMR operator   no coaser AMR level*/
00105   void AMROperatorNC(LevelData<EBCellFAB>&       a_LofPhi,
00106                      const LevelData<EBCellFAB>& a_phiFine,
00107                      const LevelData<EBCellFAB>& a_phi,
00108                      bool a_homogeneousBC,
00109                      AMRLevelOp<LevelData<EBCellFAB> >* a_finerOp);
00110 
00111   ///
00112   /**
00113      If you are approaching this operator from this interface, consider backing away and using
00114      slowEBCOFactory to generate these objects.   Really.
00115      a_eblgFine,    : grid at finer  level \\
00116      a_eblg,        : grid at this  level \\
00117      a_eblgCoar,    : grid at coarser level \\
00118      a_eblgCoarMG,  : grid at intermediate multigrid level \\
00119      a_domainBC,    : domain boundary conditions at this level   \\
00120      a_ebBC:          eb boundary conditions at this level   \\
00121      a_dx:            grid spacing at this level \\
00122      a_origin:        offset to lowest corner of the domain \\
00123      a_refToFine:     refinement ratio to finer level \\
00124      a_refToCoar:     refinement ratio to coarser level \\
00125      a_hasFiner:      true if there is a finer AMR level, false otherwise. \\
00126      a_hasCoarser:    true  if there is a coarser AMR level. \\
00127      a_hasCoarserMG:    true  if there is a coarser MultiGrid level. \\
00128      a_preCondIters:  number of iterations to do for pre-conditioning \\
00129      a_relaxType:     0 means point Jacobi, 1 is Gauss-Seidel. \\
00130      a_acoef:         coefficent of identity \\
00131      a_bcoef:          coefficient of gradient.\\
00132      a_ghostCellsPhi:  Number of ghost cells in phi, correction\\
00133      a_ghostCellsRhs:  Number of ghost cells in RHS, residual, lphi\\
00134      Ghost cell arguments are there for caching reasons.  Once you set them, an error is thrown if
00135      you send in data that does not match.
00136   */
00137   slowEBCO(const EBLevelGrid &                                  a_eblgFine,
00138                    const EBLevelGrid &                                  a_eblg,
00139                    const EBLevelGrid &                                  a_eblgCoar,
00140                    const EBLevelGrid &                                  a_eblgCoarMG,
00141                    const RefCountedPtr<EBQuadCFInterp>&                 a_quadCFI,
00142                    const RefCountedPtr<ConductivityBaseDomainBC>&       a_domainBC,
00143                    const RefCountedPtr<ConductivityBaseEBBC>&           a_ebBC,
00144                    const Real    &                                      a_dx,
00145                    const Real    &                                      a_dxCoar,
00146                    const int&                                           a_refToFine,
00147                    const int&                                           a_refToCoar,
00148                    const bool&                                          a_hasFine,
00149                    const bool&                                          a_hasCoar,
00150                    const bool&                                          a_hasMGObjects,
00151                    const bool&                                          a_layoutChanged,
00152                    const Real&                                          a_alpha,
00153                    const Real&                                          a_beta,
00154                    const RefCountedPtr<LevelData<EBCellFAB> >&          a_acoef,
00155                    const RefCountedPtr<LevelData<EBFluxFAB> >&          a_bcoef,
00156                    const RefCountedPtr<LevelData<BaseIVFAB<Real> > >&   a_bcoIrreg,
00157                    const IntVect&                                       a_ghostCellsPhi,
00158                    const IntVect&                                       a_ghostCellsRHS,
00159                    const int&                                           a_relaxType);
00160 
00161   //MGOp operations.  no finer or coarser
00162 
00163   ///
00164   /**
00165    */
00166   virtual void residual(LevelData<EBCellFAB>&       a_residual,
00167                         const LevelData<EBCellFAB>& a_phi,
00168                         const LevelData<EBCellFAB>& a_rhs,
00169                         bool                        a_homogeneousPhysBC=false);
00170 
00171   ///
00172   /**
00173    */
00174   virtual void preCond(LevelData<EBCellFAB>&       a_opPhi,
00175                        const LevelData<EBCellFAB>& a_phi);
00176 
00177   ///
00178   /**
00179      This function assumes that coarse-fine boundary condtions have
00180      been dealt with.
00181   */
00182   virtual void applyOp(LevelData<EBCellFAB>&             a_opPhi,
00183                        const LevelData<EBCellFAB>&       a_phi,
00184                        const LevelData<EBCellFAB>* const a_phiCoarse,
00185                        const bool&                       a_homogeneousPhysBC,
00186                        const bool&                       a_homogeneousCFBC);
00187 
00188   /// virtual function called by LevelTGA
00189   virtual void applyOpNoBoundary(LevelData<EBCellFAB>&        a_opPhi,
00190                                  const LevelData<EBCellFAB>&  a_phi);
00191 
00192   ///
00193   /**
00194      this is the linearop function.  CFBC is set to homogeneous.  phic is null
00195   */
00196   virtual void applyOp(LevelData<EBCellFAB>&             a_opPhi,
00197                        const LevelData<EBCellFAB>&       a_phi,
00198                        bool                              a_homogeneousPhysBC);
00199 
00200   ///
00201   /**
00202    */
00203   virtual void create(LevelData<EBCellFAB>&       a_lhs,
00204                       const LevelData<EBCellFAB>& a_rhs);
00205 
00206   ///
00207   virtual void createCoarsened(LevelData<EBCellFAB>&       a_lhs,
00208                                const LevelData<EBCellFAB>& a_rhs,
00209                                const int&                  a_refRat);
00210 
00211   Real
00212   AMRNorm(const LevelData<EBCellFAB>& a_coarResid,
00213           const LevelData<EBCellFAB>& a_fineResid,
00214           const int& a_refRat,
00215           const int& a_ord);
00216 
00217   ///
00218   /**
00219    */
00220   virtual void assign(LevelData<EBCellFAB>&       a_lhs,
00221                       const LevelData<EBCellFAB>& a_rhs);
00222 
00223   ///
00224   /**
00225    */
00226   virtual Real dotProduct(const LevelData<EBCellFAB>& a_1,
00227                           const LevelData<EBCellFAB>& a_2);
00228 
00229   ///
00230   /**
00231    */
00232   virtual void incr(LevelData<EBCellFAB>&       a_lhs,
00233                     const LevelData<EBCellFAB>& a_x,
00234                     Real                        a_scale);
00235 
00236   ///
00237   /**
00238    */
00239   virtual void axby(LevelData<EBCellFAB>&       a_lhs,
00240                     const LevelData<EBCellFAB>& a_x,
00241                     const LevelData<EBCellFAB>& a_y,
00242                     Real                        a_a,
00243                     Real                        a_b);
00244 
00245   ///
00246   /**
00247    */
00248   virtual void scale(LevelData<EBCellFAB>& a_lhs,
00249                      const Real&           a_scale);
00250 
00251   ///
00252   /**
00253    */
00254   virtual Real norm(const LevelData<EBCellFAB>& a_rhs,
00255                     int                         a_ord);
00256 
00257   ///
00258   /**
00259    */
00260   virtual void setToZero(LevelData<EBCellFAB>& a_lhs);
00261 
00262   ///
00263   /**
00264    */
00265   virtual void setVal(LevelData<EBCellFAB>& a_lhs, const Real& a_value);
00266 
00267   ///
00268   /**
00269    */
00270   virtual void createCoarser(LevelData<EBCellFAB>&       a_coarse,
00271                              const LevelData<EBCellFAB>& a_fine,
00272                              bool                        a_ghosted);
00273 
00274   ///
00275   /**
00276    */
00277   virtual void relax(LevelData<EBCellFAB>&       a_e,
00278                      const LevelData<EBCellFAB>& a_residual,
00279                      int                         a_iterations);
00280 
00281 
00282   virtual void relaxGauSai(LevelData<EBCellFAB>&       a_e,
00283                           const LevelData<EBCellFAB>& a_residual,
00284                           int                         a_iterations);
00285 
00286   virtual void relaxPoiJac(LevelData<EBCellFAB>&       a_e,
00287                            const LevelData<EBCellFAB>& a_residual,
00288                            int                         a_iterations);
00289 
00290   ///
00291   /**
00292      Calculate restricted residual:
00293      a_resCoarse[2h] = I[h->2h] (a_rhsFine[h] - L[h](a_phiFine[h]))
00294   */
00295   virtual void restrictResidual(LevelData<EBCellFAB>&       a_resCoarse,
00296                                 LevelData<EBCellFAB>&       a_phiFine,
00297                                 const LevelData<EBCellFAB>& a_rhsFine);
00298 
00299   ///
00300   /**
00301      Correct the fine solution based on coarse correction:
00302      a_phiThisLevel += I[2h->h] (a_correctCoarse)
00303   */
00304   virtual void prolongIncrement(LevelData<EBCellFAB>&       a_phiThisLevel,
00305                                 const LevelData<EBCellFAB>& a_correctCoarse);
00306 
00307   ///
00308   /** Refinement ratio between this level and coarser level.
00309       Returns 1 when there are no coarser AMRLevelOp objects */
00310   virtual int refToCoarser();
00311 
00312   ///
00313   /** Refinement ratio between this level and coarser level.
00314       Returns 1 when there are no coarser AMRLevelOp objects */
00315   virtual int refToFiner();
00316 
00317   ///
00318   /** a_residual = a_rhs - L(a_phi, a_phiFine, a_phiCoarse) */
00319   virtual void AMRResidual(LevelData<EBCellFAB>& a_residual,
00320                            const LevelData<EBCellFAB>& a_phiFine,
00321                            const LevelData<EBCellFAB>& a_phi,
00322                            const LevelData<EBCellFAB>& a_phiCoarse,
00323                            const LevelData<EBCellFAB>& a_rhs,
00324                            bool a_homogeneousBC,
00325                            AMRLevelOp<LevelData<EBCellFAB> >* a_finerOp);
00326 
00327   ///
00328   /** a_residual = a_rhs - L(a_phi, a_phiCoarse)  */
00329   virtual void AMRResidualNF(LevelData<EBCellFAB>& a_residual,
00330                              const LevelData<EBCellFAB>& a_phi,
00331                              const LevelData<EBCellFAB>& a_phiCoarse,
00332                              const LevelData<EBCellFAB>& a_rhs,
00333                              bool a_homogeneousBC);
00334 
00335 
00336   ///
00337   /** a_residual = a_rhs - L(a_phi, a_phiFine, a_phiCoarse) */
00338   virtual void AMROperator(LevelData<EBCellFAB>& a_LofPhi,
00339                            const LevelData<EBCellFAB>& a_phiFine,
00340                            const LevelData<EBCellFAB>& a_phi,
00341                            const LevelData<EBCellFAB>& a_phiCoarse,
00342                            bool a_homogeneousBC,
00343                            AMRLevelOp<LevelData<EBCellFAB> >* a_finerOp);
00344 
00345   ///
00346   /** a_residual = a_rhs - L(a_phi, a_phiCoarse)  */
00347   virtual void AMROperatorNF(LevelData<EBCellFAB>& a_LofPhi,
00348                              const LevelData<EBCellFAB>& a_phi,
00349                              const LevelData<EBCellFAB>& a_phiCoarse,
00350                              bool a_homogeneousBC);
00351 
00352   ///
00353   /** a_resCoarse = I[h-2h] (a_residual - L(a_correction, a_coarseCorrection)) */
00354   virtual void AMRRestrict(LevelData<EBCellFAB>& a_resCoarse,
00355                            const LevelData<EBCellFAB>& a_residual,
00356                            const LevelData<EBCellFAB>& a_correction,
00357                            const LevelData<EBCellFAB>& a_coarseCorrection, 
00358                            bool a_skip_res = false );
00359 
00360   ///
00361   /** a_correction += I[2h->h](a_coarseCorrection) */
00362   virtual void AMRProlong(LevelData<EBCellFAB>&       a_correction,
00363                           const LevelData<EBCellFAB>& a_coarseCorrection);
00364 
00365   ///
00366   /** a_residual = a_residual - L(a_correction, a_coarseCorrection) */
00367   virtual void AMRUpdateResidual(LevelData<EBCellFAB>&       a_residual,
00368                                  const LevelData<EBCellFAB>& a_correction,
00369                                  const LevelData<EBCellFAB>& a_coarseCorrection);
00370 
00371   void reflux(LevelData<EBCellFAB>& a_residual,
00372               const LevelData<EBCellFAB>& a_phiFine,
00373               const LevelData<EBCellFAB>& a_phi,
00374               AMRLevelOp<LevelData<EBCellFAB> >* a_finerOp);
00375 
00376   void gsrbColor(LevelData<EBCellFAB>&       a_phi,
00377                  const LevelData<EBCellFAB>& a_lph,
00378                  const LevelData<EBCellFAB>& a_rhs,
00379                  const IntVect&              a_color);
00380 
00381   void getDivFStencil(VoFStencil&      a_vofStencil,
00382                       const VolIndex&  a_vof,
00383                       const DataIndex& a_dit);
00384 
00385   void getFluxStencil(VoFStencil&      a_fluxStencil,
00386                       const FaceIndex& a_face,
00387                       const DataIndex& a_dit);
00388 
00389   void getFaceCenteredFluxStencil(VoFStencil&      a_fluxStencil,
00390                                   const FaceIndex& a_face,
00391                                   const DataIndex& a_dit);
00392 
00393   void incrOpRegularDir(EBCellFAB&             a_lhs,
00394                         const EBCellFAB&       a_phi,
00395                         const bool&            a_homogeneous,
00396                         const int&             a_dir,
00397                         const DataIndex&       a_datInd);
00398   void applyOpIrregular(EBCellFAB&             a_lhs,
00399                         const EBCellFAB&       a_phi,
00400                         const bool&            a_homogeneous,
00401                         const DataIndex&       a_datInd);
00402 protected:
00403 
00404   static bool                     s_turnOffBCs;
00405   void defineStencils();
00406   int                             m_relaxType;
00407   const IntVect                   m_ghostCellsPhi;
00408   const IntVect                   m_ghostCellsRHS;
00409 
00410   RefCountedPtr<EBQuadCFInterp>   m_quadCFIWithCoar;
00411 
00412   EBLevelGrid                     m_eblg;
00413   EBLevelGrid                     m_eblgFine;
00414   EBLevelGrid                     m_eblgCoar;
00415   EBLevelGrid                     m_eblgCoarMG;
00416   EBLevelGrid                     m_eblgCoarsenedFine;
00417 
00418   RefCountedPtr<ConductivityBaseDomainBC>     m_domainBC;
00419   RefCountedPtr<ConductivityBaseEBBC>         m_ebBC;
00420 
00421   Real                            m_dxFine;
00422   Real                            m_dx;
00423   Real                            m_dxCoar;
00424   RefCountedPtr<LevelData<EBCellFAB> >          m_acoef;
00425   RefCountedPtr<LevelData<EBFluxFAB> >          m_bcoef;
00426   RefCountedPtr<LevelData<BaseIVFAB<Real> > >   m_bcoIrreg;
00427 
00428   Real                            m_alpha;
00429   Real                            m_beta;
00430   //weights that get multiplied by alpha
00431   LayoutData<BaseIVFAB<Real> >       m_alphaDiagWeight;
00432   //weights that get multiplied by beta
00433   LayoutData<BaseIVFAB<Real> >       m_betaDiagWeight;
00434   int                             m_refToFine;
00435   int                             m_refToCoar;
00436   bool                            m_hasFine;
00437   bool                            m_hasInterpAve;
00438   bool                            m_hasCoar;
00439 
00440   //restriction object
00441   EBMGAverage                    m_ebAverage;
00442   //prolongation object
00443   EBMGInterp                     m_ebInterp;
00444 
00445   //stencils for operator evaluation
00446   //LayoutData<RefCountedPtr<EBStencil> >  m_opEBStencil;
00447   LayoutData<BaseIVFAB<VoFStencil> >  m_opStencil;
00448   //stencils for operator evaluation on gauss-seidel colors
00449   LevelData<EBCellFAB>       m_relCoef;
00450 
00451   //cache the vofiterators
00452   //for irregular cell iteration (includes buffer around multivalued cells)
00453   LayoutData<VoFIterator >                     m_vofIterIrreg;
00454   LayoutData<VoFIterator >                     m_vofIterMulti;
00455   //for domain boundary conditions at ir regular cells
00456   LayoutData<VoFIterator >                     m_vofIterDomLo[CH_SPACEDIM];
00457   LayoutData<VoFIterator >                     m_vofIterDomHi[CH_SPACEDIM];
00458 
00459   void calculateRelaxationCoefficient();
00460 
00461   // Coarse-fine stencils for homogeneous CFInterp
00462   LayoutData<CFIVS> m_loCFIVS[SpaceDim];
00463   LayoutData<CFIVS> m_hiCFIVS[SpaceDim];
00464 
00465   //flux register with finer level
00466   EBFastFR       m_fastFR;
00467 
00468   //special mg objects for when we do not have
00469   //a coarser level or when the refinement to coarse
00470   //is greater than two
00471   //flag for when we need special MG objects
00472   bool                        m_hasMGObjects;
00473   bool                        m_layoutChanged;
00474   //stuff below is only defined if m_hasMGObjects==true
00475   EBMGAverage                 m_ebAverageMG;
00476   EBMGInterp                  m_ebInterpMG;
00477   DisjointBoxLayout           m_dblCoarMG;
00478   EBISLayout                  m_ebislCoarMG;
00479   ProblemDomain               m_domainCoarMG;
00480 
00481   Vector<IntVect> m_colors;
00482 
00483 private:
00484 
00485   void incrementFRCoar(EBFastFR&                   a_fluxReg,
00486                        const LevelData<EBCellFAB>& a_phiFine,
00487                        const LevelData<EBCellFAB>& a_phi);
00488 
00489   void incrementFRFine(EBFastFR&                   a_fluxReg,
00490                        const LevelData<EBCellFAB>& a_phiFine,
00491                        const LevelData<EBCellFAB>& a_phi,
00492                        AMRLevelOp<LevelData<EBCellFAB> >* a_finerOp);
00493 
00494   void getFlux(FArrayBox&                    a_flux,
00495                const FArrayBox&              a_phi,
00496                const Box&                    a_faceBox,
00497                const int&                    a_idir,
00498                const Real&                   a_dx,
00499                const DataIndex&              a_datInd);
00500 
00501 
00502 
00503   void applyCFBCs(LevelData<EBCellFAB>&             a_phi,
00504                   const LevelData<EBCellFAB>* const a_phiCoarse,
00505                   bool a_homogeneousCFBC);
00506 
00507   void getOpVoFStencil(VoFStencil&     a_stencil,
00508                        const EBISBox&  a_ebisbox,
00509                        const VolIndex& a_vof);
00510 
00511   void getOpVoFStencil(VoFStencil&             a_stencil,
00512                        const int&              a_dir,
00513                        const Vector<VolIndex>& a_allMonotoneVoFs,
00514                        const EBISBox&          a_ebisbox,
00515                        const VolIndex&         a_vof,
00516                        const bool&             a_lowOrder);
00517 
00518 
00519   void getOpFaceStencil(VoFStencil&             a_stencil,
00520                         const Vector<VolIndex>& a_allMonotoneVofs,
00521                         const EBISBox&          a_ebisbox,
00522                         const VolIndex&         a_vof,
00523                         int                     a_dir,
00524                         const Side::LoHiSide&   a_side,
00525                         const FaceIndex&        a_face,
00526                         const bool&             a_lowOrder);
00527 
00528   void levelJacobi(LevelData<EBCellFAB>&       a_phi,
00529                    const LevelData<EBCellFAB>& a_rhs);
00530 
00531   void applyHomogeneousCFBCs(LevelData<EBCellFAB>&   a_phi);
00532 
00533   void applyHomogeneousCFBCs(EBCellFAB&            a_phi,
00534                              const DataIndex&      a_datInd,
00535                              int                   a_idir,
00536                              Side::LoHiSide        a_hiorlo);
00537 private:
00538   //copy constructor and operator= disallowed for all the usual reasons
00539   slowEBCO(const slowEBCO& a_opin)
00540   {
00541     MayDay::Error("invalid operator");
00542   }
00543 
00544   void operator=(const slowEBCO& a_opin)
00545   {
00546     MayDay::Error("invalid operator");
00547   }
00548 };
00549 
00550 
00551 #include "NamespaceFooter.H"
00552 #endif