Chombo + EB  3.2
BaseLevelCrankNicolson.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 _BASELEVELCRANKNICOLSON_H_
12 #define _BASELEVELCRANKNICOLSON_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 BaseLevelCrankNicolson
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 BaseLevelCrankNicolson : public BaseLevelHeatSolver<LevelDataType, FluxDataType, FluxRegisterType>
45 {
46 
47  public:
48 
49  //! Initializes the base class of a CrankNicolson time integrator. This must be called
50  //! by any subclass of BaseLevelCrankNicolson.
51  //! \param a_grids The DisjointBoxLayout on which the CrankNicolson 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 CrankNicolson 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 BaseLevelCrankNicolson 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  LevelDataType coarseData;
146  if ((a_crsePhiOldPtr != NULL) && (a_level > 0))
147  {
148  this->m_ops[a_level-1]->create(coarseData, *a_crsePhiOldPtr);
149  this->timeInterp(coarseData, *a_crsePhiOldPtr, *a_crsePhiNewPtr,
150  a_oldTime, a_crseOldTime, a_crseNewTime, a_level-1);
151  }
152 
153  //this makes phit = (k a I + dt/2 k L) phiOld
154  //inhomogeneous BC for solution
155  this->applyHelm(phit, a_phiOld, &coarseData, a_level, 0.5, a_dt, false);
156 
157  //set rhst = dt*src
158  this->m_ops[a_level]->incr(rhst, a_src , a_dt);
159 
160  //multiply phi old by kappa*acoef
161  this->m_ops[a_level]->diagonalScale(phit, true);
162 
163  //multiply rhs by kappa (but NOT by a) rhs = kappa dt rhs
164  if (!a_rhsAlreadyKappaWeighted)
165  this->m_ops[a_level]->kappaScale(rhst);
166 
167  //add both together to make rhs for crank solve = (kappa a I + kappa dt/2 L) phiold + kappa rhs dt
168  this->m_ops[a_level]->incr(rhst, phit, 1.0);
169 
170  this->solveHelm(a_phiNew, coarseData, rhst, a_level, 0.5, a_dt, a_zeroPhi);
171  this->incrementFlux(a_flux, a_phiNew, a_level, 0.5, a_dt, -1.0, false);
172 
173  // now increment flux registers -- note that because of the way
174  // we defined the fluxes, the dt multiplier is already in the
175  // flux
176  if ((a_fineFluxRegPtr != NULL) && (a_level < this->m_grids.size()-1))
177  {
178  Real fluxMult = 1.0;
179  for (DataIterator dit = this->m_grids[a_level].dataIterator(); dit.ok(); ++dit)
180  {
181  FluxDataType& thisFlux = a_flux[dit];
182  for (int dir=0; dir<SpaceDim; ++dir)
183  {
184  a_fineFluxRegPtr->incrementCoarse(thisFlux[dir],
185  fluxMult, dit(),
186  intervBase, // source
187  intervFlux, // dest
188  dir);
189  }
190  }
191  } // end if there is a finer-level
192 
193  if ((a_crseFluxRegPtr != NULL) && (a_level > 0))
194  {
195  Real fluxMult = 1.0;
196 
197  for (DataIterator dit = this->m_grids[a_level].dataIterator(); dit.ok(); ++dit)
198  {
199  FluxDataType& thisFlux = a_flux[dit];
200  for (int dir=0; dir<SpaceDim; ++dir)
201  {
202  a_crseFluxRegPtr->incrementFine(thisFlux[dir], fluxMult, dit(),
203  intervBase, // source
204  intervFlux, // dest
205  dir);
206  }
207  }
208  } // end if there is a coarser level
209  }
210 
211 private:
212  // Disallowed operators.
216 };
217 
218 #include "NamespaceFooter.H"
219 #endif
A reference-counting handle class.
Definition: RefCountedPtr.H:173
virtual ~BaseLevelCrankNicolson()
Destructor, called after destructors of BaseLevelCrankNicolson subclasses.
Definition: BaseLevelCrankNicolson.H:69
#define CH_assert(cond)
Definition: CHArray.H:37
A class to facilitate interaction with physical boundary conditions.
Definition: ProblemDomain.H:141
Definition: BaseLevelHeatSolver.H:43
BaseLevelCrankNicolson(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: BaseLevelCrankNicolson.H:59
Definition: DataIterator.H:190
void applyHelm(LevelDataType &a_ans, const LevelDataType &a_phi, const LevelDataType *a_phiC, int a_level, Real a_mu, Real a_dt, bool a_homogeneousBC)
Definition: BaseLevelHeatSolver.H:291
const int SpaceDim
Definition: SPACE.H:38
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: BaseLevelCrankNicolson.H:105
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
BaseLevelCrankNicolson & operator=(const BaseLevelCrankNicolson &)
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
Definition: BaseLevelCrankNicolson.H:44