Chombo + EB: AMRTGA.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 _AMRTGA_H_
00012 #define _AMRTGA_H_
00013 
00014 #include <cmath>
00015 #include <cstdlib>
00016 #include <cstdio>
00017 #include <iostream>
00018 #include <iomanip>
00019 #include <fstream>
00020 #include <string>
00021 #include "AMRMultiGrid.H"
00022 #include "NamespaceHeader.H"
00023 
00025 
00031 template <class T>
00032 class TGAHelmOp : public AMRLevelOp<T>
00033 {
00034 public:
00035   TGAHelmOp()
00036   {;}
00037   virtual ~TGAHelmOp()
00038   {;}
00039 
00041   virtual void setAlphaAndBeta(const Real& a_alpha,
00042                                const Real& a_beta) = 0;
00043 
00045 
00052   virtual void diagonalScale(T& a_rhs) = 0;
00053 
00055   virtual void applyOpNoBoundary(T& a_ans, const T& a_phi) = 0;
00056 
00057 
00058 };
00059 
00061 
00063 template <class T, class TFlux>
00064 class LevelTGAHelmOp : public TGAHelmOp< T >
00065 {
00066 public:
00068 
00070   virtual void fillGrad(const T& a_phi) = 0;
00071 
00073   virtual void getFlux(TFlux&           a_flux,
00074                        const T&         a_data,
00075                        const Box&       a_grid,
00076                        const DataIndex& a_dit,
00077                        Real             a_scale) = 0;
00078 
00079   virtual ~LevelTGAHelmOp()
00080   {;}
00081 };
00082 
00083 
00084 
00086 
00090 template <class T>
00091 class AMRTGA
00092 {
00093 public:
00094 
00096   TGAHelmOp<T>* newOp(const ProblemDomain&                                  a_indexSpace,
00097                       const AMRLevelOpFactory<T>&                           a_opFact)
00098   {
00099     AMRLevelOpFactory<T>& opFact = (AMRLevelOpFactory<T>&) a_opFact;
00100     TGAHelmOp<T>* retval = (TGAHelmOp<T>*) opFact.AMRnewOp(a_indexSpace);
00101     return retval;
00102   }
00103 
00105   ~AMRTGA();
00106 
00108 
00110   AMRTGA(const RefCountedPtr<AMRMultiGrid<T> >     &           a_solver,
00111          const AMRLevelOpFactory<T>&                           a_factory,
00112          const ProblemDomain&                                  a_level0Domain,
00113          const Vector<int>&                                    a_refRat,
00114          int a_numLevels = -1,
00115          int a_verbosity = 0);
00116 
00117 
00119 
00122   void oneStep(Vector<T*>&             a_phiNew,
00123                Vector<T*>&             a_phiOld,
00124                Vector<T*>&             a_source,
00125                const Real&             a_dt,
00126                int                     a_lbase,
00127                int                     a_lmax);
00128 
00130   void resetAlphaAndBeta(const Real& a_alpha,
00131                          const Real& a_beta);
00132 
00133 protected:
00134   void solveHelm(Vector<T*>&       a_ans,
00135                  Vector<T*>&       a_rhs,
00136                  int               a_lbase,
00137                  int               a_lmax,
00138                  Real              a_mu,
00139                  Real              a_dt);
00140 
00141   void applyHelm(Vector<T*>&       a_ans,
00142                  Vector<T*>&       a_source,
00143                  int               a_lbase,
00144                  int               a_lmax,
00145                  Real              a_mu,
00146                  Real              a_dt,
00147                  bool              a_homogeneousBC);
00148 
00149   void setMu();
00150 
00151   void createData(Vector<T* >&       a_source,
00152                   int                a_lbase,
00153                   int                a_lmax);
00154 
00155 private:
00156   ProblemDomain                                  m_level0Domain;
00157   Vector<int>                                    m_refRat;
00158   RefCountedPtr<AMRMultiGrid<T> >                m_solver;
00159   Vector<RefCountedPtr<TGAHelmOp<T> > >          m_ops;
00160   Real m_mu1, m_mu2, m_mu3, m_mu4;
00161   int m_verbosity, m_numLevels;
00162   bool                                           m_dataCreated;
00163   Vector<T*>                                     m_rhst;
00164 
00165   //copy constructor and operator= disallowed for all the usual reasons
00166   AMRTGA(const AMRTGA<T>& a_opin)
00167   {
00168     MayDay::Error("invalid operator");
00169   }
00170 
00171   void operator=(const AMRTGA<T>& a_opin)
00172   {
00173     MayDay::Error("invalid operator");
00174   }
00175 
00177   AMRTGA()
00178   {
00179     MayDay::Error("invalid operator");
00180   }
00181 
00182 
00183 };
00184 
00185 /*****/
00186 template <class T>
00187 void AMRTGA<T>::
00188 resetAlphaAndBeta(const Real& a_alpha,
00189                   const Real& a_beta)
00190 {
00191   Vector<MGLevelOp<T>* > ops = m_solver->getAllOperators();
00192   for(int iop = 0; iop < ops.size(); iop++)
00193     {
00194 
00195       TGAHelmOp<T>* helmop = (TGAHelmOp<T>*) ops[iop];
00196       helmop->setAlphaAndBeta(a_alpha, a_beta);
00197     }
00198 }
00199 
00200 /*****/
00201 template <class T>
00202 AMRTGA<T>::
00203 AMRTGA(const RefCountedPtr<AMRMultiGrid<T> >&                a_solver,
00204        const AMRLevelOpFactory<T> &                          a_opFact,
00205        const ProblemDomain&                                  a_level0Domain,
00206        const Vector<int>&                                    a_refRat,
00207        int a_numLevels,
00208        int a_verbosity)
00209 {
00210   //cast to remove const because base class definition is weird
00211   AMRLevelOpFactory<T>& opFact  = (AMRLevelOpFactory<T>&) a_opFact;
00212 
00213   m_verbosity = a_verbosity;
00214   m_level0Domain = a_level0Domain;
00215   m_refRat = a_refRat;
00216   m_solver  = a_solver;
00217   m_numLevels = a_numLevels;
00218   if(m_numLevels < 0)
00219     {
00220       m_numLevels = a_refRat.size();
00221     }
00222   m_ops.resize(m_numLevels);
00223 
00224   ProblemDomain curDom = m_level0Domain;
00225   for(int ilev = 0; ilev < m_numLevels; ilev++)
00226     {
00227       m_ops[ilev] = RefCountedPtr<TGAHelmOp<T> >(newOp(curDom, opFact));
00228       curDom.refine(a_refRat[ilev]);
00229     }
00230   setMu();
00231   m_dataCreated = false;
00232 }
00233 template <class T>
00234 void AMRTGA<T>::
00235 setMu()
00236 {
00237   Real tgaEpsilon = 1.e-12;
00238   Real a = 2.0 - sqrt(2.0) - tgaEpsilon;
00239   m_mu1 = (a - sqrt(pow(a,2) - 4.0*a + 2.0))/2.0 ;
00240   m_mu2 = (a + sqrt(pow(a,2) - 4.0*a + 2.0))/2.0 ;
00241   m_mu3 = (1.0 - a);
00242   m_mu4 = 0.5 - a;
00243   if(m_verbosity > 4)
00244     {
00245       pout() << "   AMRTGA:: epsilon = " << tgaEpsilon << std::endl;
00246     }
00247 }
00248 
00249 template <class T>
00250 void AMRTGA<T>::
00251 createData(Vector<T* >&       a_source,
00252            int                a_lbase,
00253            int                a_lmax)
00254 {
00255   m_dataCreated = true;
00256   m_rhst.resize(a_source.size(), NULL);
00257   for(int ilev = a_lbase; ilev <= a_lmax; ilev++)
00258     {
00259       m_rhst[ilev] = new T();
00260       m_ops[ilev]->create(*m_rhst[ilev], *a_source[ilev]);
00261     }
00262 }
00263 /*****/
00264 template <class T>
00265 AMRTGA<T>::~AMRTGA()
00266 {
00267   for(int ilev = 0; ilev < m_rhst.size(); ilev++)
00268     {
00269       if(m_rhst[ilev] != NULL)
00270         {
00271           delete m_rhst[ilev];
00272           m_rhst[ilev] = NULL;
00273         }
00274     }
00275 }
00276 template <class T>
00277 void AMRTGA<T>::
00278 oneStep(Vector<T* >&       a_phiNew,
00279         Vector<T* >&       a_phiOld,
00280         Vector<T* >&       a_source,
00281         const Real&        a_dt,
00282         int                a_lbase,
00283         int                a_lmax)
00284 {
00285   if(!m_dataCreated)
00286     {
00287       createData(a_source, a_lbase, a_lmax);
00288     }
00289 
00290   if(m_verbosity > 3)
00291     {
00292       pout() << "  AMRTGA:: starting mu4 operation" << std::endl;
00293     }
00294   applyHelm(m_rhst,a_source, a_lbase, a_lmax, m_mu4, a_dt, true);
00295   for(int ilev = a_lbase; ilev <= a_lmax; ilev++)
00296     {
00297       m_ops[ilev]->scale(*m_rhst[ilev], a_dt);
00298     }
00299 
00300   if(m_verbosity > 3)
00301     {
00302       pout() << "  AMRTGA:: starting mu3 operation" << std::endl;
00303     }
00304   applyHelm(a_phiNew,a_phiOld, a_lbase, a_lmax, m_mu3, a_dt, false);
00305 
00306   for(int ilev = a_lbase; ilev <= a_lmax; ilev++)
00307     {
00308       m_ops[ilev]->incr(*m_rhst[ilev], *a_phiNew[ilev], 1.0);
00309     }
00310 
00311   if(m_verbosity > 2)
00312     {
00313       pout() << "  AMRTGA:: starting mu2 operation" << std::endl;
00314     }
00315   solveHelm(a_phiNew,m_rhst, a_lbase, a_lmax, m_mu2, a_dt);
00316   for(int ilev = a_lbase; ilev <= a_lmax; ilev++)
00317     {
00318       m_ops[ilev]->assign(*m_rhst[ilev], *a_phiNew[ilev]);
00319     }
00320 
00321   if(m_verbosity > 2)
00322     {
00323       pout() << "  AMRTGA:: starting mu1 operation" << std::endl;
00324     }
00325   for(int ilev = a_lbase; ilev <= a_lmax; ilev++)
00326     {
00327       m_ops[ilev]->diagonalScale(*m_rhst[ilev]);
00328     }
00329   solveHelm(a_phiNew,m_rhst, a_lbase, a_lmax, m_mu1, a_dt);
00330 }
00331 /*******/
00332 template <class T>
00333 void AMRTGA<T>::
00334 applyHelm(Vector<T*>&      a_ans,
00335           Vector<T*>&      a_phi,
00336           int              a_lbase,
00337           int              a_lmax,
00338           Real             a_mu,
00339           Real             a_dt,
00340           bool             a_homogeneousBC)
00341 {
00342   Real factor  = a_mu*a_dt;
00343 
00344   resetAlphaAndBeta(1.0, factor);
00345 
00346   m_solver->computeAMROperator(a_ans,
00347                                a_phi,
00348                                a_lmax,
00349                                a_lbase,
00350                                a_homogeneousBC);
00351 
00352   for(int ilev = a_lbase; ilev <= a_lmax; ilev++)
00353     {
00354       m_ops[ilev]->scale(*a_ans[ilev], -1);
00355     }
00356 }
00357 /*******/
00358 template <class T>
00359 void AMRTGA<T>::
00360 solveHelm(Vector<T*>&       a_ans,
00361           Vector<T*>&       a_rhs,
00362           int               a_lbase,
00363           int               a_lmax,
00364           Real              a_mu,
00365           Real              a_dt)
00366 {
00367   Real factor  = -a_mu*a_dt;
00368 
00369   resetAlphaAndBeta(1.0, factor);
00370 
00371   m_solver->solveNoInit(a_ans,
00372                         a_rhs,
00373                         a_lmax, a_lbase);
00374   if(m_solver->m_exitStatus==1 && m_verbosity>3)
00375     {
00376       pout() << "AMRTGA:: WARNING: solver exitStatus == 1" << std::endl;
00377     }
00378 }
00379 
00380 #include "NamespaceFooter.H"
00381 #endif