Chombo + EB  3.2
BaseLevelBackwardEuler.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 #ifndef _BASELEVELBACKWARDEULER_H_
12 #define _BASELEVELBACKWARDEULER_H_
13 
14 #include <iostream>
15 #include <math.h>
16 #include "SPACE.H"
17 #include <stdlib.h>
18 #include <REAL.H>
19 #include <Box.H>
20 #include <DisjointBoxLayout.H>
21 #include <LevelData.H>
22 #include <ProblemDomain.H>
23 #include "AMRTGA.H"
24 #include "BaseLevelHeatSolver.H"
25 #include "NamespaceHeader.H"
26 
27 //! \class BaseLevelBackwardEuler
28 //! This base class implements the 1st-order implicit L0-stable time
29 //! integration algorithm for
30 //! solving elliptic equations. It relies upon linear algebraic operations
31 //! defined in the underlying Helmholtz operators.
32 //! \tparam LevelDataType The type used to store data at a grid level.
33 //! This is usually LevelData<T>, where T is some
34 //! cell-centered FArrayBox type.
35 //! \tparam FluxDataType The type used to store flux data at a grid
36 //! level. This is usually an array box clas that stores
37 //! fluxes.
38 //! \tparam FluxRegisterType The type used to store flux register data for
39 //! interactions between different grid levels.
40 //! This is usually a flux register class.
41 template <class LevelDataType,
42  class FluxDataType,
43  class FluxRegisterType>
44 class BaseLevelBackwardEuler : public BaseLevelHeatSolver<LevelDataType, FluxDataType, FluxRegisterType>
45 {
46 
47  public:
48 
49  //! Initializes the base class of a BackwardEuler time integrator. This must be called
50  //! by any subclass of BaseLevelBackwardEuler.
51  //! \param a_grids The DisjointBoxLayout on which the BackwardEuler scheme is to operate.
52  //! \param a_refRatio An array containing the refinement ratios between the
53  //! hierarchical AMR grid levels for the domain.
54  //! \param a_level0Domain The domain at the coarsest AMR grid level.
55  //! \param a_opFact A factory typename LevelDataTypehat is used to generate Helmholtz
56  //! operators to be used by the scheme.
57  //! \param a_solver An AMR Multigrid solver for solving the linear systems
58  //! at each stage of the BackwardEuler integration scheme.
60  const Vector<int>& a_refRat,
61  const ProblemDomain& a_level0Domain,
64  :BaseLevelHeatSolver<LevelDataType,FluxDataType,FluxRegisterType>(a_grids, a_refRat, a_level0Domain, a_opFact, a_solver)
65  {
66  }
67 
68  //! Destructor, called after destructors of BaseLevelBackwardEuler subclasses.
70  {
71  }
72 
73  //! Integrates the helmholtz equation represented by this object, placing
74  //! the new solution in \a a_phiNew.
75  //! \param a_phiNew The new solution (the value of phi at time n + 1) will
76  //! be stored here.
77  //! \param a_phiOld The old solution (the value of phi at time n).
78  //! \param a_src The source term on the right hand side of the Helmholtz
79  //! equation.
80  //! \param a_flux This will store the flux computed at the current grid
81  //! level during the solution of the Helmholtz equation.
82  //! \param a_fineFluxRegPtr A pointer to the flux register representing the
83  //! finer grid level adjacent to this one, or NULL
84  //! if there is no finer grid level.
85  //! \param a_crseFluxRegPtr A pointer to the flux register representing the
86  //! coarser grid level adjacent to this one, or NULL
87  //! if there is no coarser grid level.
88  //! \param a_oldTime The time at the beginning of the integration step at
89  //! the current grid level.
90  //! \param a_crseOldTime The time at the beginning of the integration step
91  //! at the coarser adjacent grid level. This parameter
92  //! is ignored if there is no coarser grid level.
93  //! \param a_crseNewTime The time at the end of the integration step
94  //! at the coarser adjacent grid level. This parameter
95  //! is ignored if there is no coarser grid level.
96  //! \param a_dt The size of the integration step at the current grid level.
97  //! \param a_level The current grid level.
98  //! \param a_zeroPhi If set to true, \a a_phiNew will be set to zero before
99  //! the integration takes place. Otherwise, a_phiNew is
100  //! assumed to be an initial estimate for the solution in
101  //! the iterative linear solve.
102  //! \param a_fluxStartComponent An index identifying the component at which
103  //! flux data begins within \a a_fineFluxRegPtr
104  //! and \a a_crseFluxRegPtr.
105  void updateSoln(LevelDataType& a_phiNew,
106  LevelDataType& a_phiOld,
107  LevelDataType& a_src,
108  LevelData<FluxDataType>& a_flux,
109  FluxRegisterType* a_fineFluxRegPtr,
110  FluxRegisterType* a_crseFluxRegPtr,
111  const LevelDataType* a_crsePhiOldPtr,
112  const LevelDataType* a_crsePhiNewPtr,
113  Real a_oldTime,
114  Real a_crseOldTime,
115  Real a_crseNewTime,
116  Real a_dt,
117  int a_level,
118  bool a_zeroPhi = true,
119  bool a_rhsAlreadyKappaWeighted = false,
120  int a_fluxStartComponent = 0)
121  {
122  CH_assert(!this->m_ops[a_level]->isTimeDependent());
123  int ncomp = a_phiNew.nComp();
124  Interval intervBase(0, ncomp-1);
125  Interval intervFlux(a_fluxStartComponent, a_fluxStartComponent + ncomp-1);
126 
127  CH_assert(a_level >= 0);
128  CH_assert(a_level < this->m_grids.size());
129  CH_assert((a_level == 0) || (a_crsePhiOldPtr != NULL));
130  CH_assert((a_level == 0) || (a_crsePhiNewPtr != NULL));
131  CH_assert(a_crseNewTime >= a_crseOldTime);
132  CH_assert(a_dt >= 0.);
133 
134  LevelDataType rhst, phit;
135  this->m_ops[a_level]->create(rhst, a_src);
136  this->m_ops[a_level]->create(phit, a_phiNew);
137 
138  this->m_ops[a_level]->setToZero(phit);
139  this->m_ops[a_level]->setToZero(rhst);
140  if (a_zeroPhi)
141  {
142  this->m_ops[a_level]->setToZero(a_phiNew);
143  }
144 
145  //set phit to a_phiOld, rhst to a_src * a_dt
146  this->m_ops[a_level]->incr(phit, a_phiOld, 1.0);
147  this->m_ops[a_level]->incr(rhst, a_src , a_dt);
148 
149  //multiply phi old by kappa*acoef
150  this->m_ops[a_level]->diagonalScale(phit, true);
151 
152  //multiply rhs by kappa (but NOT by a)
153  if (!a_rhsAlreadyKappaWeighted)
154  this->m_ops[a_level]->kappaScale(rhst);
155 
156  //add both together to make rhs for euler solve = kappa a phiold + kappa rhs dt
157  this->m_ops[a_level]->incr(rhst, phit, 1.0);
158 
159  //solve for phi new = (kappa a I - kappa dt L) phinew = kappa a phiold + kappa rhs
160  //this makes phinew = (k*a I - mu2 L)^-1 (rhs)
161  LevelDataType coarseData;
162  if ((a_crsePhiOldPtr != NULL) && (a_level > 0))
163  {
164  this->m_ops[a_level-1]->create(coarseData, *a_crsePhiOldPtr);
165  this->m_ops[a_level-1]->setToZero(coarseData);
166  // Backward Euler solves at the new time
167  Real newTime = a_oldTime + a_dt;
168  this->timeInterp(coarseData, *a_crsePhiOldPtr, *a_crsePhiNewPtr,
169  newTime, a_crseOldTime, a_crseNewTime, a_level-1);
170 
171  }
172 
173  this->solveHelm(a_phiNew, coarseData, rhst, a_level, 1.0, a_dt, a_zeroPhi);
174  this->incrementFlux(a_flux, a_phiNew, a_level, 1.0, a_dt, -1.0, true);
175 
176  // now increment flux registers -- note that because of the way
177  // we defined the fluxes, the dt multiplier is already in the
178  // flux
179  if ((a_fineFluxRegPtr != NULL) && (a_level < this->m_grids.size()-1))
180  {
181  Real fluxMult = 1.0;
182  for (DataIterator dit = this->m_grids[a_level].dataIterator(); dit.ok(); ++dit)
183  {
184  FluxDataType& thisFlux = a_flux[dit];
185  for (int dir=0; dir<SpaceDim; ++dir)
186  {
187  a_fineFluxRegPtr->incrementCoarse(thisFlux[dir],
188  fluxMult, dit(),
189  intervBase, // source
190  intervFlux, // dest
191  dir);
192  }
193  }
194  } // end if there is a finer-level
195 
196  if ((a_crseFluxRegPtr != NULL) && (a_level > 0))
197  {
198  Real fluxMult = 1.0;
199 
200  for (DataIterator dit = this->m_grids[a_level].dataIterator(); dit.ok(); ++dit)
201  {
202  FluxDataType& thisFlux = a_flux[dit];
203  for (int dir=0; dir<SpaceDim; ++dir)
204  {
205  a_crseFluxRegPtr->incrementFine(thisFlux[dir], fluxMult, dit(),
206  intervBase, // source
207  intervFlux, // dest
208  dir);
209  }
210  }
211  } // end if there is a coarser level
212  }
213 
214 private:
215  // Disallowed operators.
219 };
220 
221 #include "NamespaceFooter.H"
222 
223 #endif
BaseLevelBackwardEuler(const Vector< DisjointBoxLayout > &a_grids, const Vector< int > &a_refRat, const ProblemDomain &a_level0Domain, RefCountedPtr< AMRLevelOpFactory< LevelDataType > > &a_opFact, const RefCountedPtr< AMRMultiGrid< LevelDataType > > &a_solver)
Definition: BaseLevelBackwardEuler.H:59
A reference-counting handle class.
Definition: RefCountedPtr.H:173
BaseLevelBackwardEuler & operator=(const BaseLevelBackwardEuler &)
#define CH_assert(cond)
Definition: CHArray.H:37
A class to facilitate interaction with physical boundary conditions.
Definition: ProblemDomain.H:141
void updateSoln(LevelDataType &a_phiNew, LevelDataType &a_phiOld, LevelDataType &a_src, LevelData< FluxDataType > &a_flux, FluxRegisterType *a_fineFluxRegPtr, FluxRegisterType *a_crseFluxRegPtr, const LevelDataType *a_crsePhiOldPtr, const LevelDataType *a_crsePhiNewPtr, Real a_oldTime, Real a_crseOldTime, Real a_crseNewTime, Real a_dt, int a_level, bool a_zeroPhi=true, bool a_rhsAlreadyKappaWeighted=false, int a_fluxStartComponent=0)
Definition: BaseLevelBackwardEuler.H:105
Definition: BaseLevelHeatSolver.H:43
Definition: BaseLevelBackwardEuler.H:44
Definition: DataIterator.H:190
const int SpaceDim
Definition: SPACE.H:38
Structure for passing component ranges in code.
Definition: Interval.H:23
Vector< DisjointBoxLayout > m_grids
The disjoint box layouts at every AMR grid level.
Definition: BaseLevelHeatSolver.H:504
new code
Definition: BoxLayoutData.H:170
virtual ~BaseLevelBackwardEuler()
Destructor, called after destructors of BaseLevelBackwardEuler subclasses.
Definition: BaseLevelBackwardEuler.H:69
double Real
Definition: REAL.H:33
void timeInterp(LevelDataType &a_data, const LevelDataType &a_oldData, const LevelDataType &a_newData, Real a_time, Real a_oldTime, Real a_newTime, int a_level)
Definition: BaseLevelHeatSolver.H:478
void solveHelm(LevelDataType &a_phi, LevelDataType &a_phiC, LevelDataType &a_rhs, int a_level, Real a_mu, Real a_dt, bool a_zeroPhi=true)
Definition: BaseLevelHeatSolver.H:383
size_t size() const
Definition: Vector.H:192
Definition: AMRMultiGrid.H:233
void incrementFlux(LevelData< FluxDataType > &a_diffusiveFlux, LevelDataType &a_phi, int a_level, Real a_mu, Real a_dt, Real a_sign, bool a_setToZero)
Definition: BaseLevelHeatSolver.H:335
Vector< LevelTGAHelmOp< LevelDataType, FluxDataType > * > m_ops
Definition: BaseLevelHeatSolver.H:515