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