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