Chombo + EB + MF: 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 
00228   /** a_correction += I[h->h](a_coarseCorrection) */
00229   virtual void AMRProlong(LevelData<FArrayBox>& a_correction,
00230                           const LevelData<FArrayBox>& a_coarseCorrection);
00231 
00232   /** a_residual = a_residual - L(a_correction, a_coarseCorrection) */
00233   virtual void AMRUpdateResidual(LevelData<FArrayBox>& a_residual,
00234                                  const LevelData<FArrayBox>& a_correction,
00235                                  const LevelData<FArrayBox>& a_coarseCorrection);
00236 
00237   ///
00238   /**
00239      compute norm over all cells on coarse not covered by finer
00240   */
00241   virtual Real AMRNorm(const LevelData<FArrayBox>& a_coarseResid,
00242                        const LevelData<FArrayBox>& a_fineResid,
00243                        const int&                  a_refRat,
00244                        const int&                  a_ord);
00245 
00246   void homogeneousCFInterp(LevelData<FArrayBox>& a_phif);
00247  ///
00248   /**
00249      does homogeneous coarse/fine interpolation for phi
00250   */
00251   void
00252   homogeneousCFInterpPhi(LevelData<FArrayBox>& a_phif,
00253                          const DataIndex& a_datInd,
00254                          int a_idir,
00255                          Side::LoHiSide a_hiorlo);
00256 
00257   ///
00258   /** does homogeneous coarse/fine interpolation for tangential gradient
00259       (needs phi ghost cells to be filled in first, so should call
00260       homogeneousCFInterpPhi first)
00261   */
00262   void homogeneousCFInterpTanGrad(LevelData<FArrayBox>& a_tanGrad,
00263                                   const LevelData<FArrayBox>& a_phi,
00264                                   const DataIndex& a_datInd,
00265                                   int a_idir,
00266                                   Side::LoHiSide a_hiorlo);
00267 
00268   void interpOnIVSHomo(LevelData<FArrayBox>& a_phif,
00269                        const DataIndex& a_datInd,
00270                        const int a_idir,
00271                        const Side::LoHiSide a_hiorlo,
00272                        const IntVectSet& a_interpIVS);
00273   ///
00274   /**
00275      Apply the AMR operator, including coarse-fine matching
00276   */
00277   void AMROperator(      LevelData<FArrayBox>& a_LofPhi,
00278                    const LevelData<FArrayBox>& a_phiFine,
00279                    const LevelData<FArrayBox>& a_phi,
00280                    const LevelData<FArrayBox>& a_phiCoarse,
00281                    bool a_homogeneousDomBC,
00282                    AMRLevelOp<LevelData<FArrayBox> >*  a_finerOp);
00283 
00284   void AMROperatorNF(      LevelData<FArrayBox>& a_LofPhi,
00285                      const LevelData<FArrayBox>& a_phi,
00286                      const LevelData<FArrayBox>& a_phiCoarse,
00287                      bool a_homogeneousBC);
00288 
00289   virtual void AMROperatorNC(      LevelData<FArrayBox>& a_LofPhi,
00290                              const LevelData<FArrayBox>& a_phiFine,
00291                              const LevelData<FArrayBox>& a_phi,
00292                              bool a_homogeneousBC,
00293                              AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00294 
00295   void cellGrad(FArrayBox&             a_gradPhi,
00296                 const  FArrayBox&      a_phi,
00297                 const Box&             a_grid);
00298 
00299   void cfinterp(const LevelData<FArrayBox>& a_phiFine,
00300                 const LevelData<FArrayBox>& a_phiCoarse);
00301 
00302   virtual void fillGrad(const LevelData<FArrayBox>& a_phiFine);
00303 
00304   void loHiCenterFace(Box&                 a_loBox,
00305                       int&                 a_hasLo,
00306                       Box&                 a_hiBox,
00307                       int&                 a_hasHi,
00308                       Box&                 a_centerBox,
00309                       const ProblemDomain& a_eblg,
00310                       const Box&           a_inBox,
00311                       const int&           a_dir);
00312 
00313   void getFaceDivAndGrad(FArrayBox&             a_faceDiv,
00314                          FArrayBox&             a_faceGrad,
00315                          const FArrayBox&       a_data,
00316                          const FArrayBox&       a_gradData,
00317                          const ProblemDomain&   a_domain,
00318                          const Box&             a_faceBox,
00319                          const int&             a_faceDir,
00320                          const Real             a_dx);
00321 
00322   ///take cell centered divergence of the inputs.
00323   /**
00324       not part of the operator evaluation.  this is to
00325       test whether the divergence of b converges at h^2
00326       as b does if b is analytically divergence-free.
00327    */
00328   void divergenceCC(LevelData<FArrayBox>&       a_div,
00329                     const LevelData<FArrayBox>& a_phi,
00330                     const LevelData<FArrayBox>* a_phiC);
00331 protected:
00332   void setLambda();
00333   RefCountedPtr<LevelData<FluxBox> >         m_eta;
00334   DisjointBoxLayout                          m_grids;
00335   Real                                       m_alpha;
00336   Real                                       m_beta;
00337   int                                        m_refToCoar;
00338   int                                        m_refToFine;
00339   BCFunc                                     m_bc;
00340 
00341   //b is a scalar in 2d, vector in 3d
00342   static const int        s_nComp; //==3
00343   static const int        s_nGradComp; //==3*SpaceDim
00344   Real                    m_dx;
00345   Real                    m_dxCrse;
00346   ProblemDomain           m_domain;
00347   //relaxation coef
00348   LevelData<FArrayBox>    m_lambda;
00349   LevelData<FArrayBox>    m_grad;
00350   LevelDataOps<FArrayBox> m_levelOps;
00351   Copier                  m_exchangeCopier;
00352   TensorCFInterp          m_interpWithCoarser;
00353 
00354   LayoutData<CFIVS> m_loCFIVS[SpaceDim];
00355   LayoutData<CFIVS> m_hiCFIVS[SpaceDim];
00356   // will need these to do tangential gradient computations
00357   LayoutData<TensorFineStencilSet> m_hiTanStencilSets[SpaceDim];
00358   LayoutData<TensorFineStencilSet> m_loTanStencilSets[SpaceDim];
00359   Vector<IntVect> m_colors;
00360 private:
00361   ///weak construction is bad
00362   ResistivityOp()
00363   {
00364     MayDay::Error("invalid operator");
00365   }
00366   //copy constructor and operator= disallowed for all the usual reasons
00367   ResistivityOp(const ResistivityOp& a_opin)
00368   {
00369     MayDay::Error("invalid operator");
00370   }
00371   void operator=(const ResistivityOp& a_opin)
00372   {
00373     MayDay::Error("invalid operator");
00374   }
00375 };
00376 
00377 ///
00378 /**
00379    Factory to create ResistivityOps
00380  */
00381 class ResistivityOpFactory: public AMRLevelOpFactory<LevelData< FArrayBox> >
00382 {
00383 public:
00384   virtual ~ResistivityOpFactory()
00385   {
00386   }
00387 
00388   ResistivityOpFactory(const Vector<DisjointBoxLayout>&                     a_grids,
00389                        const Vector<RefCountedPtr<LevelData<FluxBox> > >&   a_eta,
00390                        Real                                                 a_alpha,
00391                        Real                                                 a_beta,
00392                        const Vector<int>&                                   a_refRatios,
00393                        const ProblemDomain&                                 a_domainCoar,
00394                        const Real&                                          a_dxCoar,
00395                        BCFunc                                               a_bc);
00396 
00397   ///
00398   virtual MGLevelOp< LevelData<FArrayBox> >*
00399   MGnewOp(const ProblemDomain& a_FineindexSpace,
00400           int depth,
00401           bool homoOnly = true);
00402 
00403   ///
00404   virtual AMRLevelOp< LevelData<FArrayBox> >*
00405   AMRnewOp(const ProblemDomain& a_indexSpace);
00406 
00407   ///
00408   virtual int refToFiner(const ProblemDomain&) const;
00409 
00410 private:
00411   Vector<RefCountedPtr<LevelData<FluxBox> > >    m_eta;
00412   LevelData<FArrayBox>                           m_phic;
00413   Vector<ProblemDomain>                          m_domains;
00414   Vector<DisjointBoxLayout>                      m_boxes;
00415   Vector<Real>                                   m_dx;
00416   // refinement to next coarser level
00417   Vector<int>                                    m_refRatios;
00418   BCFunc  m_bc;
00419   Real m_alpha, m_beta;
00420 
00421   ///weak construction is bad
00422   ResistivityOpFactory()
00423   {
00424     MayDay::Error("invalid operator");
00425   }
00426   //copy constructor and operator= disallowed for all the usual reasons
00427   ResistivityOpFactory(const ResistivityOpFactory& a_opin)
00428   {
00429     MayDay::Error("invalid operator");
00430   }
00431 
00432   void operator=(const ResistivityOpFactory& a_opin)
00433   {
00434     MayDay::Error("invalid operator");
00435   }
00436 };
00437 
00438 #include "NamespaceFooter.H"
00439 #endif