Chombo + EB  3.0
AMRTGA.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 _AMRTGA_H_
12 #define _AMRTGA_H_
13 
14 #include <cmath>
15 #include <cstdlib>
16 #include <cstdio>
17 #include <iostream>
18 #include <iomanip>
19 #include <fstream>
20 #include <string>
21 #include "AMRMultiGrid.H"
22 #include "NamespaceHeader.H"
23 
24 //! \class TGAHelmOp
25 //! This operator is meant to represent the general helmholtz operator that
26 //! AMRTGA will be solving. This operator is of the form
27 //! alpha A(x) phi + beta div B(x) grad phi = rho.
28 //! AMRTGA needs to reset the constants alpha and beta often.
29 //! \tparam T The LevelData class that holds solution data for the operator.
30 template <class T>
31 class TGAHelmOp : public AMRLevelOp<T>
32 {
33  public:
34 
35  //! Base class default constructor. This constructs a time-independent
36  //! operator.
38  :m_isTimeDependent(false),
40  {
41  }
42 
43  //! Base class constructor which specifies explicitly whether this
44  //! operator is time-dependent.
45  //! \param a_isTimeDependent If set to true, this Helmholtz operator will
46  //! be treated as a time-dependent operator.
47  explicit TGAHelmOp(bool a_isTimeDependent)
48  :m_isTimeDependent(a_isTimeDependent),
50  {
51  }
52 
53  //! Destructor.
54  virtual ~TGAHelmOp()
55  {
56  }
57 
58  //! Sets the scaling constants in the Helmholtz equation.
59  //! \param a_alpha The scaling constant that multiplies the identity term.
60  //! \param a_beta The scaling constant that multiplies the derivative term.
61  virtual void setAlphaAndBeta(const Real& a_alpha,
62  const Real& a_beta) = 0;
63 
64  //! Scales the right hand side of the Helmholtz equation by the
65  //! identity term in the operator. If you are solving
66  //! rho(x) dphi/dt = L(phi), this means multiply by rho (or kappa * rho in
67  //! embedded boundary calculations.
68  //! \param a_rhs The right hand side of the equation to be scaled.
69  //! \param a_kappaWeighted If set to true, \a a_rhs will be scaled by the
70  //! volume fraction in addition to the identity term.
71  virtual void diagonalScale(T& a_rhs, bool a_kappaWeighted)
72  {
73  MayDay::Error("Two argument version of diagonalScale called but not implemented!");
74  }
75 
76  /// for eb only. kappa weight the rhs but do not multiply by the identity coefficient
77  virtual void kappaScale(T& a_rhs)
78  {
79  // default implementation does nothing because only relevant in eb.
80  }
81 
82  //! Scales the right hand side of the Helmholtz equation by the
83  //! identity term in the operator. This version assumes volume fraction
84  //! weighting.
85  //! \param a_rhs The right hand side of the equation to be scaled.
86  virtual void diagonalScale(T& a_rhs)
87  {
88  diagonalScale(a_rhs, true);
89  }
90 
91  //! Divides the right hand side of the Helmholtz equation by the
92  //! identity coefficient rho(x) in the equation rho(x) dphi/dt = L(phi).
93  //! \param a_rhs The right hand side of the equation to be scaled.
94  virtual void divideByIdentityCoef(T& a_rhs) = 0;
95 
96  //! Apply the differential operator without any boundary or coarse-fine
97  //! boundary conditions and no finer level
98  //! \param a_ans The result of the application of the operator to \a a_phi.
99  //! \param a_phi The data to which the operator is applied.
100  virtual void applyOpNoBoundary(T& a_ans, const T& a_phi) = 0;
101 
102  //! Sets the time-dependent state of the operator. The default implementation
103  //! does nothing and is appropriate for time-independent operators.
104  //! \param a_time The time to be used to update the time-dependent operator.
105  virtual void setTime(Real a_time)
106  {
107  ;// most operators are time-independent.
108  }
109 
110  //! Sets the time-dependent state of the operator. This version of setTime
111  //! allows one to linearly interpolate coefficients across an integration
112  //! step, since it accepts arguments that define where in the step it is
113  //! to be updated. The default implementation calls
114  //! setTime(\a a_oldTime + a_mu * a_dt).
115  //! \param a_oldTime The time at the beginning of the current step.
116  //! \param a_mu The fraction of the current step that has elapsed.
117  //! \param a_dt The size of the current step.
118  virtual void setTime(Real a_oldTime, Real a_mu, Real a_dt)
119  {
120  setTime(a_oldTime + a_mu * a_dt);
121  }
122 
123  //! Returns true if the operator is time-dependent, false otherwise.
124  bool isTimeDependent() const
125  {
126  return m_isTimeDependent;
127  }
128 
129  //! Allows access to the identity coefficient data for the operator.
130  virtual T& identityCoef()
131  {
132  MayDay::Error("Access to identity coefficient not allowed for this operator");
133  // Yes, this is a dereference to a null pointer. No, it shouldn't be
134  // called. I'm putting this here to make the compiler shaddup. -JNJ
135  return *m_phonyIdentityCoef;
136  }
137 
138  private:
139 
140  // Disallowed operations.
141  TGAHelmOp(const TGAHelmOp&);
142  TGAHelmOp& operator=(const TGAHelmOp&);
143 
144  //! This flag is set if the Helmholtz operator's coefficiens are
145  //! time-dependent.
147 
148  // Phony identity coefficient data to avoid compiler warnings. This
149  // should go away soon.
151 };
152 
153 //! \class LevelTGAHelmOp
154 //! This subclass of TGAHelmOp proves additional functionality for computing
155 //! fluxes at refinement level boundaries.
156 //! \tparam T The LevelData class that holds solution data.
157 //! \tparam TFlux The FArrayBox class that holds flux data within a LevelData.
158 template <class T, class TFlux>
159 class LevelTGAHelmOp : public TGAHelmOp< T >
160 {
161  public:
162 
163  //! Base class default constructor. This constructs a time-independent
164  //! LevelTGAHelmOp.
166  :TGAHelmOp<T>::TGAHelmOp(false)
167  {
168  }
169 
170  //! Base class constructor which specifies explicitly whether this
171  //! operator is time-dependent.
172  //! \param a_isTimeDependent If set to true, this Helmholtz operator will
173  //! be treated as a time-dependent operator.
174  LevelTGAHelmOp(bool a_isTimeIndependent)
175  :TGAHelmOp<T>::TGAHelmOp(a_isTimeIndependent)
176  {
177  }
178 
179  //! Destructor.
180  virtual ~LevelTGAHelmOp()
181  {
182  }
183 
184  /// These functions are part of the LevelTGA interface......
185 
186  ///
187  virtual void fillGrad(const T& a_phi) = 0;
188 
189  ///
190  virtual void getFlux(TFlux& a_flux,
191  const T& a_data,
192  const Box& a_grid,
193  const DataIndex& a_dit,
194  Real a_scale) = 0;
195 
196 };
197 
198 
199 
200 ///
201 /**
202  Template implementation of the TGA algorithm
203  to solve the heat equation.
204  **/
205 template <class T>
206 class AMRTGA
207 {
208 public:
209 
210  ///
211  TGAHelmOp<T>* newOp(const ProblemDomain& a_indexSpace,
212  const AMRLevelOpFactory<T>& a_opFact)
213  {
214  AMRLevelOpFactory<T>& opFact = (AMRLevelOpFactory<T>&) a_opFact;
215  TGAHelmOp<T>* retval = (TGAHelmOp<T>*) opFact.AMRnewOp(a_indexSpace);
216  return retval;
217  }
218 
219  ///
220  ~AMRTGA();
221 
222  ///
223  /**
224  **/
225  AMRTGA(const RefCountedPtr<AMRMultiGrid<T> > & a_solver,
226  const AMRLevelOpFactory<T>& a_factory,
227  const ProblemDomain& a_level0Domain,
228  const Vector<int>& a_refRat,
229  int a_numLevels = -1,
230  int a_verbosity = 0);
231 
232 
233  ///
234  /**
235  This advances a parabolic pde from a_phiOld to a_phiNew using TGA
236  on a non-moving domain with source term a_source
237  If your operator is time-dependent, be sure to send the time.
238  **/
239  void oneStep(Vector<T*>& a_phiNew,
240  Vector<T*>& a_phiOld,
241  Vector<T*>& a_source,
242  const Real& a_dt,
243  int a_lbase,
244  int a_lmax,
245  Real a_timeOld = 0);
246 
247  ///
248  void resetAlphaAndBeta(const Real& a_alpha,
249  const Real& a_beta);
250 
251  void setTime(Real time);
252 
253 protected:
254  void solveHelm(Vector<T*>& a_ans,
255  Vector<T*>& a_rhs,
256  int a_lbase,
257  int a_lmax,
258  Real a_mu,
259  Real a_dt);
260 
261  void applyHelm(Vector<T*>& a_ans,
262  Vector<T*>& a_source,
263  int a_lbase,
264  int a_lmax,
265  Real a_mu,
266  Real a_dt,
267  bool a_homogeneousBC);
268 
269  void setMu();
270 
271  void createData(Vector<T* >& a_source,
272  int a_lbase,
273  int a_lmax);
274 
275 private:
276 
277 
278  //You do not own these operators!! don't delete it. the memory is
279  //owned by the solver
289 
290  //copy constructor and operator= disallowed for all the usual reasons
291  AMRTGA(const AMRTGA<T>& a_opin)
292  {
293  MayDay::Error("invalid operator");
294  }
295 
296  void operator=(const AMRTGA<T>& a_opin)
297  {
298  MayDay::Error("invalid operator");
299  }
300 
301  /// weak construction is bad. Ref Counted pointers are your friends.
303  {
304  MayDay::Error("invalid operator");
305  }
306 
307 
308 };
309 
310 /*****/
311 template <class T>
312 void AMRTGA<T>::
313 resetAlphaAndBeta(const Real& a_alpha,
314  const Real& a_beta)
315 {
316  Vector<MGLevelOp<T>* > ops = m_solver->getAllOperators();
317  for (int iop = 0; iop < ops.size(); iop++)
318  {
319 
320  TGAHelmOp<T>* helmop = (TGAHelmOp<T>*) ops[iop];
321  helmop->setAlphaAndBeta(a_alpha, a_beta);
322  }
323 }
324 
325 /*****/
326 template <class T>
327 void AMRTGA<T>::
328 setTime(Real a_time)
329 {
330  Vector<MGLevelOp<T>* > ops = m_solver->getAllOperators();
331  for (int iop = 0; iop < ops.size(); iop++)
332  {
333 
334  TGAHelmOp<T>* helmop = (TGAHelmOp<T>*) ops[iop];
335  helmop->setTime(a_time);
336  }
337 }
338 
339 /*****/
340 template <class T>
343  const AMRLevelOpFactory<T> & a_opFact,
344  const ProblemDomain& a_level0Domain,
345  const Vector<int>& a_refRat,
346  int a_numLevels,
347  int a_verbosity)
348 {
349  //cast to remove const because base class definition is weird
350  m_verbosity = a_verbosity;
351  m_level0Domain = a_level0Domain;
352  m_refRat = a_refRat;
353  m_solver = a_solver;
354  m_numLevels = a_numLevels;
355  if (m_numLevels < 0)
356  {
357  m_numLevels = a_refRat.size();
358  }
359 
360  m_ops.resize(m_numLevels);
361  Vector< AMRLevelOp<T> * >& amrops = m_solver->getAMROperators();
362  for (int ilev = 0; ilev < m_numLevels; ilev++)
363  {
364  m_ops[ilev] = dynamic_cast<TGAHelmOp<T>* >(amrops[ilev]);
365  if (m_ops[ilev]==NULL)
366  {
367  MayDay::Error("dynamic cast failed---is that operator really a TGAHelmOp?");
368  }
369  }
370 
371  setMu();
372  m_dataCreated = false;
373 }
374 template <class T>
375 void AMRTGA<T>::
377 {
378  Real tgaEpsilon = 1.e-12;
379 #ifdef CH_USE_FLOAT
380  tgaEpsilon = sqrt(tgaEpsilon);
381 #endif
382  Real a = 2.0 - sqrt(2.0) - tgaEpsilon;
383  m_mu1 = (a - sqrt( a*a - 4.0*a + 2.0))/2.0 ;
384  m_mu2 = (a + sqrt( a*a - 4.0*a + 2.0))/2.0 ;
385  m_mu3 = (1.0 - a);
386  m_mu4 = 0.5 - a;
387  if (m_verbosity > 4)
388  {
389  pout() << " AMRTGA:: epsilon = " << tgaEpsilon << std::endl;
390  }
391 }
392 
393 template <class T>
394 void AMRTGA<T>::
396  int a_lbase,
397  int a_lmax)
398 {
399  m_dataCreated = true;
400  m_rhst.resize(a_source.size(), NULL);
401  m_srct.resize(a_source.size(), NULL);
402  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
403  {
404  m_rhst[ilev] = new T();
405  m_srct[ilev] = new T();
406  m_ops[ilev]->create(*m_rhst[ilev], *a_source[ilev]);
407  m_ops[ilev]->create(*m_srct[ilev], *a_source[ilev]);
408  }
409 }
410 /*****/
411 template <class T>
413 {
414  for (int ilev = 0; ilev < m_rhst.size(); ilev++)
415  {
416  if (m_rhst[ilev] != NULL)
417  {
418  delete m_rhst[ilev];
419  delete m_srct[ilev];
420  m_rhst[ilev] = NULL;
421  m_srct[ilev] = NULL;
422  }
423  }
424 }
425 template <class T>
426 void AMRTGA<T>::
428  Vector<T* >& a_phiOld,
429  Vector<T* >& a_source,
430  const Real& a_dt,
431  int a_lbase,
432  int a_lmax,
433  Real a_told)
434 {
435  if (!m_dataCreated)
436  {
437  createData(a_source, a_lbase, a_lmax);
438  }
439 
440  if (m_verbosity > 3)
441  {
442  pout() << " AMRTGA:: starting mu4 operation" << std::endl;
443  }
444 
445  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
446  {
447  m_ops[ilev]->setToZero(*m_srct[ilev]);
448  m_ops[ilev]->incr(*m_srct[ilev], *a_source[ilev], 1.0);
449 
450  //this sets srct to S/a
451  m_ops[ilev]->divideByIdentityCoef(*m_srct[ilev]);
452  }
453 
454  //this operation is homogeneous and therefore does not need time set
455  //this makes rhs hold (k*a I + mu4 L) S/a
456  applyHelm(m_rhst, m_srct, a_lbase, a_lmax, m_mu4, a_dt, true);
457  //from here on k is kappa and L is kappa L
458  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
459  {
460  //this makes rhs hold dt(k*a I + mu4 L) S/a
461  m_ops[ilev]->scale( *m_rhst[ilev], a_dt);
462  }
463 
464  if (m_verbosity > 3)
465  {
466  pout() << " AMRTGA:: starting mu3 operation" << std::endl;
467  }
468  //set time to t_old
469  setTime(a_told);
470  //this makes a_phiNew hold (k*a I + mu3 L) phi^n
471  applyHelm(a_phiNew, a_phiOld, a_lbase, a_lmax, m_mu3, a_dt, false);
472 
473 
474 
475  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
476  {
477  //this makes rhs hold [(k*a I + mu3 L) phi^n + dt(k*a I + mu4 L) S/a]
478  m_ops[ilev]->incr(*m_rhst[ilev], *a_phiNew[ilev], 1.0);
479  }
480 
481  if (m_verbosity > 2)
482  {
483  pout() << " AMRTGA:: starting mu2 operation" << std::endl;
484  }
485  //set time to t_intermediate
486  setTime(a_told + a_dt*(m_mu2 + m_mu3));
487  //this makes phinew = phiold if using phiold as init guess
488  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
489  {
490  m_ops[ilev]->assign(*a_phiNew[ilev], *a_phiOld[ilev]);
491  }
492  solveHelm(a_phiNew, m_rhst, a_lbase, a_lmax, m_mu2, a_dt);
493  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
494  {
495  //this puts the answer into rhst so we can do the final solve
496  m_ops[ilev]->assign(*m_rhst[ilev], *a_phiNew[ilev]);
497  }
498 
499  //this makes rhst hold k*a[(k*a I - mu2 L)^-1 (rhs)]
500  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
501  {
502  m_ops[ilev]->diagonalScale(*m_rhst[ilev]);
503  }
504 
505  if (m_verbosity > 2)
506  {
507  pout() << " AMRTGA:: starting mu1 operation" << std::endl;
508  }
509  //this happens at tnew
510  setTime(a_told + a_dt);
511  //this makes phinew = phiold if using phiold as init guess
512  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
513  {
514  m_ops[ilev]->assign(*a_phiNew[ilev], *a_phiOld[ilev]);
515  }
516  solveHelm(a_phiNew, m_rhst, a_lbase, a_lmax, m_mu1, a_dt);
517 }
518 /*******/
519 template <class T>
520 void AMRTGA<T>::
522  Vector<T*>& a_phi,
523  int a_lbase,
524  int a_lmax,
525  Real a_mu,
526  Real a_dt,
527  bool a_homogeneousBC)
528 {
529  Real factor = a_mu*a_dt;
530 
531  resetAlphaAndBeta(1.0, factor);
532 
533  m_solver->computeAMROperator(a_ans,
534  a_phi,
535  a_lmax,
536  a_lbase,
537  a_homogeneousBC);
538 
539 // for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
540 // {
541 // m_ops[ilev]->scale(*a_ans[ilev], -1);
542 // }
543 }
544 /*******/
545 template <class T>
546 void AMRTGA<T>::
548  Vector<T*>& a_rhs,
549  int a_lbase,
550  int a_lmax,
551  Real a_mu,
552  Real a_dt)
553 {
554  Real factor = -a_mu*a_dt;
555 
556  resetAlphaAndBeta(1.0, factor);
557 
558  m_solver->solveNoInit(a_ans,
559  a_rhs,
560  a_lmax,
561  a_lbase,
562  false);
563  if (m_solver->m_exitStatus==1 && m_verbosity>3)
564  {
565  pout() << "AMRTGA:: WARNING: solver exitStatus == 1" << std::endl;
566  }
567 }
568 
569 #include "NamespaceFooter.H"
570 #endif
LevelTGAHelmOp(bool a_isTimeIndependent)
Definition: AMRTGA.H:174
std::ostream & pout()
Use this in place of std::cout for program output.
virtual AMRLevelOp< T > * AMRnewOp(const ProblemDomain &a_indexSpace)=0
A reference-counting handle class.
Definition: RefCountedPtr.H:66
TGAHelmOp & operator=(const TGAHelmOp &)
virtual void kappaScale(T &a_rhs)
for eb only. kappa weight the rhs but do not multiply by the identity coefficient ...
Definition: AMRTGA.H:77
A class to facilitate interaction with physical boundary conditions.
Definition: ProblemDomain.H:130
Vector< T * > m_srct
Definition: AMRTGA.H:288
virtual void diagonalScale(T &a_rhs, bool a_kappaWeighted)
Definition: AMRTGA.H:71
Real m_mu1
Definition: AMRTGA.H:284
void oneStep(Vector< T * > &a_phiNew, Vector< T * > &a_phiOld, Vector< T * > &a_source, const Real &a_dt, int a_lbase, int a_lmax, Real a_timeOld=0)
Definition: AMRTGA.H:427
int m_verbosity
Definition: AMRTGA.H:285
Real m_mu4
Definition: AMRTGA.H:284
Vector< TGAHelmOp< T > * > m_ops
Definition: AMRTGA.H:280
TGAHelmOp()
Definition: AMRTGA.H:37
virtual void getFlux(TFlux &a_flux, const T &a_data, const Box &a_grid, const DataIndex &a_dit, Real a_scale)=0
virtual void setTime(Real a_time)
Definition: AMRTGA.H:105
AMRTGA()
weak construction is bad. Ref Counted pointers are your friends.
Definition: AMRTGA.H:302
void resetAlphaAndBeta(const Real &a_alpha, const Real &a_beta)
Definition: AMRTGA.H:313
Definition: AMRMultiGrid.H:306
void solveHelm(Vector< T * > &a_ans, Vector< T * > &a_rhs, int a_lbase, int a_lmax, Real a_mu, Real a_dt)
Definition: AMRTGA.H:547
virtual void diagonalScale(T &a_rhs)
Definition: AMRTGA.H:86
virtual void divideByIdentityCoef(T &a_rhs)=0
Definition: AMRTGA.H:31
virtual void setAlphaAndBeta(const Real &a_alpha, const Real &a_beta)=0
Definition: AMRMultiGrid.H:35
RefCountedPtr< AMRMultiGrid< T > > m_solver
Definition: AMRTGA.H:283
virtual ~LevelTGAHelmOp()
Destructor.
Definition: AMRTGA.H:180
void createData(Vector< T * > &a_source, int a_lbase, int a_lmax)
Definition: AMRTGA.H:395
bool m_isTimeDependent
Definition: AMRTGA.H:146
void setTime(Real time)
Definition: AMRTGA.H:328
virtual void applyOpNoBoundary(T &a_ans, const T &a_phi)=0
Real m_mu3
Definition: AMRTGA.H:284
void applyHelm(Vector< T * > &a_ans, Vector< T * > &a_source, int a_lbase, int a_lmax, Real a_mu, Real a_dt, bool a_homogeneousBC)
Definition: AMRTGA.H:521
double Real
Definition: REAL.H:33
TGAHelmOp< T > * newOp(const ProblemDomain &a_indexSpace, const AMRLevelOpFactory< T > &a_opFact)
Definition: AMRTGA.H:211
Vector< T * > m_rhst
Definition: AMRTGA.H:287
LevelTGAHelmOp()
Definition: AMRTGA.H:165
AMRTGA(const AMRTGA< T > &a_opin)
Definition: AMRTGA.H:291
void operator=(const AMRTGA< T > &a_opin)
Definition: AMRTGA.H:296
TGAHelmOp(bool a_isTimeDependent)
Definition: AMRTGA.H:47
void createData(hid_t &a_dataset, hid_t &a_dataspace, HDF5Handle &handle, const std::string &name, hid_t type, hsize_t size)
void setMu()
Definition: AMRTGA.H:376
static void Error(const char *const a_msg=m_nullString, int m_exitCode=CH_DEFAULT_ERROR_CODE)
Print out message to cerr and exit with the specified exit code.
Vector< int > m_refRat
Definition: AMRTGA.H:282
bool isTimeDependent() const
Returns true if the operator is time-dependent, false otherwise.
Definition: AMRTGA.H:124
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:465
Definition: DataIndex.H:112
ProblemDomain m_level0Domain
Definition: AMRTGA.H:281
Real m_mu2
Definition: AMRTGA.H:284
bool m_dataCreated
Definition: AMRTGA.H:286
Definition: AMRTGA.H:159
virtual void fillGrad(const T &a_phi)=0
These functions are part of the LevelTGA interface......
size_t size() const
Definition: Vector.H:177
~AMRTGA()
Definition: AMRTGA.H:412
virtual T & identityCoef()
Allows access to the identity coefficient data for the operator.
Definition: AMRTGA.H:130
virtual void setTime(Real a_oldTime, Real a_mu, Real a_dt)
Definition: AMRTGA.H:118
Definition: AMRMultiGrid.H:231
int m_numLevels
Definition: AMRTGA.H:285
virtual ~TGAHelmOp()
Destructor.
Definition: AMRTGA.H:54
Definition: AMRTGA.H:206
T * m_phonyIdentityCoef
Definition: AMRTGA.H:150