Chombo + EB: AMRPoissonOp.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 _AMRPOISSONOP_H_
00012 #define _AMRPOISSONOP_H_
00013 
00014 #include "REAL.H"
00015 #include "Box.H"
00016 #include "FArrayBox.H"
00017 #include "FluxBox.H"
00018 #include "QuadCFInterp.H"
00019 #include "LevelFluxRegister.H"
00020 #include "LevelDataOps.H"
00021 #include "AMRMultiGrid.H"
00022 #include "BCFunc.H"
00023 #include "BaseLevelTGA.H"
00024 #include "CFRegion.H"
00025 #include "AMRIO.H"
00026 #include "CornerCopier.H"
00027 
00028 #include "NamespaceHeader.H"
00029 
00030 ///
00031 /**
00032    Operator for solving (alpha I + beta*Laplacian)(phi) = rho
00033    over an AMR hierarchy.
00034 */
00035 class AMRPoissonOp : public LevelTGAHelmOp<LevelData<FArrayBox> , FluxBox>
00036 {
00037 public:
00038 
00039   /**
00040      \name AMRPoissonOp functions */
00041   /*@{*/
00042 
00043   ///
00044   /**
00045    */
00046   AMRPoissonOp()
00047   {
00048   }
00049 
00050   ///
00051   /**
00052    */
00053   virtual ~AMRPoissonOp()
00054   {
00055   }
00056 
00057   ///
00058   /** full define function for AMRLevelOp with both coarser and finer levels */
00059   void define(const DisjointBoxLayout&   a_grids,
00060               const DisjointBoxLayout&   a_gridsFiner,
00061               const DisjointBoxLayout&   a_gridsCoarser,
00062               const Real&                a_dxLevel,
00063               int                        a_refRatio,
00064               int                        a_refRatioFiner,
00065               const ProblemDomain&       a_domain,
00066               BCHolder                   a_bc,
00067               const Copier&              a_exchange,
00068               const CFRegion&            a_cfregion,
00069               const int                  a_nComp = 1);
00070 
00071   /** full define function for AMRLevelOp with finer levels, but no coarser */
00072   void define(const DisjointBoxLayout&   a_grids,
00073               const DisjointBoxLayout&   a_gridsFiner,
00074               const Real&                a_dxLevel,
00075               int                        a_refRatio, // dummy arg, send in 1
00076               int                        a_refRatioFiner,
00077               const ProblemDomain&       a_domain,
00078               BCHolder                   a_bc,
00079               const Copier&              a_exchange,
00080               const CFRegion&            a_cfregion,
00081               const int                  a_nComp = 1);
00082 
00083   ///
00084   /**
00085      define function for AMRLevelOp which has no finer AMR level
00086      Jamie Parkinson added a_numComp parameter so that we define CF interp objects correctly
00087   */
00088   void define(const DisjointBoxLayout&   a_grids,
00089               const DisjointBoxLayout&   a_baseBAPtr,
00090               const Real&                a_dxLevel,
00091               int                        a_refRatio,
00092               const ProblemDomain&       a_domain,
00093               BCHolder                   a_bc,
00094               const Copier&              a_exchange,
00095               const CFRegion&            a_cfregion,
00096               int                        a_numComp=1);
00097 
00098   ///
00099   /**
00100      define function for AMRLevelOp which has no finer or coarser AMR level
00101    */
00102   void define(const DisjointBoxLayout&   a_grids,
00103               const Real&                a_dx,
00104               const ProblemDomain&       a_domain,
00105               BCHolder                   a_bc,
00106               const Copier&              a_exchange,
00107               const CFRegion&            a_cfregion);
00108 
00109   ///
00110   /**
00111      define function for AMRLevelOp which has no finer or coarser AMR level
00112    */
00113   void define(const DisjointBoxLayout&   a_grids,
00114               const Real&                a_dx,
00115               const ProblemDomain&       a_domain,
00116               BCHolder                   a_bc);
00117 
00118   /// Full define function that mimics the old PoissonOp.
00119   /**
00120       Makes all coarse-fine information and sets internal variables
00121   */
00122   void define(const DisjointBoxLayout& a_grids,
00123               const DisjointBoxLayout* a_baseBAPtr,
00124               Real                     a_dxLevel,
00125               int                      a_refRatio,
00126               const ProblemDomain&     a_domain,
00127               BCHolder                 a_bc);
00128 
00129   virtual void residual(LevelData<FArrayBox>&       a_lhs,
00130                         const LevelData<FArrayBox>& a_phi,
00131                         const LevelData<FArrayBox>& a_rhs,
00132                         bool                        a_homogeneous = false);
00133 
00134   /// Residual which uses coarse-level boundary conditions, but ignores finer levels
00135   /** useful for full-multigrid types of schemes
00136    */
00137   virtual void residualNF(LevelData<FArrayBox>& a_lhs,
00138                           LevelData<FArrayBox>& a_phi,
00139                           const LevelData<FArrayBox>* a_phiCoarse,
00140                           const LevelData<FArrayBox>& a_rhs,
00141                           bool a_homogeneous = false);
00142 
00143  /// despite what you might think, the "I" here means "Ignore the coarse-fine boundary"
00144   virtual void residualI(LevelData<FArrayBox>&       a_lhs,
00145                          const LevelData<FArrayBox>& a_phi,
00146                          const LevelData<FArrayBox>& a_rhs,
00147                          bool                        a_homogeneous = false);
00148 
00149   virtual void preCond(LevelData<FArrayBox>&       a_correction,
00150                        const LevelData<FArrayBox>& a_residual);
00151 
00152   virtual void applyOp(LevelData<FArrayBox>&       a_lhs,
00153                        const LevelData<FArrayBox>& a_phi,
00154                        bool                        a_homogeneous = false);
00155 
00156   /// despite what you might think, the "I" here means "Ignore the coarse-fine boundary"
00157   virtual void applyOpI(LevelData<FArrayBox>&       a_lhs,
00158                         const LevelData<FArrayBox>& a_phi,
00159                         bool                        a_homogeneous = false);
00160 
00161   virtual void applyOpNoBoundary(LevelData<FArrayBox>& a_lhs,
00162                                  const LevelData<FArrayBox>& a_rhs) ;
00163 
00164   virtual void create(LevelData<FArrayBox>&       a_lhs,
00165                       const LevelData<FArrayBox>& a_rhs);
00166 
00167   virtual void createCoarsened(LevelData<FArrayBox>&       a_lhs,
00168                                const LevelData<FArrayBox>& a_rhs,
00169                                const int&                  a_refRat);
00170 
00171   virtual void assign(LevelData<FArrayBox>&       a_lhs,
00172                       const LevelData<FArrayBox>& a_rhs);
00173 
00174   virtual void assignLocal(LevelData<FArrayBox>&       a_lhs,
00175                            const LevelData<FArrayBox>& a_rhs);
00176 
00177   virtual void buildCopier(Copier& a_copier,
00178                            const LevelData<FArrayBox>& a_lhs,
00179                            const LevelData<FArrayBox>& a_rhs);
00180 
00181   virtual void assignCopier(LevelData<FArrayBox>&       a_lhs,
00182                             const LevelData<FArrayBox>& a_rhs,
00183                             const Copier&               a_copier);
00184 
00185   virtual void zeroCovered(LevelData<FArrayBox>& a_lhs,
00186                            LevelData<FArrayBox>& a_rhs,
00187                            const Copier&         a_copier);
00188 
00189   virtual Real dotProduct(const LevelData<FArrayBox>& a_1,
00190                           const LevelData<FArrayBox>& a_2);
00191   /* multiple dot products (for GMRES) */
00192   virtual void mDotProduct(const LevelData<FArrayBox>& a_1,
00193                            const int a_sz,
00194                            const LevelData<FArrayBox> a_2[],
00195                            Real a_mdots[]);
00196 
00197   virtual void incr(LevelData<FArrayBox>&       a_lhs,
00198                     const LevelData<FArrayBox>& a_x,
00199                     Real                        a_scale);
00200 
00201   virtual void axby(LevelData<FArrayBox>&       a_lhs,
00202                     const LevelData<FArrayBox>& a_x,
00203                     const LevelData<FArrayBox>& a_y,
00204                     Real                        a_a,
00205                     Real                        a_b);
00206 
00207   virtual void scale(LevelData<FArrayBox>& a_lhs,
00208                      const Real&           a_scale);
00209 
00210   virtual Real norm(const LevelData<FArrayBox>& a_x,
00211                     int                         a_ord);
00212 
00213   virtual Real localMaxNorm(const LevelData<FArrayBox>& a_x);
00214 
00215   virtual void setToZero( LevelData<FArrayBox>& a_x);
00216   /*@}*/
00217 
00218   /**
00219      \name MGLevelOp functions */
00220   /*@{*/
00221 
00222   virtual void relax(LevelData<FArrayBox>&       a_e,
00223                      const LevelData<FArrayBox>& a_residual,
00224                      int                         a_iterations);
00225 
00226   /// Do CF interpolation then relax as normal
00227   /** useful for full-multigrid/FAS schemes
00228    */
00229   virtual void relaxNF(LevelData<FArrayBox>& a_phi,
00230                        const LevelData<FArrayBox>* a_phiCoarse,
00231                        const LevelData<FArrayBox>& a_rhs,
00232                        int a_iterations);
00233 
00234   virtual void createCoarser(LevelData<FArrayBox>&       a_coarse,
00235                              const LevelData<FArrayBox>& a_fine,
00236                              bool                        a_ghosted);
00237   /**
00238      calculate restricted residual
00239      a_resCoarse[2h] = I[h->2h] (rhsFine[h] - L[h](phiFine[h])
00240   */
00241   virtual void restrictResidual(LevelData<FArrayBox>&       a_resCoarse,
00242                                 LevelData<FArrayBox>&       a_phiFine,
00243                                 const LevelData<FArrayBox>& a_rhsFine);
00244 
00245   /**
00246      correct the fine solution based on coarse correction
00247      a_phiThisLevel += I[2h->h](a_correctCoarse)
00248   */
00249   virtual void prolongIncrement(LevelData<FArrayBox>&       a_phiThisLevel,
00250                                 const LevelData<FArrayBox>& a_correctCoarse);
00251 
00252   /*@}*/
00253 
00254   /**
00255      \name AMRLevelOp functions */
00256   /*@{*/
00257 
00258   /** returns 1 when there are no coarser AMRLevelOp objects */
00259   virtual int refToCoarser()
00260   {
00261     return m_refToCoarser;
00262   }
00263 
00264   /** a_residual = a_rhs - L(a_phi, a_phiFine, a_phiCoarse) */
00265   virtual void AMRResidual(LevelData<FArrayBox>&              a_residual,
00266                            const LevelData<FArrayBox>&        a_phiFine,
00267                            const LevelData<FArrayBox>&        a_phi,
00268                            const LevelData<FArrayBox>&        a_phiCoarse,
00269                            const LevelData<FArrayBox>&        a_rhs,
00270                            bool                               a_homogeneousPhysBC,
00271                            AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00272 
00273   /** residual assuming no more coarser AMR levels */
00274 
00275   virtual void AMRResidualNC(LevelData<FArrayBox>&              a_residual,
00276                              const LevelData<FArrayBox>&        a_phiFine,
00277                              const LevelData<FArrayBox>&        a_phi,
00278                              const LevelData<FArrayBox>&        a_rhs,
00279                              bool                               a_homogeneousPhysBC,
00280                              AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00281 
00282   /** a_residual = a_rhs - L(a_phi, a_phiCoarse)  */
00283   virtual void AMRResidualNF(LevelData<FArrayBox>&       a_residual,
00284                              const LevelData<FArrayBox>& a_phi,
00285                              const LevelData<FArrayBox>& a_phiCoarse,
00286                              const LevelData<FArrayBox>& a_rhs,
00287                              bool                        a_homogeneousPhysBC);
00288 
00289   ///
00290   /**
00291      Apply the AMR operator, including coarse-fine matching
00292   */
00293   virtual void AMROperator(LevelData<FArrayBox>&              a_LofPhi,
00294                            const LevelData<FArrayBox>&        a_phiFine,
00295                            const LevelData<FArrayBox>&        a_phi,
00296                            const LevelData<FArrayBox>&        a_phiCoarse,
00297                            bool                               a_homogeneousDomBC,
00298                            AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00299 
00300   ///
00301   /**
00302       Apply the AMR operator, including coarse-fine matching
00303       assume no coarser AMR level
00304   */
00305   virtual void AMROperatorNC(LevelData<FArrayBox>&              a_LofPhi,
00306                              const LevelData<FArrayBox>&        a_phiFine,
00307                              const LevelData<FArrayBox>&        a_phi,
00308                              bool                               a_homogeneousBC,
00309                              AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00310 
00311   ///
00312   /**
00313       Apply the AMR operator, including coarse-fine matching.
00314       assume no finer AMR level
00315   */
00316   virtual void AMROperatorNF(LevelData<FArrayBox>&       a_LofPhi,
00317                              const LevelData<FArrayBox>& a_phi,
00318                              const LevelData<FArrayBox>& a_phiCoarse,
00319                              bool                        a_homogeneousBC);
00320 
00321   /**
00322       a_resCoarse = I[h-2h]( a_residual - L(a_correction, a_coarseCorrection))
00323       it is assumed that a_resCoarse has already been filled in with the coarse
00324       version of AMRResidualNF and that this operation is free to overwrite
00325       in the overlap regions.
00326   */
00327   virtual void AMRRestrict(LevelData<FArrayBox>&       a_resCoarse,
00328                            const LevelData<FArrayBox>& a_residual,
00329                            const LevelData<FArrayBox>& a_correction,
00330                            const LevelData<FArrayBox>& a_coarseCorrection,
00331                            bool a_skip_res = false );
00332 
00333   virtual void AMRRestrictS(LevelData<FArrayBox>&       a_resCoarse,
00334                             const LevelData<FArrayBox>& a_residual,
00335                             const LevelData<FArrayBox>& a_correction,
00336                             const LevelData<FArrayBox>& a_coarseCorrection,
00337                             LevelData<FArrayBox>&       a_scratch,
00338                             bool a_skip_res = false );
00339 
00340   /**
00341       a_correction += I[h->h](a_coarseCorrection)
00342   */
00343   virtual void AMRProlong(LevelData<FArrayBox>&       a_correction,
00344                           const LevelData<FArrayBox>& a_coarseCorrection);
00345 
00346   /**
00347       optimization of AMRProlong that sends in the existing temporary and copier
00348   */
00349   virtual void AMRProlongS(LevelData<FArrayBox>&       a_correction,
00350                            const LevelData<FArrayBox>& a_coarseCorrection,
00351                            LevelData<FArrayBox>&       a_temp,
00352                            const Copier&               a_copier);
00353   /**
00354       optimization of AMRProlong that sends in the existing temporary and copier -- higher order
00355   */
00356   virtual void AMRProlongS_2(LevelData<FArrayBox>&       a_correction,
00357                              const LevelData<FArrayBox>& a_coarseCorrection,
00358                              LevelData<FArrayBox>&       a_temp,
00359                              const Copier&               a_copier,
00360                              const Copier&         a_cornerCopier,
00361                              const AMRLevelOp<LevelData<FArrayBox> >*  a_crsOp
00362                              );
00363 
00364   /**
00365       a_residual = a_residual - L(a_correction, a_coarseCorrection)
00366   */
00367   virtual void AMRUpdateResidual(LevelData<FArrayBox>&       a_residual,
00368                                  const LevelData<FArrayBox>& a_correction,
00369                                  const LevelData<FArrayBox>& a_coarseCorrection);
00370 
00371   ///
00372   /**
00373       compute norm over all cells on coarse not covered by finer
00374   */
00375   virtual Real AMRNorm(const LevelData<FArrayBox>& a_coarseResid,
00376                        const LevelData<FArrayBox>& a_fineResid,
00377                        const int&                  a_refRat,
00378                        const int&                  a_ord);
00379 
00380   /// For tga to reset stuff
00381   /**
00382   */
00383   virtual void setAlphaAndBeta(const Real& a_alpha,
00384                                const Real& a_beta);
00385 
00386   /// Change boundary conditions
00387   virtual void setBC(const BCHolder& a_bc);
00388 
00389   /// For tga stuff---in this case a noop
00390   virtual void        diagonalScale(LevelData<FArrayBox>& a_rhs,
00391                                     bool a_kappaWeighted)
00392   {
00393   }
00394 
00395   /// For tga stuff---in this case a noop
00396   virtual void divideByIdentityCoef(LevelData<FArrayBox>& a_rhs)
00397   {
00398   }
00399 
00400   virtual void fillGrad(const LevelData<FArrayBox>& a_phi)
00401   {
00402     ;//does not apply here
00403   }
00404 
00405   virtual void reflux(const LevelData<FArrayBox>&        a_phiFine,
00406                       const LevelData<FArrayBox>&        a_phi,
00407                       LevelData<FArrayBox>&              a_residual,
00408                       AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00409 
00410   virtual void getFlux(FluxBox&                    a_flux,
00411                        const LevelData<FArrayBox>& a_data,
00412                        const Box&                  a_grid,
00413                        const DataIndex&            a_dit,
00414                        Real                        a_scale)
00415   {
00416     const FArrayBox& data = a_data[a_dit];
00417     a_flux.define(a_grid, a_data.nComp());
00418     for (int idir = 0; idir < SpaceDim; idir++)
00419       {
00420         getFlux(a_flux[idir], data, a_flux[idir].box(), idir, 1);
00421         a_flux[idir] *= a_scale;
00422       }
00423   }
00424 
00425   virtual void write(const LevelData<FArrayBox>* a_data,
00426                      const char*                 a_filename);
00427 
00428   /*@}*/
00429 
00430   /// public constants
00431   Real m_alpha, m_beta, m_aCoef, m_bCoef;
00432 
00433   // needed for homogeneous interpolation
00434   // set by the factory
00435   Real  m_dxCrse;
00436 
00437   Vector<IntVect> m_colors;
00438   static int s_exchangeMode;
00439   static int s_relaxMode;
00440   static int s_maxCoarse;
00441   static int s_prolongType;
00442 
00443   virtual Real dx() const
00444   {
00445     return m_dx;
00446   }
00447 
00448 protected:
00449   Real                    m_dx;
00450   ProblemDomain           m_domain;
00451 
00452   LevelDataOps<FArrayBox> m_levelOps;
00453 
00454   BCHolder                m_bc;
00455 
00456   CFRegion                m_cfregion;
00457   Copier                  m_exchangeCopier;
00458   QuadCFInterp            m_interpWithCoarser;
00459 
00460   LevelFluxRegister       m_levfluxreg;
00461 
00462   DisjointBoxLayout       m_coarsenedMGrids;
00463 
00464   int                     m_refToCoarser;
00465   int                     m_refToFiner;
00466 
00467   virtual void levelGSRB(LevelData<FArrayBox>&       a_phi,
00468                          const LevelData<FArrayBox>& a_rhs);
00469 
00470   virtual void levelMultiColor(LevelData<FArrayBox>&       a_phi,
00471                                const LevelData<FArrayBox>& a_rhs);
00472 
00473   virtual void looseGSRB(LevelData<FArrayBox>&       a_phi,
00474                          const LevelData<FArrayBox>& a_rhs);
00475 
00476   virtual void overlapGSRB(LevelData<FArrayBox>&       a_phi,
00477                            const LevelData<FArrayBox>& a_rhs);
00478 
00479   virtual void levelGSRBLazy(LevelData<FArrayBox>&       a_phi,
00480                              const LevelData<FArrayBox>& a_rhs);
00481 
00482   virtual void levelJacobi(LevelData<FArrayBox>&       a_phi,
00483                            const LevelData<FArrayBox>& a_rhs);
00484 
00485   virtual void homogeneousCFInterp(LevelData<FArrayBox>& a_phif);
00486 
00487   virtual void homogeneousCFInterp(LevelData<FArrayBox>& a_phif,
00488                                    const DataIndex&      a_datInd,
00489                                    int                   a_idir,
00490                                    Side::LoHiSide        a_hiorlo);
00491 
00492   virtual void singleBoxCFInterp(FArrayBox& a_phi);
00493 
00494   virtual void interpOnIVSHomo(LevelData<FArrayBox>& a_phif,
00495                                const DataIndex&      a_datInd,
00496                                const int             a_idir,
00497                                const Side::LoHiSide  a_hiorlo,
00498                                const IntVectSet&     a_interpIVS);
00499 
00500   virtual void getFlux(FArrayBox&       a_flux,
00501                        const FArrayBox& a_data,
00502                        const Box&       a_edgebox,
00503                        int              a_dir,
00504                        int              a_ref = 1) const ;
00505 
00506   virtual void getFlux(FArrayBox&       a_flux,
00507                        const FArrayBox& a_data,
00508                        int              a_dir,
00509                        int              a_ref = 1) const ;
00510 };
00511 
00512 ///
00513 /**
00514    Factory to create AMRPoissonOps
00515  */
00516 class AMRPoissonOpFactory: public AMRLevelOpFactory<LevelData<FArrayBox> >
00517 {
00518 public:
00519   virtual ~AMRPoissonOpFactory()
00520   {
00521   }
00522 
00523   ///
00524   /**
00525      a_coarseDomain is the domain at the coarsest level.
00526      a_grids is the AMR  hierarchy.
00527      a_refRatios are the refinement ratios between levels.  The ratio lives
00528          with the coarser level so a_refRatios[ilev] is the ratio between
00529          ilev and ilev+1
00530      a_coarseDx is the grid spacing at the coarsest level.
00531      a_bc holds the boundary conditions.
00532      a_alpha is the identity coefficient
00533      a_beta is the laplacian coefficient.
00534   */
00535   void define(const ProblemDomain&             a_coarseDomain,
00536               const Vector<DisjointBoxLayout>& a_grids,
00537               const Vector<int>&               a_refRatios,
00538               const Real&                      a_coarsedx,
00539               BCHolder                         a_bc,
00540               Real                             a_alpha = 0.0,
00541               Real                             a_beta  = 1.0);
00542 
00543   // regular multigrid definition function --deprecated
00544   void define(const ProblemDomain&     a_domain,
00545               const DisjointBoxLayout& a_grid,
00546               const Real&              a_dx,
00547               BCHolder                 a_bc,
00548               int                      a_maxDepth = -1,
00549               Real                     a_alpha = 0.0,
00550               Real                     a_beta  = 1.0);
00551 
00552   ///
00553   virtual MGLevelOp<LevelData<FArrayBox> >* MGnewOp(const ProblemDomain& a_FineindexSpace,
00554                                                      int                  a_depth,
00555                                                      bool                 a_homoOnly = true);
00556 
00557   ///
00558   virtual  AMRLevelOp<LevelData<FArrayBox> >* AMRnewOp(const ProblemDomain& a_indexSpace);
00559 
00560   ///
00561   virtual int refToFiner(const ProblemDomain& a_domain) const;
00562 
00563 private:
00564   Vector<ProblemDomain>     m_domains;
00565   Vector<DisjointBoxLayout> m_boxes;
00566 
00567   Vector<Real> m_dx;
00568   Vector<int>  m_refRatios; // refinement to next coarser level
00569 
00570   BCHolder m_bc;
00571 
00572   Real m_alpha;
00573   Real m_beta;
00574 
00575   Vector<Copier>   m_exchangeCopiers;
00576   Vector<CFRegion> m_cfregion;
00577 };
00578 
00579 #include "NamespaceFooter.H"
00580 #endif