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