Chombo + EB: EBAdvectPatchIntegrator.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 _EBADVECTPATCHINTEGRATOR_H_
00012 #define _EBADVECTPATCHINTEGRATOR_H_
00013 
00014 #include "Box.H"
00015 #include "IntVectSet.H"
00016 #include "Vector.H"
00017 #include "CH_HDF5.H"
00018 #include "EBCellFAB.H"
00019 #include "EBFluxFAB.H"
00020 #include "EBISBox.H"
00021 #include "BaseIVFAB.H"
00022 #include "EBPhysIBC.H"
00023 #include "EBPhysIBCFactory.H"
00024 #include "Stencils.H"
00025 #include "AggStencil.H"
00026 #include "NamespaceHeader.H"
00027 
00028 #define EBAPI_TOL 1.e-12
00029 
00030 ///
00031 /**
00032    This does the same as EBPatchAdvect but presumably with less memory footprint. 
00033    and with a somewhat more sane interface.   
00034 
00035    I have pared down the EBPatchGodunov interface and tried to optimize for minimal 
00036    memory without completely destroying performance.
00037 */
00038 class EBAdvectPatchIntegrator 
00039 {
00040 public:
00041 
00042   ///boundary condions are set via setEBPhysIBC
00043   EBAdvectPatchIntegrator(const Box&            a_validBox,
00044                           const EBISBox&        a_ebisBox,
00045                           const IntVectSet&     a_coarseFineIVS,
00046                           const ProblemDomain&  a_domain,
00047                           const RealVect&       a_dx,
00048                           bool                  a_useSlopeLimiting);
00049   
00050   virtual ~EBAdvectPatchIntegrator()
00051   {;}
00052 
00053 
00054   void
00055   extrapolatePrim(EBFluxFAB&                       a_flux,
00056                   Vector<BaseIVFAB<Real> * >&      a_coveredFluxMinu,
00057                   Vector<BaseIVFAB<Real> * >&      a_coveredFluxPlus,
00058                   const Vector<IntVectSet >&       a_coveredSetsMinu,
00059                   const Vector<IntVectSet >&       a_coveredSetsPlus,
00060                   const Vector<Vector<VolIndex> >& a_coveredFaceMinu,
00061                   const Vector<Vector<VolIndex> >& a_coveredFacePlus,
00062                   const EBCellFAB&                 a_consState,
00063                   const EBCellFAB&                 a_source,
00064                   const DataIndex&                 a_dit,
00065                   const Real&                      a_time,
00066                   const Real&                      a_dt);
00067 
00068   ///
00069   /**
00070    */
00071   void
00072   advectiveDerivative(EBCellFAB&                       a_uDotDelRho,
00073                       const EBFluxFAB&                 a_faceRho,
00074                       const EBFluxFAB&                 a_faceVel,
00075                       const Vector<BaseIVFAB<Real>*> & a_coveredRhoLo,
00076                       const Vector<BaseIVFAB<Real>*> & a_coveredRhoHi,
00077                       const Vector<BaseIVFAB<Real>*> & a_coveredVelLo,
00078                       const Vector<BaseIVFAB<Real>*> & a_coveredVelHi,
00079                       const Vector<Vector<VolIndex> >& a_coveredFaceLo,
00080                       const Vector<Vector<VolIndex> >& a_coveredFaceHi,
00081                       const Box&                       a_box);
00082 
00083 
00084   ///
00085   /**
00086      Version that leaves out the covered face stuff.  This is wrong near the EB but
00087      the codes that use it overwrite the EB stuff anyway.
00088    */
00089   void
00090   advectiveDerivative(EBCellFAB&                       a_uDotDelRho,
00091                       const EBFluxFAB&                 a_faceRho,
00092                       const EBFluxFAB&                 a_faceVel,
00093                       const Box&                       a_box);
00094   ///
00095   /**
00096    */
00097   void
00098   setEBPhysIBC(const EBPhysIBCFactory& a_bc)
00099   {
00100     m_bc = RefCountedPtr<EBPhysIBC>(a_bc.create());
00101     m_bc->define(m_domain, m_dx);
00102     m_isBCSet = true;
00103   }
00104 
00105   void setVelocities(const EBCellFAB& a_normalVel,
00106                      const EBFluxFAB& a_advectionVel)
00107   {
00108     m_normalVelPtr    = &a_normalVel;
00109     m_advectionVelPtr = &a_advectionVel;
00110     m_isVelSet = true;
00111   }
00112   
00113 
00114   /**
00115      For when EBFlux is always zero.
00116    */
00117   void
00118   consUndividedDivergence(BaseIVFAB<Real>&       a_divF,
00119                           const BaseIFFAB<Real>  a_centroidFlux[SpaceDim],
00120                           const IntVectSet&      a_ivs);
00121 
00122 
00123   void
00124   consUndividedDivergence(BaseIVFAB<Real>&       a_divF,
00125                           const BaseIFFAB<Real>  a_centroidFlux[SpaceDim],
00126                           const BaseIVFAB<Real>& a_ebIrregFlux,
00127                           const IntVectSet&      a_ivs);
00128 
00129   void
00130   interpolateFluxToCentroids(BaseIFFAB<Real>              a_centroidFlux[SpaceDim],
00131                              const BaseIFFAB<Real>* const a_fluxInterpolant[SpaceDim],
00132                              const IntVectSet&            a_irregIVS);
00133 
00134   void
00135   averageVelToCC(EBCellFAB&                        a_normalVel,
00136                  const EBFluxFAB&                  a_advectionVel,
00137                  const Vector<BaseIVFAB<Real> * >& a_coveredVeloLo,
00138                  const Vector<BaseIVFAB<Real> * >& a_coveredVeloHi,
00139                  const Vector<Vector<VolIndex> >&  a_coveredFaceLo,
00140                  const Vector<Vector<VolIndex> >&  a_coveredFaceHi,
00141                  const Box&                        a_box) const;
00142 
00143 
00144   void
00145   extrapolateBCG(EBFluxFAB&                       a_flux,
00146                  const EBCellFAB&                 a_consState,
00147                  const EBCellFAB&                 a_source,
00148                  const DataIndex&                 a_dit,
00149                  const Real&                      a_time,
00150                  const Real&                      a_dt);
00151 
00152 
00153   void
00154   mvExtrapolateBCG(EBFluxFAB                                      &  a_flux,
00155                    const EBCellFAB                                &  a_consState,
00156                    const EBFluxFAB                                &  a_advectionVel,
00157                    const EBCellFAB                                &  a_normalVel,
00158                    const EBCellFAB                                &  a_source,
00159                    const Vector<RefCountedPtr<EBPhysIBCFactory> > &  a_allAdvectBC,
00160                    const DataIndex                                &  a_dit,
00161                    const Real                                     &  a_time,
00162                    const Real                                     &  a_dt,
00163                    const int                                      &  a_doingVel);
00164 
00165   ///here are a couple of awful hooks necessary to get the minutiae of the algorithm correct
00166   static void setCurComp(int a_curComp)
00167   {
00168     s_curComp = a_curComp;
00169   }
00170   ///
00171   static void setDoingVel(int a_yesorno)
00172   {
00173     s_doingVel = a_yesorno;
00174   }
00175 
00176   static int getDoingVel()
00177   {
00178     return s_doingVel;
00179   }
00180   static int getCurComp()
00181   {
00182     return s_curComp;
00183   }
00184 
00185   ///
00186   /**
00187      This is called by EBAMRNoSubCycle.   The insane version that uses cell-centered
00188      data holders for face data (the plus-minus stuff) is called internally.
00189   */
00190   void
00191   extrapToCoveredFaces(BaseIVFAB<Real>&        a_extendedPrim,
00192                        const EBFaceFAB&        a_primFace,
00193                        const EBCellFAB&        a_primState,
00194                        const Vector<VolIndex>& a_coveredFaces,
00195                        const int&              a_faceDir,
00196                        const Side::LoHiSide&   a_sd,
00197                        const Box&              a_box);
00198 
00199   void setMaxMin(const Real& a_maxVal,
00200                  const Real& a_minVal)
00201   {
00202     m_isMaxMinSet = true;
00203     m_maxVal = a_maxVal;
00204     m_minVal = a_minVal;
00205   }
00206 private:
00207   const EBFluxFAB* m_advectionVelPtr;
00208   const EBCellFAB* m_normalVelPtr;
00209   bool m_isVelSet;
00210   bool m_isMaxMinSet;
00211   bool m_useLimiting;
00212   Real m_maxVal;
00213   Real m_minVal;
00214 
00215   RefCountedPtr<EBPhysIBC> m_bc;
00216   bool             m_isBCSet;
00217   Box              m_validBox;
00218   ProblemDomain    m_domain;
00219   EBISBox          m_ebisBox;
00220   RealVect         m_dx;
00221   IntVectSet       m_cfivs;
00222   Box              m_validBoxG4;
00223   Vector<VolIndex> m_irregVoFs;
00224 
00225   IntVectSet       m_coveredSetsPlusG4[SpaceDim];
00226   IntVectSet       m_coveredSetsMinuG4[SpaceDim];
00227   Vector<VolIndex> m_coveredFacePlusG4[SpaceDim];
00228   Vector<VolIndex> m_coveredFaceMinuG4[SpaceDim];
00229 
00230   ///these exist because special things have to be done for velocity
00231   static int  s_curComp;
00232   static int  s_doingVel;
00233 
00234   /// weak construction is bad.
00235   EBAdvectPatchIntegrator();
00236 
00237   /// internal functions (all the madness below probably needs to get cleaned up)
00238 
00239   ///
00240   void extrapolatePrim3D(EBCellFAB           a_primMinu[SpaceDim],
00241                          EBCellFAB           a_primPlus[SpaceDim],
00242                          const EBCellFAB&    a_consState,
00243                          const EBCellFAB&    a_source,
00244                          const DataIndex&    a_dit,
00245                          const Real&         a_time,
00246                          const Real&         a_dt);
00247 
00248   ///
00249   void
00250   extrapolateBCG(EBCellFAB           a_primMinu[SpaceDim],
00251                  EBCellFAB           a_primPlus[SpaceDim],
00252                  const EBCellFAB&    a_consState,
00253                  const EBCellFAB&    a_source,
00254                  const DataIndex&    a_dit,
00255                  const Real&                      a_time,
00256                  const Real&                      a_dt);
00257 
00258   ///
00259   void extrapolatePrim2D(EBCellFAB           a_primMinu[SpaceDim],
00260                          EBCellFAB           a_primPlus[SpaceDim],
00261                          const EBCellFAB&    a_consState,
00262                          const EBCellFAB&    a_source,
00263                          const DataIndex&    a_dit,
00264                          const Real&         a_time,
00265                          const Real&         a_dt);
00266 
00267   ///
00268   void
00269   riemann(BaseIVFAB<Real>&        a_coveredPrim,
00270           const BaseIVFAB<Real>&  a_exteState,
00271           const EBCellFAB&        a_primState,
00272           const Vector<VolIndex>& a_vofset,
00273           const int&              a_faceDir,
00274           const Side::LoHiSide&   a_sd,
00275           const Box&       a_box);
00276 
00277   ///
00278   void
00279   riemann(EBFaceFAB&       a_primGdnv,
00280           const EBCellFAB& a_primLeft,
00281           const EBCellFAB& a_primRigh,
00282           const int&       a_faceDir,
00283           const Box&       a_box);
00284 
00285   ///
00286   FaceStencil getInterpStencil(const FaceIndex& a_face) const;
00287 
00288 
00289   void
00290   pointExtrapToCovered2D(Vector<Real>&           a_extrapVal,
00291                          const EBFaceFAB&        a_primFace,
00292                          const EBCellFAB&        a_primState,
00293                          const int&              a_faceDir,
00294                          const VolIndex&         a_vof,
00295                          const RealVect&         a_normal,
00296                          const Side::LoHiSide&   a_sd,
00297                          const int&              a_numPrim);
00298 
00299   void
00300   pointExtrapToCovered3D(Vector<Real>&           a_extrapVal,
00301                          const EBFaceFAB&        a_primFace,
00302                          const EBCellFAB&        a_primState,
00303                          const int&              a_faceDir,
00304                          const VolIndex&         a_vof,
00305                          const RealVect&         a_normal,
00306                          const Side::LoHiSide&   a_sd,
00307                          const int&              a_numPrim);
00308 
00309   /// 
00310   /**
00311      This insane version called internally.   The primMinu and primPlus stuff
00312      are really face centered data but because they *came* from cell-centered data
00313      they are left there.   
00314   */
00315   void
00316   extrapToCoveredFaces(BaseIVFAB<Real>&        a_extendedPrim,
00317                        const EBCellFAB&        a_primMinu,
00318                        const EBCellFAB&        a_primPlus,
00319                        const EBCellFAB&        a_primState,
00320                        const Vector<VolIndex>& a_coveredFaces,
00321                        const int&              a_faceDir,
00322                        const Side::LoHiSide&   a_sd,
00323                        const Box&              a_box);
00324 
00325   void
00326   pointExtrapToCovered2D(Vector<Real>&           a_extrapVal,
00327                          const EBCellFAB&        a_primMinu,
00328                          const EBCellFAB&        a_primPlus,
00329                          const EBCellFAB&        a_primState,
00330                          const int&              a_faceDir,
00331                          const VolIndex&         a_vof,
00332                          const RealVect&         a_normal,
00333                          const Side::LoHiSide&   a_sd,
00334                          const int&              a_numPrim);
00335 
00336   void
00337   pointExtrapToCovered3D(Vector<Real>&           a_extrapVal,
00338                          const EBCellFAB&        a_primMinu,
00339                          const EBCellFAB&        a_primPlus,
00340                          const EBCellFAB&        a_primState,
00341                          const int&              a_faceDir,
00342                          const VolIndex&         a_vof,
00343                          const RealVect&         a_normal,
00344                          const Side::LoHiSide&   a_sd,
00345                          const int&              a_numPrim);
00346 
00347 
00348   /// floors if m_isMaxMinSet
00349   void
00350   floorPrimitives(EBCellFAB&  a_primState,
00351                   const Box&  a_box);
00352 
00353   /// floors if m_isMaxMinSet
00354   void
00355   floorPrimitives(BaseIVFAB<Real>&   a_primState,
00356                   const IntVectSet&  a_set);
00357 
00358   void 
00359   coveredExtrapSlopes(Real&            a_dq,
00360                       const VolIndex&  a_vof,
00361                       const EBCellFAB& a_primState,
00362                       const int&       a_dir,
00363                       const int&       a_ivar);
00364 
00365   void
00366   pointGetSlopes(Real&            a_dql,
00367                  Real&            a_dqr,
00368                  Real&            a_dqc,
00369                  bool&            a_hasFacesLeft,
00370                  bool&            a_hasFacesRigh,
00371                  const VolIndex&  a_vof,
00372                  const EBCellFAB& a_primState,
00373                  const int&       a_dir,
00374                  const int&       a_ivar,
00375                  const bool&      a_verbose);
00376 
00377   Real
00378   bilinearFunc(const Real  a_WVal[2][2],
00379                const Real& a_xd1,
00380                const Real& a_xd2);
00381 
00382   void
00383   pointGetSlopesUpwind(Real&            a_dql,
00384                        Real&            a_dqr,
00385                        Real&            a_dqc,
00386                        bool&            a_hasFacesLeft,
00387                        bool&            a_hasFacesRigh,
00388                        const VolIndex&  a_vof,
00389                        const EBCellFAB& a_primState,
00390                        const int&       a_dir,
00391                        const int&       a_ivar,
00392                        const bool&      a_verbose);
00393 
00394   ///and this is the *simplified* version
00395   void
00396   doNormalDerivativeExtr2D(EBCellFAB              a_primMinu[SpaceDim],
00397                            EBCellFAB              a_primPlus[SpaceDim],
00398                            EBFaceFAB              a_fluxOne[SpaceDim],
00399                            BaseIVFAB<Real>        a_coveredFluxNormMinu[SpaceDim],
00400                            BaseIVFAB<Real>        a_coveredFluxNormPlus[SpaceDim],
00401                            Vector<VolIndex>       a_coveredFaceNormMinu[SpaceDim],
00402                            Vector<VolIndex>       a_coveredFaceNormPlus[SpaceDim],
00403                            const EBCellFAB&       a_primState,
00404                            const EBCellFAB&       a_source,
00405                            const DataIndex&       a_dit,
00406                            const Real     &       a_time,
00407                            const Real     &       a_dt);
00408 
00409   ///options for 4th ordeer slopes and flattening removed
00410   void
00411   slope(EBCellFAB&       a_slopePrim,
00412         const EBCellFAB& a_primState,
00413         const int&       a_dir,
00414         const Box&       a_box);
00415 
00416 
00417   void
00418   finalExtrap2D(EBCellFAB              a_primMinu[SpaceDim],
00419                 EBCellFAB              a_primPlus[SpaceDim],
00420                 const BaseIVFAB<Real>  a_coveredFluxNormMinu[SpaceDim],
00421                 const BaseIVFAB<Real>  a_coveredFluxNormPlus[SpaceDim],
00422                 const Vector<VolIndex> a_coveredFaceNormMinu[SpaceDim],
00423                 const Vector<VolIndex> a_coveredFaceNormPlus[SpaceDim],
00424                 const EBFaceFAB        a_fluxOne[SpaceDim],
00425                 const EBCellFAB&       a_primState,
00426                 const Real     &       a_time,
00427                 const Real     &       a_dt);
00428 
00429   //and it goes on like this...
00430   void
00431   finalExtrap3D(EBCellFAB              a_primMinu[SpaceDim],
00432                 EBCellFAB              a_primPlus[SpaceDim],
00433                 const BaseIVFAB<Real>  a_coveredFluxMinu3D[SpaceDim][SpaceDim],
00434                 const BaseIVFAB<Real>  a_coveredFluxPlus3D[SpaceDim][SpaceDim],
00435                 const EBFaceFAB        a_fluxTwo[SpaceDim][SpaceDim],
00436                 const EBCellFAB&       a_primState,
00437                 const Real     &       a_time,
00438                 const Real     &       a_dt);
00439 
00440   void
00441   do111coupling(EBFaceFAB              a_fluxTwo[SpaceDim][SpaceDim],
00442                 BaseIVFAB<Real>        a_coveredFluxMinu3D[SpaceDim][SpaceDim],
00443                 BaseIVFAB<Real>        a_coveredFluxPlus3D[SpaceDim][SpaceDim],
00444                 const EBCellFAB        a_primMinu[SpaceDim],
00445                 const EBCellFAB        a_primPlus[SpaceDim],
00446                 const BaseIVFAB<Real>  a_coveredFluxNormMinu[SpaceDim],
00447                 const BaseIVFAB<Real>  a_coveredFluxNormPlus[SpaceDim],
00448                 const Vector<VolIndex> a_coveredFaceNormMinu[SpaceDim],
00449                 const Vector<VolIndex> a_coveredFaceNormPlus[SpaceDim],
00450                 const EBFaceFAB        a_fluxOne[SpaceDim],
00451                 const EBCellFAB&       a_primState,
00452                 const DataIndex&       a_dit,
00453                 const Real     &       a_time,
00454                 const Real     &       a_dt);
00455 
00456   void
00457   doNormalDerivativeExtr3D(EBCellFAB              a_primMinu[SpaceDim],
00458                            EBCellFAB              a_primPlus[SpaceDim],
00459                            EBFaceFAB              a_fluxOne[SpaceDim],
00460                            BaseIVFAB<Real>        a_coveredFluxNormMinu[SpaceDim],
00461                            BaseIVFAB<Real>        a_coveredFluxNormPlus[SpaceDim],
00462                            Vector<VolIndex>       a_coveredFaceNormMinu[SpaceDim],
00463                            Vector<VolIndex>       a_coveredFaceNormPlus[SpaceDim],
00464                            const EBCellFAB&       a_primState,
00465                            const EBCellFAB&       a_source,
00466                            const DataIndex&       a_dit,
00467                            const Real     &       a_time,
00468                            const Real     &       a_dt);
00469   
00470 
00471   void
00472   normalPred(EBCellFAB&       a_primLo,
00473              EBCellFAB&       a_primHi,
00474              const EBCellFAB& a_primState,
00475              const EBCellFAB& a_slopePrim,
00476              const Real&      a_scale,
00477              const int&       a_dir,
00478              const Box&       a_box);
00479 
00480   void
00481   incrementWithSource(EBCellFAB&       a_primState,
00482                       const EBCellFAB& a_source,
00483                       const Real&      a_scale,
00484                       const Box&       a_box) ;
00485 
00486   void
00487   updatePrim(EBCellFAB&              a_primMinu,
00488              EBCellFAB&              a_primPlus,
00489              const EBFaceFAB&        a_primFace,
00490              const BaseIVFAB<Real>&  a_coveredPrimMinu,
00491              const BaseIVFAB<Real>&  a_coveredPrimPlus,
00492              const Vector<VolIndex>& a_coveredFaceMinu,
00493              const Vector<VolIndex>& a_coveredFacePlus,
00494              const int&              a_faceDir,
00495              const Box&              a_box,
00496              const Real&             a_scale);
00497 
00498   void
00499   upwindSlope(EBCellFAB&       a_slopeUpWi,
00500               const EBCellFAB& a_primState,
00501               const int&       a_dir,
00502               const Box&       a_box);
00503 
00504   void
00505   transversePred(EBCellFAB&       a_rhoLo,
00506                  EBCellFAB&       a_rhoHi,
00507                  const EBCellFAB& a_rho,
00508                  const EBCellFAB& a_dRho,
00509                  const Real&      a_dtbydx,
00510                  const int&       a_dir,
00511                  const Box&       a_box);
00512 };
00513 
00514 #include "NamespaceFooter.H"
00515 #endif