Chombo + EB  3.0
BaseLevelTGA.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 _BASELEVELTGA_H__
12 #define _BASELEVELTGA_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 BaseLevelTGA
29 //! This base class implements the 2nd-order implicit L0-stable time
30 //! integration algorithm developed by Twizell, Gumel, and Arigu 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 BaseLevelTGA : public BaseLevelHeatSolver<LevelDataType, FluxDataType, FluxRegisterType>
46 {
47 
48  public:
49 
50  //! Initializes the base class of a TGA time integrator. This must be called
51  //! by any subclass of BaseLevelTGA.
52  //! \param a_grids The DisjointBoxLayout on which the TGA 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 TGA 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  m_mu1(0.0),
67  m_mu2(0.0),
68  m_mu3(0.0),
69  m_mu4(0.0),
70  m_r1(0.0)
71  {
72  Real tgaEpsilon = 1.e-12;
73 #ifdef CH_USE_FLOAT
74  tgaEpsilon = sqrt(tgaEpsilon);
75 #endif
76  Real a = 2.0 - sqrt(2.0) - tgaEpsilon;
77  m_mu1 = (a - sqrt( a*a - 4.0*a + 2.0))/2.0 ;
78  m_mu2 = (a + sqrt( a*a - 4.0*a + 2.0))/2.0 ;
79  m_mu3 = (1.0 - a);
80  m_mu4 = 0.5 - a;
81 
82  Real discr = sqrt(a*a - 4.0*a + 2.0);
83  m_r1 = (2.0*a - 1.0)/(a + discr);
84  }
85 
86  //! Destructor, called after destructors of BaseLevelTGA subclasses.
87  virtual ~BaseLevelTGA()
88  {
89  }
90 
91  //! Integrates the helmholtz equation represented by this object, placing
92  //! the new solution in \a a_phiNew.
93  //! \param a_phiNew The new solution (the value of phi at time n + 1) will
94  //! be stored here.
95  //! \param a_phiOld The old solution (the value of phi at time n).
96  //! \param a_src The source term on the right hand side of the Helmholtz
97  //! equation.
98  //! \param a_flux This will store the flux computed at the current grid
99  //! level during the solution of the Helmholtz equation.
100  //! \param a_fineFluxRegPtr A pointer to the flux register representing the
101  //! finer grid level adjacent to this one, or NULL
102  //! if there is no finer grid level.
103  //! \param a_crseFluxRegPtr A pointer to the flux register representing the
104  //! coarser grid level adjacent to this one, or NULL
105  //! if there is no coarser grid level.
106  //! \param a_oldTime The time at the beginning of the integration step at
107  //! the current grid level.
108  //! \param a_crseOldTime The time at the beginning of the integration step
109  //! at the coarser adjacent grid level. This parameter
110  //! is ignored if there is no coarser grid level.
111  //! \param a_crseNewTime The time at the end of the integration step
112  //! at the coarser adjacent grid level. This parameter
113  //! is ignored if there is no coarser grid level.
114  //! \param a_dt The size of the integration step at the current grid level.
115  //! \param a_level The current grid level.
116  //! \param a_zeroPhi If set to true, \a a_phiNew will be set to zero before
117  //! the integration takes place. Otherwise, a_phiNew is
118  //! assumed to be an initial estimate for the solution in
119  //! the iterative linear solve.
120  //! \param a_fluxStartComponent An index identifying the component at which
121  //! flux data begins within \a a_fineFluxRegPtr
122  //! and \a a_crseFluxRegPtr.
123  void updateSoln(LevelDataType& a_phiNew,
124  LevelDataType& a_phiOld,
125  LevelDataType& a_src,
126  LevelData<FluxDataType>& a_flux,
127  FluxRegisterType* a_fineFluxRegPtr,
128  FluxRegisterType* a_crseFluxRegPtr,
129  const LevelDataType* a_crsePhiOldPtr,
130  const LevelDataType* a_crsePhiNewPtr,
131  Real a_oldTime,
132  Real a_crseOldTime,
133  Real a_crseNewTime,
134  Real a_dt,
135  int a_level,
136  bool a_zeroPhi = true,
137  bool a_rhsAlreadyKappaWeighted = false,
138  int a_fluxStartComponent = 0)
139  {
140  // If our operators are time-independent, do the "easy" thing.
141  if (!this->m_ops[a_level]->isTimeDependent())
142  {
143  this->updateSolnWithTimeIndependentOp(a_phiNew, a_phiOld, a_src, a_flux,
144  a_fineFluxRegPtr, a_crseFluxRegPtr,
145  a_crsePhiOldPtr, a_crsePhiNewPtr,
146  a_oldTime, a_crseOldTime, a_crseNewTime,
147  a_dt, a_level, a_zeroPhi,
148  a_fluxStartComponent);
149  return;
150  }
151  else
152  {
153  // We have to assume that the operator and its coefficients are
154  // time-dependent.
155  this->updateSolnWithTimeDependentOp(a_phiNew, a_phiOld, a_src, a_flux,
156  a_fineFluxRegPtr, a_crseFluxRegPtr,
157  a_crsePhiOldPtr, a_crsePhiNewPtr,
158  a_oldTime, a_crseOldTime, a_crseNewTime,
159  a_dt, a_level, a_zeroPhi,
160  a_fluxStartComponent);
161  return;
162  }
163  }
164 
165  //! Computes the time-centered diffusion term L(phi). This can be used to
166  //! find contributions to the solution from diffusion. The diffusion term
167  //! is computed by computing a finite difference approximation for d phi/dt
168  //! using the updated and original values of phi and the time step. Most of
169  //! the arguments given here are passed along to updateSoln and retain their
170  //! significance therein.
171  //! \param a_diffusiveTerm The diffusion term L(phi) will be stored here.
172  //! \param a_phiOld The old solution (the value of phi at time n).
173  //! \param a_src The source term on the right hand side of the Helmholtz
174  //! equation (used to fine the new value of phi).
175  //! \param a_fineFluxRegPtr A pointer to the flux register representing the
176  //! finer grid level adjacent to this one, or NULL
177  //! if there is no finer grid level.
178  //! \param a_crseFluxRegPtr A pointer to the flux register representing the
179  //! coarser grid level adjacent to this one, or NULL
180  //! if there is no coarser grid level.
181  //! \param a_crsePhiOldTime A pointer to the value of phi at the beginning
182  //! of the integration step at the coarser grid
183  //! level adjacent to this one, or NULL if there
184  //! is no coarser grid level.
185  //! \param a_crsePhiNewTime A pointer to the value of phi at the end
186  //! of the integration step at the coarser grid
187  //! level adjacent to this one, or NULL if there
188  //! is no coarser grid level.
189  //! \param a_oldTime The time at the beginning of the integration step at
190  //! the current grid level.
191  //! \param a_crseOldTime The time at the beginning of the integration step
192  //! at the coarser adjacent grid level. This parameter
193  //! is ignored if there is no coarser grid level.
194  //! \param a_crseNewTime The time at the end of the integration step
195  //! at the coarser adjacent grid level. This parameter
196  //! is ignored if there is no coarser grid level.
197  //! \param a_dt The size of the integration step at the current grid level.
198  //! \param a_level The current grid level.
199  //! \param a_zeroPhi If set to true, the new value of phi will be set to
200  //! zero before the integration takes place. Otherwise, it
201  //! will be set to the value in \a a_diffusiveTerm.
202  void computeDiffusion(LevelDataType& a_diffusiveTerm,
203  LevelDataType& a_phiOld,
204  LevelDataType& a_src,
205  LevelData<FluxDataType>& a_flux,
206  FluxRegisterType* a_fineFluxRegPtr,
207  FluxRegisterType* a_crseFluxRegPtr,
208  const LevelDataType* a_crsePhiOldPtr,
209  const LevelDataType* a_crsePhiNewPtr,
210  Real a_oldTime,
211  Real a_crseOldTime,
212  Real a_crseNewTime,
213  Real a_dt,
214  int a_level,
215  bool a_zeroPhi = true,
216  bool a_rhsAlreadyKappaWeighted = false)
217  {
218  // in this function, we call the updateSoln function, compute
219  // the update, then subtract off the extra pieces to return the
220  // diffusive part of the update
221 
222  if (!this->m_ops[a_level]->isTimeDependent())
223  {
224  // The operator has no time-dependent parameters. Life is easier.
225 
226  // first compute updated solution
227  LevelDataType phiNew;
228 
229  this->m_ops[a_level]->create(phiNew, a_phiOld);
230  this->m_ops[a_level]->setToZero(phiNew);
231  if (!a_zeroPhi)
232  {
233  this->m_ops[a_level]->assign(phiNew, a_phiOld);
234  }
235 
236  updateSoln(phiNew, a_phiOld, a_src, a_flux,
237  a_fineFluxRegPtr, a_crseFluxRegPtr,
238  a_crsePhiOldPtr, a_crsePhiNewPtr,
239  a_oldTime, a_crseOldTime,
240  a_crseNewTime, a_dt, a_level, a_zeroPhi, a_rhsAlreadyKappaWeighted);
241 
242  // now subtract everything off to leave us with diffusive term
243  this->m_ops[a_level]->incr(phiNew, a_phiOld, -1.0);
244  this->m_ops[a_level]->scale(phiNew, 1.0/a_dt);
245 
246  //now multiply by a if there is an a
247  this->m_ops[a_level]->diagonalScale(phiNew, false);
248 
249  // and finally, subtract off a_src
250  this->m_ops[a_level]->incr(phiNew, a_src, -1.0);
251 
252  // what's left should be the time-centered diffusive part of the update
253  this->m_ops[a_level]->assign(a_diffusiveTerm, phiNew);
254  }
255  else
256  {
257  // The operator has time-dependent coefficients. We must be more careful!
258 
259  // n+1 n
260  // n (a - a )
261  // There's an extra source term (phi (da/dt) = phi * -----------
262  // dt
263  // from the time-changing density that we need to subtract from the RHS.
264  LevelDataType phidadt, aOld, aNew, rhs;
265  this->m_ops[a_level]->create(phidadt, a_phiOld);
266  this->m_ops[a_level]->create(aOld, a_phiOld);
267  this->m_ops[a_level]->create(aNew, a_phiOld);
268  this->m_ops[a_level]->create(rhs, a_phiOld);
269  for (DataIterator dit = a_phiOld.disjointBoxLayout().dataIterator(); dit.ok(); ++dit)
270  {
271  // Set the old and new a coefficients.
272  this->m_ops[a_level]->setTime(a_oldTime, 0.0, a_dt);
273  aOld[dit()].setVal(1.);
274  this->m_ops[a_level]->diagonalScale(aOld);
275  this->m_ops[a_level]->setTime(a_oldTime, 1.0, a_dt);
276  aNew[dit()].setVal(1.);
277  this->m_ops[a_level]->diagonalScale(aNew);
278 
279  // Compute the above term.
280  phidadt[dit()].axby(aNew[dit()], aOld[dit()], 1.0/a_dt, -1.0/a_dt);
281  phidadt[dit()] *= a_phiOld[dit()];
282 
283  // Make a new right hand side out of the difference between the
284  // source term and phidadt.
285  rhs[dit()].axby(a_src[dit()], phidadt[dit()], 1.0, -1.0);
286  }
287 
288  // Compute the updated solution.
289  LevelDataType phiNew;
290  this->m_ops[a_level]->create(phiNew, a_phiOld);
291  this->m_ops[a_level]->setToZero(phiNew);
292  if (!a_zeroPhi)
293  {
294  this->m_ops[a_level]->assign(phiNew, a_phiOld);
295  }
296  updateSoln(phiNew, a_phiOld, rhs, a_flux,
297  a_fineFluxRegPtr, a_crseFluxRegPtr,
298  a_crsePhiOldPtr, a_crsePhiNewPtr,
299  a_oldTime, a_crseOldTime,
300  a_crseNewTime, a_dt, a_level, a_zeroPhi);
301 
302  // n+1 n
303  // ([a phi] - [a phi] )
304  // ---------------------- - a_src -> a_diffusiveTerm.
305  // dt
306  for (DataIterator dit = a_phiOld.disjointBoxLayout().dataIterator(); dit.ok(); ++dit)
307  {
308  aNew[dit()] *= phiNew[dit()];
309  aOld[dit()] *= a_phiOld[dit()];
310  a_diffusiveTerm[dit()].axby(aNew[dit()], aOld[dit()], 1.0/a_dt, -1.0/a_dt);
311  a_diffusiveTerm[dit()] -= a_src[dit()];
312  }
313  }
314  }
315 
316 
317 
318 
319  protected:
320 
321  //! Sets the value of the source term in the Helmholtz equation on ghost
322  //! cells. This method should be overridden by subclasses of BaseLevelTGA.
323  //! \param a_src The source term in the Helmholtz equation whose ghost cell
324  //! values are to be set.
325  //! \param a_dbl The disjoint box layout that indexes \a a_src.
326  //! \param a_level The grid level at which the ghost cell values are to be
327  //! set.
328  virtual void setSourceGhostCells(LevelDataType& a_src,
329  const DisjointBoxLayout& a_dbl,
330  int a_level) = 0;
331 
332  //! The times within the integration step at which the operators are
333  //! evaluated.
335 
336 private:
337 
338  // Update the solution assuming that the operator's coefficients are
339  // independent of time. Same arguments as updateSoln.
340  void updateSolnWithTimeIndependentOp(LevelDataType& a_phiNew,
341  LevelDataType& a_phiOld,
342  LevelDataType& a_src,
343  LevelData<FluxDataType>& a_flux,
344  FluxRegisterType* a_fineFluxRegPtr,
345  FluxRegisterType* a_crseFluxRegPtr,
346  const LevelDataType* a_crsePhiOldPtr,
347  const LevelDataType* a_crsePhiNewPtr,
348  Real a_oldTime,
349  Real a_crseOldTime,
350  Real a_crseNewTime,
351  Real a_dt,
352  int a_level,
353  bool a_zeroPhi = true,
354  int a_fluxStartComponent = 0)
355  {
356 
357  CH_assert(!this->m_ops[a_level]->isTimeDependent());
358  int ncomp = a_phiNew.nComp();
359  Interval intervBase(0, ncomp-1);
360  Interval intervFlux(a_fluxStartComponent, a_fluxStartComponent + ncomp-1);
361 
362  CH_assert(a_level >= 0);
363  CH_assert(a_level < this->m_grids.size());
364  CH_assert((a_level == 0) || (a_crsePhiOldPtr != NULL));
365  CH_assert((a_level == 0) || (a_crsePhiNewPtr != NULL));
366  CH_assert(a_crseNewTime >= a_crseOldTime);
367  CH_assert(a_dt >= 0.);
368 
369  LevelDataType rhst, srct, phis;
370  this->m_ops[a_level]->create(rhst, a_src);
371  this->m_ops[a_level]->create(srct, a_phiNew);
372  this->m_ops[a_level]->create(phis, a_phiNew);
373 
374  this->m_ops[a_level]->setToZero(srct);
375  this->m_ops[a_level]->setToZero(rhst);
376  //this makes srct = a_src*dt
377  this->m_ops[a_level]->incr(srct, a_src, a_dt);
378 
379  //save input guess if we are not using zero
380  if (!a_zeroPhi)
381  {
382  this->m_ops[a_level]->setToZero(phis);
383  this->m_ops[a_level]->incr(phis, a_phiNew, 1.0);
384  }
385 
386  // Divide the source S by the identity coefficient a. srct = rhs*dt/a
387  this->m_ops[a_level]->divideByIdentityCoef(srct);
388 
389  LevelDataType coarseData;
390 
391  if ((a_crsePhiOldPtr != NULL) && (a_level > 0))
392  {
393  this->m_ops[a_level-1]->create(coarseData, *a_crsePhiOldPtr);
394  setSourceGhostCells(srct, this->m_grids[a_level], a_level);
395  }
396 
397  //from here on k is kappa and L is kappa L
398  //this makes rhs hold (k*a I + mu4 k L) (S/a) dt
399  this->applyHelm(rhst, srct, NULL, a_level, m_mu4, a_dt, true); //true means the application is homogeneous
400  this->incrementFlux(a_flux, srct, a_level, m_mu4, a_dt, -1.0, true);
401 
402  // first need to compute coarse-level BC at this level's old time
403  if (a_level > 0)
404  {
405  this->timeInterp(coarseData, *a_crsePhiOldPtr, *a_crsePhiNewPtr,
406  a_oldTime, a_crseOldTime, a_crseNewTime, a_level-1);
407  }
408 
409  //this makes a_phiNew hold (k*a I + mu3 k L) phi^n
410  //'true' apply CF and domain BC
411  this->applyHelm(a_phiNew, a_phiOld, &coarseData, a_level, m_mu3, a_dt, true);
412  this->incrementFlux(a_flux, a_phiOld, a_level, m_mu3, a_dt, -1., false);
413 
414  //this makes rhs hold (k*a I + mu3 L) phi^n + dt(k*a I +mu4 L) S/a
415  this->m_ops[a_level]->incr(rhst, a_phiNew, 1.0);
416 
417  // now construct coarse-level BC for intermediate solve
418  // intermediate solution will be at time = oldTime + (1-r1)dt
419  if (a_level > 0)
420  {
421  Real intermediateTime = a_oldTime + (1.0-m_r1)*a_dt;
422 
423  this->timeInterp(coarseData, *a_crsePhiOldPtr, *a_crsePhiNewPtr,
424  intermediateTime, a_crseOldTime, a_crseNewTime, a_level-1);
425  }
426 
427  //if user thought her original guess was good, use it.
428  if (!a_zeroPhi)
429  {
430  this->m_ops[a_level]->setToZero(a_phiNew);
431  this->m_ops[a_level]->incr(a_phiNew, phis, 1.0);
432  }
433 
434  //by now rhs = ((k*a I + mu3 L) phi^n + dt(k*a I +mu4 L) S/a )
435  //this makes phinew = (k*a I - mu2 k L)^-1 (rhs)
436  this->solveHelm(a_phiNew, coarseData, rhst, a_level, m_mu2, a_dt, a_zeroPhi);
437  this->incrementFlux(a_flux, a_phiNew, a_level, m_mu2, a_dt, -1.0, false);
438 
439  //this puts the answer into rhst so we can do the final solve
440  this->m_ops[a_level]->assign(rhst, a_phiNew);
441 
442 
443  // now construct CF-BC for final solve
444  if (a_level > 0)
445  {
446  Real newTime = a_oldTime + a_dt;
447  this->timeInterp(coarseData, *a_crsePhiOldPtr, *a_crsePhiNewPtr,
448  newTime, a_crseOldTime, a_crseNewTime, a_level-1);
449  }
450 
451  //by now rhs = ((k*a I + mu3 L) phi^n + dt(k*a I +mu4 L) S/a )
452  //this makes rhs hold k*a[(k*a I - mu2 L)^-1 (rhs)]
453  this->m_ops[a_level]->diagonalScale(rhst, true);
454 
455  //if user thought her original guess was good, use it again.
456  if (!a_zeroPhi)
457  {
458  this->m_ops[a_level]->setToZero(a_phiNew);
459  this->m_ops[a_level]->incr(a_phiNew, phis, 1.0);
460  }
461 
462  //this makes phinew = (k*a I - mu1 L)^-1 [ka ((k*a I - mu2 L)^-1 (orig rhs))]
463  this->solveHelm(a_phiNew, coarseData, rhst, a_level, m_mu1, a_dt, a_zeroPhi);
464  this->incrementFlux(a_flux, a_phiNew, a_level, m_mu1, a_dt, -1.0, false);
465 
466  // now increment flux registers -- note that because of the way
467  // we defined the fluxes, the dt multiplier is already in the
468  // flux
469  if ((a_fineFluxRegPtr != NULL) && (a_level < this->m_grids.size()-1))
470  {
471  Real fluxMult = 1.0;
472  for (DataIterator dit = this->m_grids[a_level].dataIterator(); dit.ok(); ++dit)
473  {
474  FluxDataType& thisFlux = a_flux[dit];
475  for (int dir=0; dir<SpaceDim; ++dir)
476  {
477  a_fineFluxRegPtr->incrementCoarse(thisFlux[dir],
478  fluxMult, dit(),
479  intervBase, // source
480  intervFlux, // dest
481  dir);
482  }
483  }
484  } // end if there is a finer-level
485 
486  if ((a_crseFluxRegPtr != NULL) && (a_level > 0))
487  {
488  Real fluxMult = 1.0;
489 
490  for (DataIterator dit = this->m_grids[a_level].dataIterator(); dit.ok(); ++dit)
491  {
492  FluxDataType& thisFlux = a_flux[dit];
493  for (int dir=0; dir<SpaceDim; ++dir)
494  {
495  a_crseFluxRegPtr->incrementFine(thisFlux[dir], fluxMult, dit(),
496  intervBase, // source
497  intervFlux, // dest
498  dir);
499  }
500  }
501  } // end if there is a coarser level
502 
503  }
504 
505 
506  // Update the solution assuming that the operator's coefficients change
507  // with time. Same arguments as updateSoln.
508  void updateSolnWithTimeDependentOp(LevelDataType& a_phiNew,
509  LevelDataType& a_phiOld,
510  LevelDataType& a_src,
511  LevelData<FluxDataType>& a_flux,
512  FluxRegisterType* a_fineFluxRegPtr,
513  FluxRegisterType* a_crseFluxRegPtr,
514  const LevelDataType* a_crsePhiOldPtr,
515  const LevelDataType* a_crsePhiNewPtr,
516  Real a_oldTime,
517  Real a_crseOldTime,
518  Real a_crseNewTime,
519  Real a_dt,
520  int a_level,
521  bool a_zeroPhi = true,
522  int a_fluxStartComponent = 0)
523  {
524  int ncomp = a_phiNew.nComp();
525  Interval intervBase(0, ncomp-1);
526  Interval intervFlux(a_fluxStartComponent, a_fluxStartComponent + ncomp-1);
527 
528  CH_assert(a_level >= 0);
529  CH_assert(a_level < this->m_grids.size());
530  CH_assert((a_level == 0) || (a_crsePhiOldPtr != NULL));
531  CH_assert((a_level == 0) || (a_crsePhiNewPtr != NULL));
532  CH_assert(a_crseNewTime >= a_crseOldTime);
533  CH_assert(a_dt >= 0.);
534 
535  LevelDataType rhst, srct, phis;
536  this->m_ops[a_level]->create(rhst, a_src);
537  this->m_ops[a_level]->create(srct, a_phiNew);
538  this->m_ops[a_level]->create(phis, a_phiNew);
539 
540  this->m_ops[a_level]->setToZero(srct);
541  this->m_ops[a_level]->setToZero(rhst);
542  this->m_ops[a_level]->incr(srct, a_src, 1.0);
543 
544  //save input guess if we are not using zero
545  if (!a_zeroPhi)
546  {
547  this->m_ops[a_level]->setToZero(phis);
548  this->m_ops[a_level]->incr(phis, a_phiNew, 1.0);
549  }
550 
551  // In the comments, we use superscripts to denote different time centerings
552  // in the operator L, the identity coefficient a, and the solution phi.
553  // A '0' superscript means the quantity at time n.
554  // A '1/2' superscript means the quantity at time n + 1/2.
555  // A '*' superscript means the quantity at the "intermediate" time used
556  // by the TGA scheme.
557  // A '1' superscript means the quantity at time n + 1.
558  //
559  // In embedded boundary problems we use k for kappa and L means kappa L.
560 
561  //-----------------------------------------------------
562  // 1/2
563  // dt 1/2 1/2 S
564  // Compute ---- x (a I + mu4 L ) x --- -> rhst
565  // 1/2 1/2
566  // a a
567  //-----------------------------------------------------
568 
569  // First, feed the half-step time to the operator to set its coefficients.
570  this->m_ops[a_level]->setTime(a_oldTime, 0.5, a_dt);
571 
572  // 1/2
573  // Divide the source S by the identity coefficient a .
574  this->m_ops[a_level]->divideByIdentityCoef(srct);
575 
576  // Set the ghost cells in the source at the coarser level.
577  LevelDataType coarseData;
578  if ((a_crsePhiOldPtr != NULL) && (a_level > 0))
579  {
580  this->m_ops[a_level-1]->create(coarseData, *a_crsePhiOldPtr);
581  setSourceGhostCells(srct, this->m_grids[a_level], a_level);
582  }
583 
584  // 1/2 1/2
585  // (k a I + mu4 k L) (S / a ) -> rhst
586  // (The false parameter tells the operator to use the source's ghost data
587  // instead of applying boundary conditions).
588  //homogeneous BC for L(rhs)
589  this->applyHelm(rhst, srct, NULL, a_level, m_mu4, a_dt, true);
590  this->incrementFlux(a_flux, srct, a_level, a_dt*m_mu4, a_dt, -1.0, true);
591 
592  //--------------------------------------------------------
593  // 1 0 0 0
594  // Compute --- x (k a I + mu4 k L ) x phi and add it to rhst.
595  // 0
596  // a
597  //--------------------------------------------------------
598 
599  // First we compute the coarse-level boundary condition data at time n.
600  if (a_level > 0)
601  {
602  this->timeInterp(coarseData, *a_crsePhiOldPtr, *a_crsePhiNewPtr,
603  a_oldTime, a_crseOldTime, a_crseNewTime, a_level-1);
604  }
605 
606  // Feed the old time to the operator.
607  this->m_ops[a_level]->setTime(a_oldTime, 0.0, a_dt);
608 
609  // 0 0 0
610  // Compute (k a I + mu3 L ) phi and store it in a_phiNew for now.
611  // (The true parameter means apply the coarse-fine and domain boundary
612  // conditions).
613  //inhommogeneous bc for solution
614  this->applyHelm(a_phiNew, a_phiOld, &coarseData, a_level, m_mu3, a_dt, false);
615  this->incrementFlux(a_flux, a_phiOld, a_level, m_mu3, a_dt, -1., false);
616 
617  // n
618  // Divide a_phiNew by a .
619  this->m_ops[a_level]->divideByIdentityCoef(a_phiNew);
620 
621  // Add a_phiNew to rhst to obtain the right hand side of our
622  // TGA update equation.
623  this->m_ops[a_level]->incr(rhst, a_phiNew, 1.0);
624 
625  //-----------------------------------------------------------------
626  // * * n *
627  // Solve (k a I + mu4 k L ) phi = a rhst and place the answer in rhst,
628  //-----------------------------------------------------------------
629 
630  // Construct the coarse-level boundary condition data for our TGA
631  // * n
632  // intermediate solve. The intermediate time is t = t + (1-r1)dt.
633  if (a_level > 0)
634  {
635  Real intermediateTime = a_oldTime + (1.0-m_r1)*a_dt;
636 
637  this->timeInterp(coarseData, *a_crsePhiOldPtr, *a_crsePhiNewPtr,
638  intermediateTime, a_crseOldTime, a_crseNewTime, a_level-1);
639  }
640 
641  // If the user thought her original guess was good, use it.
642  if (!a_zeroPhi)
643  {
644  this->m_ops[a_level]->setToZero(a_phiNew);
645  this->m_ops[a_level]->incr(a_phiNew, phis, 1.0);
646  }
647 
648  // Feed the intermediate time to the operator to get the time-centered
649  // coefficients.
650  this->m_ops[a_level]->setTime(a_oldTime, 1.0 - m_r1, a_dt);
651 
652  // *
653  // Multiply rhst by a .
654  this->m_ops[a_level]->diagonalScale(rhst, false);
655 
656  // * * -1
657  // This computes (k a I - mu2 L ) (rhst) and stores it in a_phiNew.
658  this->solveHelm(a_phiNew, coarseData, rhst, a_level, m_mu2, a_dt, a_zeroPhi);
659  this->incrementFlux(a_flux, a_phiNew, a_level, m_mu2, a_dt, -1.0, false);
660 
661  // This puts the answer into rhst.
662  this->m_ops[a_level]->assign(rhst, a_phiNew);
663 
664  //--------------------------------------------------
665  // 1 1 1 1
666  // Solve (k a I + k mu4 L ) phi = rhst and scale by a .
667  //--------------------------------------------------
668 
669  // Compute the coarse-fine boundary condition data for the final solve.
670  if (a_level > 0)
671  {
672  Real newTime = a_oldTime + a_dt;
673  this->timeInterp(coarseData, *a_crsePhiOldPtr, *a_crsePhiNewPtr,
674  newTime, a_crseOldTime, a_crseNewTime, a_level-1);
675  }
676 
677  // Feed the new time to the operator.
678  this->m_ops[a_level]->setTime(a_oldTime, 1.0, a_dt);
679 
680  // 1
681  // Scale rhst by k a .
682  this->m_ops[a_level]->diagonalScale(rhst);
683 
684  // If the user thought her original guess was good, use it again.
685  if (!a_zeroPhi)
686  {
687  this->m_ops[a_level]->setToZero(a_phiNew);
688  this->m_ops[a_level]->incr(a_phiNew, phis, 1.0);
689  }
690 
691  // 1 1 -1 1 * * -1 *
692  // Compute [k a I - mu1 k L )] k a [(k a I - mu2 k L ) a (orig rhs)].
693  this->solveHelm(a_phiNew, coarseData, rhst, a_level, m_mu1, a_dt, a_zeroPhi);
694  this->incrementFlux(a_flux, a_phiNew, a_level, m_mu1, a_dt, -1.0, false);
695 
696  // Now increment the flux registers -- note that because of the way
697  // we defined the fluxes, the dt multiplier is already in the
698  // flux.
699  if ((a_fineFluxRegPtr != NULL) && (a_level < this->m_grids.size()-1))
700  {
701  Real fluxMult = 1.0;
702  for (DataIterator dit = this->m_grids[a_level].dataIterator(); dit.ok(); ++dit)
703  {
704  FluxDataType& thisFlux = a_flux[dit];
705  for (int dir=0; dir<SpaceDim; ++dir)
706  {
707  a_fineFluxRegPtr->incrementCoarse(thisFlux[dir],
708  fluxMult, dit(),
709  intervBase, // source
710  intervFlux, // dest
711  dir);
712  }
713  }
714  } // end if there is a finer-level
715 
716  if ((a_crseFluxRegPtr != NULL) && (a_level > 0))
717  {
718  Real fluxMult = 1.0;
719 
720  for (DataIterator dit = this->m_grids[a_level].dataIterator(); dit.ok(); ++dit)
721  {
722  FluxDataType& thisFlux = a_flux[dit];
723  for (int dir=0; dir<SpaceDim; ++dir)
724  {
725  a_crseFluxRegPtr->incrementFine(thisFlux[dir], fluxMult, dit(),
726  intervBase, // source
727  intervFlux, // dest
728  dir);
729  }
730  }
731  } // end if there is a coarser level
732  }
733 
734 
735 private:
736 
737  // Disallowed operators.
739  BaseLevelTGA(const BaseLevelTGA& a_opin);
740  BaseLevelTGA();
741 };
742 
743 
744 #include "NamespaceFooter.H"
745 #endif
Real m_mu1
Definition: BaseLevelTGA.H:334
A reference-counting handle class.
Definition: RefCountedPtr.H:66
#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
void updateSolnWithTimeIndependentOp(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, int a_fluxStartComponent=0)
Definition: BaseLevelTGA.H:340
virtual void setSourceGhostCells(LevelDataType &a_src, const DisjointBoxLayout &a_dbl, int a_level)=0
void computeDiffusion(LevelDataType &a_diffusiveTerm, 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)
Definition: BaseLevelTGA.H:202
void updateSolnWithTimeDependentOp(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, int a_fluxStartComponent=0)
Definition: BaseLevelTGA.H:508
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: BaseLevelTGA.H:123
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
Vector< T > & assign(const T &inval)
assign a scalar to every element of the vector
Definition: Vector.H:138
Real m_mu2
Definition: BaseLevelTGA.H:334
double Real
Definition: REAL.H:33
Definition: BaseLevelTGA.H:45
A BoxLayout that has a concept of disjointedness.
Definition: DisjointBoxLayout.H:31
Real m_mu4
Definition: BaseLevelTGA.H:334
virtual ~BaseLevelTGA()
Destructor, called after destructors of BaseLevelTGA subclasses.
Definition: BaseLevelTGA.H:87
BaseLevelTGA(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: BaseLevelTGA.H:60
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
virtual bool ok() const
return true if this iterator is still in its Layout
Definition: LayoutIterator.H:110
Real m_mu3
Definition: BaseLevelTGA.H:334
BaseLevelTGA & operator=(const BaseLevelTGA &)
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
Real m_r1
Definition: BaseLevelTGA.H:334