Chombo + EB  3.0
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 
28 //! \class BaseLevelBackwardEuler
29 //! This base class implements the 1st-order implicit L0-stable time
30 //! integration algorithm for
31 //! solving elliptic equations. It relies upon linear algebraic operations
32 //! defined in the underlying Helmholtz operators.
33 //! \tparam LevelDataType The type used to store data at a grid level.
34 //! This is usually LevelData<T>, where T is some
35 //! cell-centered FArrayBox type.
36 //! \tparam FluxDataType The type used to store flux data at a grid
37 //! level. This is usually an array box clas that stores
38 //! fluxes.
39 //! \tparam FluxRegisterType The type used to store flux register data for
40 //! interactions between different grid levels.
41 //! This is usually a flux register class.
42 template <class LevelDataType,
43  class FluxDataType,
44  class FluxRegisterType>
45 class BaseLevelBackwardEuler : public BaseLevelHeatSolver<LevelDataType, FluxDataType, FluxRegisterType>
46 {
47 
48  public:
49 
50  //! Initializes the base class of a BackwardEuler time integrator. This must be called
51  //! by any subclass of BaseLevelBackwardEuler.
52  //! \param a_grids The DisjointBoxLayout on which the BackwardEuler scheme is to operate.
53  //! \param a_refRatio An array containing the refinement ratios between the
54  //! hierarchical AMR grid levels for the domain.
55  //! \param a_level0Domain The domain at the coarsest AMR grid level.
56  //! \param a_opFact A factory typename LevelDataTypehat is used to generate Helmholtz
57  //! operators to be used by the scheme.
58  //! \param a_solver An AMR Multigrid solver for solving the linear systems
59  //! at each stage of the BackwardEuler integration scheme.
61  const Vector<int>& a_refRat,
62  const ProblemDomain& a_level0Domain,
65  :BaseLevelHeatSolver<LevelDataType,FluxDataType,FluxRegisterType>(a_grids, a_refRat, a_level0Domain, a_opFact, a_solver)
66  {
67  }
68 
69  //! Destructor, called after destructors of BaseLevelBackwardEuler subclasses.
71  {
72  }
73 
74  //! Integrates the helmholtz equation represented by this object, placing
75  //! the new solution in \a a_phiNew.
76  //! \param a_phiNew The new solution (the value of phi at time n + 1) will
77  //! be stored here.
78  //! \param a_phiOld The old solution (the value of phi at time n).
79  //! \param a_src The source term on the right hand side of the Helmholtz
80  //! equation.
81  //! \param a_flux This will store the flux computed at the current grid
82  //! level during the solution of the Helmholtz equation.
83  //! \param a_fineFluxRegPtr A pointer to the flux register representing the
84  //! finer grid level adjacent to this one, or NULL
85  //! if there is no finer grid level.
86  //! \param a_crseFluxRegPtr A pointer to the flux register representing the
87  //! coarser grid level adjacent to this one, or NULL
88  //! if there is no coarser grid level.
89  //! \param a_oldTime The time at the beginning of the integration step at
90  //! the current grid level.
91  //! \param a_crseOldTime The time at the beginning of the integration step
92  //! at the coarser adjacent grid level. This parameter
93  //! is ignored if there is no coarser grid level.
94  //! \param a_crseNewTime The time at the end of the integration step
95  //! at the coarser adjacent grid level. This parameter
96  //! is ignored if there is no coarser grid level.
97  //! \param a_dt The size of the integration step at the current grid level.
98  //! \param a_level The current grid level.
99  //! \param a_zeroPhi If set to true, \a a_phiNew will be set to zero before
100  //! the integration takes place. Otherwise, a_phiNew is
101  //! assumed to be an initial estimate for the solution in
102  //! the iterative linear solve.
103  //! \param a_fluxStartComponent An index identifying the component at which
104  //! flux data begins within \a a_fineFluxRegPtr
105  //! and \a a_crseFluxRegPtr.
106  void updateSoln(LevelDataType& a_phiNew,
107  LevelDataType& a_phiOld,
108  LevelDataType& a_src,
109  LevelData<FluxDataType>& a_flux,
110  FluxRegisterType* a_fineFluxRegPtr,
111  FluxRegisterType* a_crseFluxRegPtr,
112  const LevelDataType* a_crsePhiOldPtr,
113  const LevelDataType* a_crsePhiNewPtr,
114  Real a_oldTime,
115  Real a_crseOldTime,
116  Real a_crseNewTime,
117  Real a_dt,
118  int a_level,
119  bool a_zeroPhi = true,
120  bool a_rhsAlreadyKappaWeighted = false,
121  int a_fluxStartComponent = 0)
122  {
123  CH_assert(!this->m_ops[a_level]->isTimeDependent());
124  int ncomp = a_phiNew.nComp();
125  Interval intervBase(0, ncomp-1);
126  Interval intervFlux(a_fluxStartComponent, a_fluxStartComponent + ncomp-1);
127 
128  CH_assert(a_level >= 0);
129  CH_assert(a_level < this->m_grids.size());
130  CH_assert((a_level == 0) || (a_crsePhiOldPtr != NULL));
131  CH_assert((a_level == 0) || (a_crsePhiNewPtr != NULL));
132  CH_assert(a_crseNewTime >= a_crseOldTime);
133  CH_assert(a_dt >= 0.);
134 
135  LevelDataType rhst, phit;
136  this->m_ops[a_level]->create(rhst, a_src);
137  this->m_ops[a_level]->create(phit, a_phiNew);
138 
139  this->m_ops[a_level]->setToZero(phit);
140  this->m_ops[a_level]->setToZero(rhst);
141  if (a_zeroPhi)
142  {
143  this->m_ops[a_level]->setToZero(a_phiNew);
144  }
145 
146  //set phit to a_phiOld, rhst to a_src * a_dt
147  this->m_ops[a_level]->incr(phit, a_phiOld, 1.0);
148  this->m_ops[a_level]->incr(rhst, a_src , a_dt);
149 
150  //multiply phi old by kappa*acoef
151  this->m_ops[a_level]->diagonalScale(phit, true);
152 
153  //multiply rhs by kappa (but NOT by a)
154  if (!a_rhsAlreadyKappaWeighted)
155  this->m_ops[a_level]->kappaScale(rhst);
156 
157  //add both together to make rhs for euler solve = kappa a phiold + kappa rhs dt
158  this->m_ops[a_level]->incr(rhst, phit, 1.0);
159 
160  //solve for phi new = (kappa a I - kappa dt L) phinew = kappa a phiold + kappa rhs
161  //this makes phinew = (k*a I - mu2 L)^-1 (rhs)
162  LevelDataType coarseData;
163  if ((a_crsePhiOldPtr != NULL) && (a_level > 0))
164  {
165  this->m_ops[a_level-1]->create(coarseData, *a_crsePhiOldPtr);
166  this->m_ops[a_level-1]->setToZero(coarseData);
167 
168  this->timeInterp(coarseData, *a_crsePhiOldPtr, *a_crsePhiNewPtr,
169  a_oldTime, 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 
215 
216 
217 private:
218 
219  // Disallowed operators.
223 };
224 
225 
226 #include "NamespaceFooter.H"
227 #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:60
A reference-counting handle class.
Definition: RefCountedPtr.H:66
BaseLevelBackwardEuler & operator=(const BaseLevelBackwardEuler &)
#define CH_assert(cond)
Definition: CHArray.H:37
A class to facilitate interaction with physical boundary conditions.
Definition: ProblemDomain.H:130
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:106
Definition: BaseLevelHeatSolver.H:44
Definition: BaseLevelBackwardEuler.H:45
Definition: DataIterator.H:140
const int SpaceDim
Definition: SPACE.H:39
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:511
Definition: BoxLayoutData.H:136
virtual ~BaseLevelBackwardEuler()
Destructor, called after destructors of BaseLevelBackwardEuler subclasses.
Definition: BaseLevelBackwardEuler.H:70
double Real
Definition: REAL.H:33
size_t size() const
Definition: Vector.H:177
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:484
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:387
Definition: AMRMultiGrid.H:231
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:339
Vector< LevelTGAHelmOp< LevelDataType, FluxDataType > *> m_ops
Definition: BaseLevelHeatSolver.H:522