Chombo + EB: ResistivityOp.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 // dtg 2007
00012 
00013 #ifndef _RESISTIVITYOP_H_
00014 #define _RESISTIVITYOP_H_
00015 
00016 #include "AMRMultiGrid.H"
00017 #include "REAL.H"
00018 #include "Box.H"
00019 #include "LevelDataOps.H"
00020 #include "BCFunc.H"
00021 #include "FArrayBox.H"
00022 #include "FluxBox.H"
00023 #include "CFIVS.H"
00024 #include "TensorCFInterp.H"
00025 #include "CoarseAverage.H"
00026 #include "AMRTGA.H"
00027 #include "NamespaceHeader.H"
00028 
00029 ///
00030 /**
00031    operator L is defined as
00032    L(B) = alpha B + beta*(divF)
00033    where
00034    F = eta((grad B - grad B^T) + I  (div B))
00035    alpha and beta are constants.  eta is a function of space.
00036    beta and eta are incorporated into the flux F when getFlux is called.
00037    This is always 3 variables.  SpaceDim refers to only the derivatives
00038    that are meaningful in this bizzare application.
00039 */
00040 class ResistivityOp : public LevelTGAHelmOp<LevelData<FArrayBox>, FluxBox>
00041 {
00042 public:
00043 
00044   virtual void setAlphaAndBeta(const Real& a_alpha,
00045                                const Real& a_beta)
00046   {
00047     m_alpha = a_alpha;
00048     m_beta = a_beta;
00049     setLambda();
00050   }
00051 
00052   virtual void diagonalScale(LevelData<FArrayBox>& a_rhs, bool a_kappaWeighted)
00053   {
00054     diagonalScale(a_rhs);
00055   }
00056 
00057   virtual void diagonalScale(LevelData<FArrayBox>& a_rhs)
00058   {
00059     //does not apply here
00060   }
00061 
00062   virtual void divideByIdentityCoef(LevelData<FArrayBox>& a_rhs)
00063   {
00064     //does not apply here
00065   }
00066 
00067   ///
00068   /**
00069    */
00070   virtual ~ResistivityOp()
00071   {
00072   }
00073 
00074   //use at thy peril
00075   virtual void getFlux(FluxBox&             a_flux,
00076                const LevelData<FArrayBox>&  a_data,
00077                const Box&       a_grid,
00078                const DataIndex& a_dit,
00079                Real             a_scale);
00080 
00081   ///
00082   /**
00083    */
00084   ResistivityOp(const DisjointBoxLayout&                    a_grids,
00085                 const DisjointBoxLayout&                    a_gridsFine,
00086                 const DisjointBoxLayout&                    a_gridsCoar,
00087                 const RefCountedPtr<LevelData<FluxBox> >&   a_eta,
00088                 Real                                        a_alpha,
00089                 Real                                        a_beta,
00090                 int                                         a_refToFine,
00091                 int                                         a_refToCoar,
00092                 const ProblemDomain&                        a_domain,
00093                 const Real&                                 a_dxLevel,
00094                 const Real&                                 a_dxCoar,
00095                 BCFunc                                      a_bc);
00096 
00097   virtual void residual(  LevelData<FArrayBox>& a_lhs,
00098                           const LevelData<FArrayBox>& a_phi,
00099                           const LevelData<FArrayBox>& a_rhs,
00100                           bool a_homogeneous = false);
00101 
00102   virtual void preCond(   LevelData<FArrayBox>& a_correction,
00103                           const LevelData<FArrayBox>& a_residual);
00104 
00105   virtual void applyOp(   LevelData<FArrayBox>& a_lhs,
00106                           const LevelData<FArrayBox>& a_phi,
00107                           bool a_homogeneous = false);
00108 
00109   virtual void applyOpNoBoundary( LevelData<FArrayBox>& a_lhs,
00110                                   const LevelData<FArrayBox>& a_phi);
00111 
00112   virtual void create(    LevelData<FArrayBox>& a_lhs,
00113                           const LevelData<FArrayBox>& a_rhs);
00114   virtual void createCoarsened(    LevelData<FArrayBox>& a_lhs,
00115                                    const LevelData<FArrayBox>& a_rhs,
00116                                    const int& a_refRat);
00117 
00118   void reflux(const LevelData<FArrayBox>& a_phiFine,
00119               const LevelData<FArrayBox>& a_phi,
00120               LevelData<FArrayBox>& residual,
00121               AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00122 
00123   void getFlux(FArrayBox&       a_flux,
00124                const FArrayBox& a_data,
00125                const FArrayBox& a_gradData,
00126                const FArrayBox& a_etaFace,
00127                const Box&       a_facebox,
00128                int              a_dir,
00129                int ref = 1);
00130 
00131   /// utility function which computes operator after all bc's have been set
00132   void computeOperatorNoBCs(LevelData<FArrayBox>& a_lhs,
00133                             const LevelData<FArrayBox>& a_phi);
00134 
00135   virtual void assign(    LevelData<FArrayBox>& a_lhs,
00136                           const LevelData<FArrayBox>& a_rhs) ;
00137   virtual Real dotProduct(const LevelData<FArrayBox>& a_1,
00138                           const LevelData<FArrayBox>& a_2) ;
00139   virtual void incr( LevelData<FArrayBox>& a_lhs,
00140                      const LevelData<FArrayBox>& a_x,
00141                      Real a_scale) ;
00142   virtual void axby( LevelData<FArrayBox>& a_lhs,
00143                      const LevelData<FArrayBox>& a_x,
00144                      const LevelData<FArrayBox>& a_y,
00145                      Real a, Real b) ;
00146 
00147   virtual void scale(LevelData<FArrayBox>& a_lhs, const Real& a_scale) ;
00148 
00149   virtual Real norm(const LevelData<FArrayBox>& a_x, int a_ord);
00150 
00151   virtual void setToZero( LevelData<FArrayBox>& a_x);
00152   /*@}*/
00153 
00154   /**
00155      \name MGLevelOp functions */
00156   /*@{*/
00157 
00158   virtual void relax(LevelData<FArrayBox>& a_e,
00159                      const LevelData<FArrayBox>& a_residual,
00160                      int iterations);
00161 
00162   virtual void createCoarser(LevelData<FArrayBox>& a_coarse,
00163                              const LevelData<FArrayBox>& a_fine,
00164                              bool ghosted);
00165   /**
00166      calculate restricted residual
00167      a_resCoarse[2h] = I[h->2h] (rhsFine[h] - L[h](phiFine[h])
00168   */
00169   virtual void restrictResidual(LevelData<FArrayBox>& a_resCoarse,
00170                                 LevelData<FArrayBox>& a_phiFine,
00171                                 const LevelData<FArrayBox>& a_rhsFine);
00172 
00173   /**
00174      correct the fine solution based on coarse correction
00175      a_phiThisLevel += I[2h->h](a_correctCoarse)
00176   */
00177   virtual void prolongIncrement(LevelData<FArrayBox>& a_phiThisLevel,
00178                                 const LevelData<FArrayBox>& a_correctCoarse);
00179 
00180   /*@}*/
00181 
00182   /**
00183      \name AMRLevelOp functions */
00184   /*@{*/
00185 
00186   /** returns 1 when there are no coarser AMRLevelOp objects */
00187   virtual int refToCoarser()
00188   {
00189     return m_refToCoar;
00190   }
00191 
00192   /** a_residual = a_rhs - L(a_phi, a_phiFine, a_phiCoarse) */
00193   virtual void AMRResidual(LevelData<FArrayBox>& a_residual,
00194                            const LevelData<FArrayBox>& a_phiFine,
00195                            const LevelData<FArrayBox>& a_phi,
00196                            const LevelData<FArrayBox>& a_phiCoarse,
00197                            const LevelData<FArrayBox>& a_rhs,
00198                            bool a_homogeneousPhysBC,
00199                            AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00200 
00201   /** residual assuming no more coarser AMR levels */
00202 
00203   virtual void AMRResidualNC(LevelData<FArrayBox>& a_residual,
00204                              const LevelData<FArrayBox>& a_phiFine,
00205                              const LevelData<FArrayBox>& a_phi,
00206                              const LevelData<FArrayBox>& a_rhs,
00207                              bool a_homogeneousPhysBC,
00208                              AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00209 
00210   /** a_residual = a_rhs - L(a_phi, a_phiCoarse)  */
00211   virtual void AMRResidualNF(LevelData<FArrayBox>& a_residual,
00212                              const LevelData<FArrayBox>& a_phi,
00213                              const LevelData<FArrayBox>& a_phiCoarse,
00214                              const LevelData<FArrayBox>& a_rhs,
00215                              bool a_homogeneousPhysBC);
00216 
00217   /** a_resCoarse = I[h-2h]( a_residual - L(a_correction, a_coarseCorrection))
00218       it is assumed that a_resCoarse has already been filled in with the coarse
00219       version of AMRResidualNF and that this operation is free to overwrite
00220       in the overlap regions.
00221   */
00222 
00223   virtual void AMRRestrict(LevelData<FArrayBox>& a_resCoarse,
00224                            const LevelData<FArrayBox>& a_residual,
00225                            const LevelData<FArrayBox>& a_correction,
00226                            const LevelData<FArrayBox>& a_coarseCorrection, 
00227                            bool a_skip_res = false );
00228 
00229   /** a_correction += I[h->h](a_coarseCorrection) */
00230   virtual void AMRProlong(LevelData<FArrayBox>& a_correction,
00231                           const LevelData<FArrayBox>& a_coarseCorrection);
00232 
00233   /** a_residual = a_residual - L(a_correction, a_coarseCorrection) */
00234   virtual void AMRUpdateResidual(LevelData<FArrayBox>& a_residual,
00235                                  const LevelData<FArrayBox>& a_correction,
00236                                  const LevelData<FArrayBox>& a_coarseCorrection);
00237 
00238   ///
00239   /**
00240      compute norm over all cells on coarse not covered by finer
00241   */
00242   virtual Real AMRNorm(const LevelData<FArrayBox>& a_coarseResid,
00243                        const LevelData<FArrayBox>& a_fineResid,
00244                        const int&                  a_refRat,
00245                        const int&                  a_ord);
00246 
00247   void homogeneousCFInterp(LevelData<FArrayBox>& a_phif);
00248  ///
00249   /**
00250      does homogeneous coarse/fine interpolation for phi
00251   */
00252   void
00253   homogeneousCFInterpPhi(LevelData<FArrayBox>& a_phif,
00254                          const DataIndex& a_datInd,
00255                          int a_idir,
00256                          Side::LoHiSide a_hiorlo);
00257 
00258   ///
00259   /** does homogeneous coarse/fine interpolation for tangential gradient
00260       (needs phi ghost cells to be filled in first, so should call
00261       homogeneousCFInterpPhi first)
00262   */
00263   void homogeneousCFInterpTanGrad(LevelData<FArrayBox>& a_tanGrad,
00264                                   const LevelData<FArrayBox>& a_phi,
00265                                   const DataIndex& a_datInd,
00266                                   int a_idir,
00267                                   Side::LoHiSide a_hiorlo);
00268 
00269   void interpOnIVSHomo(LevelData<FArrayBox>& a_phif,
00270                        const DataIndex& a_datInd,
00271                        const int a_idir,
00272                        const Side::LoHiSide a_hiorlo,
00273                        const IntVectSet& a_interpIVS);
00274   ///
00275   /**
00276      Apply the AMR operator, including coarse-fine matching
00277   */
00278   void AMROperator(      LevelData<FArrayBox>& a_LofPhi,
00279                    const LevelData<FArrayBox>& a_phiFine,
00280                    const LevelData<FArrayBox>& a_phi,
00281                    const LevelData<FArrayBox>& a_phiCoarse,
00282                    bool a_homogeneousDomBC,
00283                    AMRLevelOp<LevelData<FArrayBox> >*  a_finerOp);
00284 
00285   void AMROperatorNF(      LevelData<FArrayBox>& a_LofPhi,
00286                      const LevelData<FArrayBox>& a_phi,
00287                      const LevelData<FArrayBox>& a_phiCoarse,
00288                      bool a_homogeneousBC);
00289 
00290   virtual void AMROperatorNC(      LevelData<FArrayBox>& a_LofPhi,
00291                              const LevelData<FArrayBox>& a_phiFine,
00292                              const LevelData<FArrayBox>& a_phi,
00293                              bool a_homogeneousBC,
00294                              AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00295 
00296   void cellGrad(FArrayBox&             a_gradPhi,
00297                 const  FArrayBox&      a_phi,
00298                 const Box&             a_grid);
00299 
00300   void cfinterp(const LevelData<FArrayBox>& a_phiFine,
00301                 const LevelData<FArrayBox>& a_phiCoarse);
00302 
00303   virtual void fillGrad(const LevelData<FArrayBox>& a_phiFine);
00304 
00305   void loHiCenterFace(Box&                 a_loBox,
00306                       int&                 a_hasLo,
00307                       Box&                 a_hiBox,
00308                       int&                 a_hasHi,
00309                       Box&                 a_centerBox,
00310                       const ProblemDomain& a_eblg,
00311                       const Box&           a_inBox,
00312                       const int&           a_dir);
00313 
00314   void getFaceDivAndGrad(FArrayBox&             a_faceDiv,
00315                          FArrayBox&             a_faceGrad,
00316                          const FArrayBox&       a_data,
00317                          const FArrayBox&       a_gradData,
00318                          const ProblemDomain&   a_domain,
00319                          const Box&             a_faceBox,
00320                          const int&             a_faceDir,
00321                          const Real             a_dx);
00322 
00323   ///take cell centered divergence of the inputs.
00324   /**
00325       not part of the operator evaluation.  this is to
00326       test whether the divergence of b converges at h^2
00327       as b does if b is analytically divergence-free.
00328    */
00329   void divergenceCC(LevelData<FArrayBox>&       a_div,
00330                     const LevelData<FArrayBox>& a_phi,
00331                     const LevelData<FArrayBox>* a_phiC);
00332 protected:
00333   void setLambda();
00334   RefCountedPtr<LevelData<FluxBox> >         m_eta;
00335   DisjointBoxLayout                          m_grids;
00336   Real                                       m_alpha;
00337   Real                                       m_beta;
00338   int                                        m_refToCoar;
00339   int                                        m_refToFine;
00340   BCFunc                                     m_bc;
00341 
00342   //b is a scalar in 2d, vector in 3d
00343   static const int        s_nComp; //==3
00344   static const int        s_nGradComp; //==3*SpaceDim
00345   Real                    m_dx;
00346   Real                    m_dxCrse;
00347   ProblemDomain           m_domain;
00348   //relaxation coef
00349   LevelData<FArrayBox>    m_lambda;
00350   LevelData<FArrayBox>    m_grad;
00351   LevelDataOps<FArrayBox> m_levelOps;
00352   Copier                  m_exchangeCopier;
00353   TensorCFInterp          m_interpWithCoarser;
00354 
00355   LayoutData<CFIVS> m_loCFIVS[SpaceDim];
00356   LayoutData<CFIVS> m_hiCFIVS[SpaceDim];
00357   // will need these to do tangential gradient computations
00358   LayoutData<TensorFineStencilSet> m_hiTanStencilSets[SpaceDim];
00359   LayoutData<TensorFineStencilSet> m_loTanStencilSets[SpaceDim];
00360   Vector<IntVect> m_colors;
00361 private:
00362   ///weak construction is bad
00363   ResistivityOp()
00364   {
00365     MayDay::Error("invalid operator");
00366   }
00367   //copy constructor and operator= disallowed for all the usual reasons
00368   ResistivityOp(const ResistivityOp& a_opin)
00369   {
00370     MayDay::Error("invalid operator");
00371   }
00372   void operator=(const ResistivityOp& a_opin)
00373   {
00374     MayDay::Error("invalid operator");
00375   }
00376 };
00377 
00378 ///
00379 /**
00380    Factory to create ResistivityOps
00381  */
00382 class ResistivityOpFactory: public AMRLevelOpFactory<LevelData< FArrayBox> >
00383 {
00384 public:
00385   virtual ~ResistivityOpFactory()
00386   {
00387   }
00388 
00389   ResistivityOpFactory(const Vector<DisjointBoxLayout>&                     a_grids,
00390                        const Vector<RefCountedPtr<LevelData<FluxBox> > >&   a_eta,
00391                        Real                                                 a_alpha,
00392                        Real                                                 a_beta,
00393                        const Vector<int>&                                   a_refRatios,
00394                        const ProblemDomain&                                 a_domainCoar,
00395                        const Real&                                          a_dxCoar,
00396                        BCFunc                                               a_bc);
00397 
00398   ///
00399   virtual MGLevelOp< LevelData<FArrayBox> >*
00400   MGnewOp(const ProblemDomain& a_FineindexSpace,
00401           int depth,
00402           bool homoOnly = true);
00403 
00404   ///
00405   virtual AMRLevelOp< LevelData<FArrayBox> >*
00406   AMRnewOp(const ProblemDomain& a_indexSpace);
00407 
00408   ///
00409   virtual int refToFiner(const ProblemDomain&) const;
00410 
00411 private:
00412   Vector<RefCountedPtr<LevelData<FluxBox> > >    m_eta;
00413   LevelData<FArrayBox>                           m_phic;
00414   Vector<ProblemDomain>                          m_domains;
00415   Vector<DisjointBoxLayout>                      m_boxes;
00416   Vector<Real>                                   m_dx;
00417   // refinement to next coarser level
00418   Vector<int>                                    m_refRatios;
00419   BCFunc  m_bc;
00420   Real m_alpha, m_beta;
00421 
00422   ///weak construction is bad
00423   ResistivityOpFactory()
00424   {
00425     MayDay::Error("invalid operator");
00426   }
00427   //copy constructor and operator= disallowed for all the usual reasons
00428   ResistivityOpFactory(const ResistivityOpFactory& a_opin)
00429   {
00430     MayDay::Error("invalid operator");
00431   }
00432 
00433   void operator=(const ResistivityOpFactory& a_opin)
00434   {
00435     MayDay::Error("invalid operator");
00436   }
00437 };
00438 
00439 #include "NamespaceFooter.H"
00440 #endif