Chombo + EB  3.2
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  Template implementation of the TGA algorithm
201  to solve the heat equation.
202  **/
203 template <class T>
204 class AMRTGA
205 {
206 public:
207 
208  ///
209  TGAHelmOp<T>* newOp(const ProblemDomain& a_indexSpace,
210  const AMRLevelOpFactory<T>& a_opFact)
211  {
212  AMRLevelOpFactory<T>& opFact = (AMRLevelOpFactory<T>&) a_opFact;
213  TGAHelmOp<T>* retval = (TGAHelmOp<T>*) opFact.AMRnewOp(a_indexSpace);
214  return retval;
215  }
216 
217  ///
218  ~AMRTGA();
219 
220  ///
221  /**
222  **/
223  AMRTGA(const RefCountedPtr<AMRMultiGrid<T> > & a_solver,
224  const AMRLevelOpFactory<T>& a_factory,
225  const ProblemDomain& a_level0Domain,
226  const Vector<int>& a_refRat,
227  int a_numLevels = -1,
228  int a_verbosity = 0);
229 
230  ///
231  /**
232  This advances a parabolic pde from a_phiOld to a_phiNew using TGA
233  on a non-moving domain with source term a_source
234  If your operator is time-dependent, be sure to send the time.
235  **/
236  void oneStep(Vector<T*>& a_phiNew,
237  Vector<T*>& a_phiOld,
238  Vector<T*>& a_source,
239  const Real& a_dt,
240  int a_lbase,
241  int a_lmax,
242  Real a_timeOld = 0);
243 
244  ///
245  void resetAlphaAndBeta(const Real& a_alpha,
246  const Real& a_beta);
247 
248  void setTime(Real time);
249 
250 protected:
251  void solveHelm(Vector<T*>& a_ans,
252  Vector<T*>& a_rhs,
253  int a_lbase,
254  int a_lmax,
255  Real a_mu,
256  Real a_dt);
257 
258  void applyHelm(Vector<T*>& a_ans,
259  Vector<T*>& a_source,
260  int a_lbase,
261  int a_lmax,
262  Real a_mu,
263  Real a_dt,
264  bool a_homogeneousBC);
265 
266  void setMu();
267 
268  void createData(Vector<T* >& a_source,
269  int a_lbase,
270  int a_lmax);
271 
272 private:
273  //You do not own these operators!! don't delete it. the memory is
274  //owned by the solver
284 
285  //copy constructor and operator= disallowed for all the usual reasons
286  AMRTGA(const AMRTGA<T>& a_opin)
287  {
288  MayDay::Error("invalid operator");
289  }
290 
291  void operator=(const AMRTGA<T>& a_opin)
292  {
293  MayDay::Error("invalid operator");
294  }
295 
296  /// weak construction is bad. Ref Counted pointers are your friends.
298  {
299  MayDay::Error("invalid operator");
300  }
301 };
302 
303 /*****/
304 template <class T>
305 void AMRTGA<T>::
306 resetAlphaAndBeta(const Real& a_alpha,
307  const Real& a_beta)
308 {
309  Vector<MGLevelOp<T>* > ops = m_solver->getAllOperators();
310  for (int iop = 0; iop < ops.size(); iop++)
311  {
312 
313  TGAHelmOp<T>* helmop = (TGAHelmOp<T>*) ops[iop];
314  helmop->setAlphaAndBeta(a_alpha, a_beta);
315  }
316 }
317 
318 /*****/
319 template <class T>
320 void AMRTGA<T>::
321 setTime(Real a_time)
322 {
323  Vector<MGLevelOp<T>* > ops = m_solver->getAllOperators();
324  for (int iop = 0; iop < ops.size(); iop++)
325  {
326 
327  TGAHelmOp<T>* helmop = (TGAHelmOp<T>*) ops[iop];
328  helmop->setTime(a_time);
329  }
330 }
331 
332 /*****/
333 template <class T>
336  const AMRLevelOpFactory<T> & a_opFact,
337  const ProblemDomain& a_level0Domain,
338  const Vector<int>& a_refRat,
339  int a_numLevels,
340  int a_verbosity)
341 {
342  //cast to remove const because base class definition is weird
343  m_verbosity = a_verbosity;
344  m_level0Domain = a_level0Domain;
345  m_refRat = a_refRat;
346  m_solver = a_solver;
347  m_numLevels = a_numLevels;
348  if (m_numLevels < 0)
349  {
350  m_numLevels = a_refRat.size();
351  }
352 
353  m_ops.resize(m_numLevels);
354  Vector< AMRLevelOp<T> * >& amrops = m_solver->getAMROperators();
355  for (int ilev = 0; ilev < m_numLevels; ilev++)
356  {
357  m_ops[ilev] = dynamic_cast<TGAHelmOp<T>* >(amrops[ilev]);
358  if (m_ops[ilev]==NULL)
359  {
360  MayDay::Error("dynamic cast failed---is that operator really a TGAHelmOp?");
361  }
362  }
363 
364  setMu();
365  m_dataCreated = false;
366 }
367 template <class T>
368 void AMRTGA<T>::
370 {
371  Real tgaEpsilon = 1.e-12;
372 #ifdef CH_USE_FLOAT
373  tgaEpsilon = sqrt(tgaEpsilon);
374 #endif
375  Real a = 2.0 - sqrt(2.0) - tgaEpsilon;
376  m_mu1 = (a - sqrt( a*a - 4.0*a + 2.0))/2.0 ;
377  m_mu2 = (a + sqrt( a*a - 4.0*a + 2.0))/2.0 ;
378  m_mu3 = (1.0 - a);
379  m_mu4 = 0.5 - a;
380  if (m_verbosity > 4)
381  {
382  pout() << " AMRTGA:: epsilon = " << tgaEpsilon << std::endl;
383  }
384 }
385 
386 template <class T>
387 void AMRTGA<T>::
389  int a_lbase,
390  int a_lmax)
391 {
392  m_dataCreated = true;
393  m_rhst.resize(a_source.size(), NULL);
394  m_srct.resize(a_source.size(), NULL);
395  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
396  {
397  m_rhst[ilev] = new T();
398  m_srct[ilev] = new T();
399  m_ops[ilev]->create(*m_rhst[ilev], *a_source[ilev]);
400  m_ops[ilev]->create(*m_srct[ilev], *a_source[ilev]);
401  }
402 }
403 /*****/
404 template <class T>
406 {
407  for (int ilev = 0; ilev < m_rhst.size(); ilev++)
408  {
409  if (m_rhst[ilev] != NULL)
410  {
411  delete m_rhst[ilev];
412  delete m_srct[ilev];
413  m_rhst[ilev] = NULL;
414  m_srct[ilev] = NULL;
415  }
416  }
417 }
418 template <class T>
419 void AMRTGA<T>::
421  Vector<T* >& a_phiOld,
422  Vector<T* >& a_source,
423  const Real& a_dt,
424  int a_lbase,
425  int a_lmax,
426  Real a_told)
427 {
428  if (!m_dataCreated)
429  {
430  createData(a_source, a_lbase, a_lmax);
431  }
432 
433  if (m_verbosity > 3)
434  {
435  pout() << " AMRTGA:: starting mu4 operation" << std::endl;
436  }
437 
438  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
439  {
440  m_ops[ilev]->setToZero(*m_srct[ilev]);
441  m_ops[ilev]->incr(*m_srct[ilev], *a_source[ilev], 1.0);
442 
443  //this sets srct to S/a
444  m_ops[ilev]->divideByIdentityCoef(*m_srct[ilev]);
445  }
446 
447  //this operation is homogeneous and therefore does not need time set
448  //this makes rhs hold (k*a I + mu4 L) S/a
449  applyHelm(m_rhst, m_srct, a_lbase, a_lmax, m_mu4, a_dt, true);
450  //from here on k is kappa and L is kappa L
451  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
452  {
453  //this makes rhs hold dt(k*a I + mu4 L) S/a
454  m_ops[ilev]->scale( *m_rhst[ilev], a_dt);
455  }
456 
457  if (m_verbosity > 3)
458  {
459  pout() << " AMRTGA:: starting mu3 operation" << std::endl;
460  }
461 
462  //set time to t_old
463  setTime(a_told);
464  //this makes a_phiNew hold (k*a I + mu3 L) phi^n
465  // false means that we're using the inhomogenous form of the BCs
466  applyHelm(a_phiNew, a_phiOld, a_lbase, a_lmax, m_mu3, a_dt, false);
467 
468  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
469  {
470  //this makes rhs hold [(k*a I + mu3 L) phi^n + dt(k*a I + mu4 L) S/a]
471  m_ops[ilev]->incr(*m_rhst[ilev], *a_phiNew[ilev], 1.0);
472  }
473 
474  if (m_verbosity > 2)
475  {
476  pout() << " AMRTGA:: starting mu2 operation" << std::endl;
477  }
478  //set time to t_intermediate
479  setTime(a_told + a_dt*(m_mu2 + m_mu3));
480  //this makes phinew = phiold if using phiold as init guess
481  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
482  {
483  m_ops[ilev]->assignLocal(*a_phiNew[ilev], *a_phiOld[ilev]);
484  }
485  solveHelm(a_phiNew, m_rhst, a_lbase, a_lmax, m_mu2, a_dt);
486  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
487  {
488  //this puts the answer into rhst so we can do the final solve
489  m_ops[ilev]->assignLocal(*m_rhst[ilev], *a_phiNew[ilev]);
490  }
491 
492  //this makes rhst hold k*a[(k*a I - mu2 L)^-1 (rhs)]
493  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
494  {
495  m_ops[ilev]->diagonalScale(*m_rhst[ilev]);
496  }
497 
498  if (m_verbosity > 2)
499  {
500  pout() << " AMRTGA:: starting mu1 operation" << std::endl;
501  }
502  //this happens at tnew
503  setTime(a_told + a_dt);
504  //this makes phinew = phiold if using phiold as init guess
505  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
506  {
507  m_ops[ilev]->assignLocal(*a_phiNew[ilev], *a_phiOld[ilev]);
508  }
509  solveHelm(a_phiNew, m_rhst, a_lbase, a_lmax, m_mu1, a_dt);
510 }
511 /*******/
512 template <class T>
513 void AMRTGA<T>::
515  Vector<T*>& a_phi,
516  int a_lbase,
517  int a_lmax,
518  Real a_mu,
519  Real a_dt,
520  bool a_homogeneousBC)
521 {
522  Real factor = a_mu*a_dt;
523 
524  resetAlphaAndBeta(1.0, factor);
525 
526  m_solver->computeAMROperator(a_ans,
527  a_phi,
528  a_lmax,
529  a_lbase,
530  a_homogeneousBC);
531 
532 // for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
533 // {
534 // m_ops[ilev]->scale(*a_ans[ilev], -1);
535 // }
536 }
537 /*******/
538 template <class T>
539 void AMRTGA<T>::
541  Vector<T*>& a_rhs,
542  int a_lbase,
543  int a_lmax,
544  Real a_mu,
545  Real a_dt)
546 {
547  Real factor = -a_mu*a_dt;
548 
549  resetAlphaAndBeta(1.0, factor);
550 
551  m_solver->solveNoInit(a_ans,
552  a_rhs,
553  a_lmax,
554  a_lbase,
555  false);//zeroPhi
556  if (m_solver->m_exitStatus==1 && m_verbosity>3)
557  {
558  pout() << "AMRTGA:: WARNING: solver exitStatus == 1" << std::endl;
559  }
560 }
561 
562 #include "NamespaceFooter.H"
563 #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:173
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:141
Vector< T * > m_srct
Definition: AMRTGA.H:283
virtual void diagonalScale(T &a_rhs, bool a_kappaWeighted)
Definition: AMRTGA.H:71
Real m_mu1
Definition: AMRTGA.H:279
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:420
int m_verbosity
Definition: AMRTGA.H:280
Real m_mu4
Definition: AMRTGA.H:279
Vector< TGAHelmOp< T > * > m_ops
Definition: AMRTGA.H:275
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:297
void resetAlphaAndBeta(const Real &a_alpha, const Real &a_beta)
Definition: AMRTGA.H:306
Definition: AMRMultiGrid.H:308
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:540
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:39
RefCountedPtr< AMRMultiGrid< T > > m_solver
Definition: AMRTGA.H:278
virtual ~LevelTGAHelmOp()
Destructor.
Definition: AMRTGA.H:180
void createData(Vector< T * > &a_source, int a_lbase, int a_lmax)
Definition: AMRTGA.H:388
bool m_isTimeDependent
Definition: AMRTGA.H:146
void setTime(Real time)
Definition: AMRTGA.H:321
virtual void applyOpNoBoundary(T &a_ans, const T &a_phi)=0
Real m_mu3
Definition: AMRTGA.H:279
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:514
double Real
Definition: REAL.H:33
TGAHelmOp< T > * newOp(const ProblemDomain &a_indexSpace, const AMRLevelOpFactory< T > &a_opFact)
Definition: AMRTGA.H:209
Vector< T * > m_rhst
Definition: AMRTGA.H:282
LevelTGAHelmOp()
Definition: AMRTGA.H:165
AMRTGA(const AMRTGA< T > &a_opin)
Definition: AMRTGA.H:286
void operator=(const AMRTGA< T > &a_opin)
Definition: AMRTGA.H:291
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:369
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:277
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:276
Real m_mu2
Definition: AMRTGA.H:279
bool m_dataCreated
Definition: AMRTGA.H:281
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:192
~AMRTGA()
Definition: AMRTGA.H:405
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:233
int m_numLevels
Definition: AMRTGA.H:280
virtual ~TGAHelmOp()
Destructor.
Definition: AMRTGA.H:54
Definition: AMRTGA.H:204
T * m_phonyIdentityCoef
Definition: AMRTGA.H:150