Chombo + EB  3.0
AMRMultiGrid.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 // BVS, June 18, 2003
12 
13 // We can assume that template class T has null construction.
14 
15 #ifndef _AMRMULTIGRID_H_
16 #define _AMRMULTIGRID_H_
17 
18 #include "MultiGrid.H"
19 #include "REAL.H"
20 #include "Box.H"
21 #include "NoOpSolver.H"
22 #include "parstream.H"
23 #include "CH_Timer.H"
24 #include "Copier.H"
25 #include "SPMD.H"
26 #include "Misc.H"
27 
28 #include "NamespaceHeader.H"
29 
30 ///
31 /**
32  Operator class for AMR Multigrid
33  */
34 template <typename T>
35 class AMRLevelOp : public MGLevelOp<T>
36 {
37 public:
38 
39  ///
40  virtual void dumpAMR(Vector<T*>& a_data, string name)
41  {
42  }
43 
44  virtual void dumpLevel(T& a_data, string name)
45  {
46  }
47 
48  //! Constructor.
50  :MGLevelOp<T>()
51  {
52  }
53 
54  virtual void dumpStuff(Vector<T*> data, string filename)
55  {
56  }
57 
58  //! Destructor.
59  virtual ~AMRLevelOp()
60  {
61  }
62 
63  virtual Real AMRNorm(const T& a_coarResid,
64  const T& a_fineResid,
65  const int& a_refRat,
66  const int& a_ord)
67  {
68  return this->norm(a_coarResid, 0);
69  }
70 
71  //debugging hook
72  virtual void outputLevel(T& a_rhs, string& a_name)
73  {
74  }
75 
76  //debugging hook
77  virtual void outputAMR(Vector<T*>& a_rhs, string& a_name)
78  {
79  }
80  ///
81  /**
82  return the refinement ratio to next coarser level.
83  return 1 when there are no coarser AMRLevelOp objects.
84  */
85  virtual int refToCoarser() = 0;
86 
87  ///
88  /**
89  a_residual = a_rhs - L(a_phiFine, a_phi, a_phiCoarse)
90  */
91  virtual void AMRResidual(T& a_residual, const T& a_phiFine, const T& a_phi,
92  const T& a_phiCoarse, const T& a_rhs,
93  bool a_homogeneousDomBC,
94  AMRLevelOp<T>* a_finerOp) = 0;
95 
96  ///
97  /**
98  a_residual = a_rhs - L^nf(a_phi, a_phiCoarse)
99  assume no finer AMR level
100  */
101  virtual void AMRResidualNF(T& a_residual, const T& a_phi, const T& a_phiCoarse,
102  const T& a_rhs, bool a_homogeneousBC) = 0;
103 
104  ///
105  /**
106  a_residual = a_rhs - L(a_phiFine, a_phi)
107  assume no coarser AMR level
108  */
109  virtual void AMRResidualNC(T& a_residual, const T& a_phiFine, const T& a_phi,
110  const T& a_rhs, bool a_homogeneousBC,
111  AMRLevelOp<T>* a_finerOp) = 0;
112 
113 
114 
115  ///
116  /**
117  Apply the AMR operator, including coarse-fine matching
118  */
119  virtual void AMROperator(T& a_LofPhi,
120  const T& a_phiFine, const T& a_phi,
121  const T& a_phiCoarse,
122  bool a_homogeneousDomBC,
123  AMRLevelOp<T>* a_finerOp) = 0;
124 
125  ///
126  /**
127  Apply the AMR operator, including coarse-fine matching.
128  assume no finer AMR level
129  */
130  virtual void AMROperatorNF(T& a_LofPhi,
131  const T& a_phi,
132  const T& a_phiCoarse,
133  bool a_homogeneousBC) = 0;
134 
135  ///
136  /**
137  Apply the AMR operator, including coarse-fine matching
138  assume no coarser AMR level
139  */
140  virtual void AMROperatorNC(T& a_LofPhi,
141  const T& a_phiFine,
142  const T& a_phi,
143  bool a_homogeneousBC,
144  AMRLevelOp<T>* a_finerOp) = 0;
145 
146 
147  ///
148  /** a_resCoarse = I[h-2h]( a_residual - L(a_correction, a_coarseCorrection)) */
149  virtual void AMRRestrict(T& a_resCoarse, const T& a_residual, const T& a_correction,
150  const T& a_coarseCorrection) = 0;
151 
152 
153 
154  ///
155  /** a_correction += I[2h->h](a_coarseCorrection) */
156  virtual void AMRProlong(T& a_correction, const T& a_coarseCorrection) = 0;
157 
158 
159 
160  ///
161  /** a_residual = a_residual - L(a_correction, a_coarseCorrection) */
162  virtual void AMRUpdateResidual(T& a_residual, const T& a_correction,
163  const T& a_coarseCorrection) = 0;
164 
165 
166  ///
167  /**
168  */
169  virtual void createCoarsened(T& a_lhs,
170  const T& a_rhs,
171  const int& a_refRat) = 0;
172 
173 
174 
175  //===================================================================
176  // optional optimizations for an AMRLevelOp. These are not pure virtual
177  // functions, since we can build the equivalent algorithmic components from
178  // pure virtual functions. The AMRMultiGrid algorithm actually calls *these*
179  // functions, which a smart operator can perform in faster ways.
180  //===================================================================
181 
182 
183 
184  virtual void buildCopier(Copier& a_copier, const T& a_lhs, const T& a_rhs)
185  {
186  }
187 
188  virtual void assignCopier(T& a_lhs, const T& a_rhs, const Copier& a_copier)
189  {
190  this->assign(a_lhs, a_rhs);
191  }
192 
193  virtual void zeroCovered(T& a_lhs, T& a_rhs, const Copier& a_copier)
194  {
195  this->setToZero(a_rhs);
196  this->assignCopier(a_lhs, a_rhs, a_copier);
197  }
198  virtual Real localMaxNorm(const T& a_phi)
199  {
200  return this->norm(a_phi, 0);
201  }
202 
203  /** optimization of AMRProlong that sends in the existing temporary and copier */
204  virtual void AMRProlongS(T& a_correction, const T& a_coarseCorrection,
205  T& a_temp, const Copier& a_copier)
206  {
207  AMRProlong(a_correction, a_coarseCorrection);
208  }
209  virtual void AMRRestrictS(T& a_resCoarse, const T& a_residual, const T& a_correction,
210  const T& a_coarseCorrection, T& scratch)
211  {
212  AMRRestrict(a_resCoarse, a_residual, a_correction, a_coarseCorrection);
213  }
214 
215  virtual unsigned int orderOfAccuracy(void) const
216  {
217  return 2;
218  }
219 
220  /// This routine is for operators with orderOfAccuracy()>2.
221  virtual void enforceCFConsistency(T& a_coarseCorrection, const T& a_correction)
222  {
223  }
224 };
225 
226 ///
227 /**
228  Factory to create AMRLevelOps
229  */
230 template <class T>
232 {
233 public:
235  {
236  }
237 
238  ///
239  /**
240  return a new operator. this is done with a new call.
241  caller is responsible for deletion
242  */
243  virtual AMRLevelOp<T>* AMRnewOp(const ProblemDomain& a_indexSpace)=0;
244 
245  ///
246  /**
247  return refinement ratio to next finer level.
248  */
249  virtual int refToFiner(const ProblemDomain& a_indexSpace) const =0;
250 
251 };
252 
253 //! \class AMRMultiGridInspector
254 //! This base class allows one to construct methods for inspecting the
255 //! multigrid algorithm at each of its steps.
256 template <class T>
258 {
259  public:
260 
261  //! Base class constructor. This must be called by all subclasses.
263  {
264  }
265 
266  //! Destructor.
268  {
269  }
270 
271  //! Override this method to keep track of a multigrid residual computed
272  //! during a multigrid iteration at the given levels.
273  //! \param a_residuals An array containing the residuals computed by the multigrid
274  //! algorithm at each level in the range [a_minLevel, a_maxLevel].
275  //! \param a_minLevel The lowest AMR level at which residuals were computed.
276  //! \param a_maxLevel The highest AMR level at which residuals were computed.
277  //! \param a_iter The multigrid iteration number.
278  virtual void recordResiduals(const Vector<T*>& a_residuals,
279  int a_minLevel,
280  int a_maxLevel,
281  int a_iter) = 0;
282 
283  //! Override this method to keep track of a multigrid correction computed
284  //! during a V cycle at the given level.
285  //! \param a_corrections An array containing the corrections computed during a V cycle
286  //! at each level in the range [a_minLevel, a_maxLevel].
287  //! \param a_minLevel The lowest AMR level at which corrections were computed.
288  //! \param a_maxLevel The highest AMR level at which corrections were computed.
289  //! \param a_iter The multigrid iteration number.
290  virtual void recordCorrections(const Vector<T*>& a_corrections,
291  int a_minLevel,
292  int a_maxLevel,
293  int a_iter) = 0;
294 
295  private:
296 
299 };
300 
301 ///
302 /**
303  Class to solve elliptic equations using the Martin and Cartwright algorithm.
304  */
305 template <class T>
307 {
308 public:
309 
310  AMRMultiGrid();
311 
312  virtual ~AMRMultiGrid();
313 
314  void outputAMR(Vector<T*>& a_data, string& a_name, int a_lmax, int a_lbase);
315 
316  ///
317  /**
318  Define the solver.
319  a_coarseDomain is the index space on the coarsest AMR level.
320  a_factory is the operator factory through which all special information is conveyed.
321  a_bottomSolver is the solver to be used at the termination of multigrid coarsening.
322  It is the client's responsibility to free up the dynamically-allocated memory.
323  a_numLevels is the number of AMR levels.
324  */
325  virtual void define(const ProblemDomain& a_coarseDomain,
326  AMRLevelOpFactory<T>& a_factory,
327  LinearSolver<T>* a_bottomSolver,
328  int a_numLevels);
329 
330  //! Add an inspector to the list of inspectors maintained by this AMRMultiGrid
331  //! instance. It will be given the opportunity to record intermediate data.
333  {
334  CH_assert(!a_inspector.isNull());
335  m_inspectors.push_back(a_inspector);
336  }
337 
338  ///
339  /**
340  Solve L(phi) = rho from l_base to l_max. To solve over all levels,
341  l_base = 0 and l_max = max_level = numLevels-1.
342  */
343  virtual void solve(Vector<T*>& a_phi, const Vector<T*>& a_rhs,
344  int l_max, int l_base, bool a_zeroPhi=true,
345  bool forceHomogeneous = false);
346 
347  ///
348  /** same as "solve" except user has taken the reponsibility of having previously
349  called "init" so solver can allocate temporary holders.
350  */
351  virtual void solveNoInit(Vector<T*>& a_phi, const Vector<T*>& a_rhs,
352  int l_max, int l_base, bool a_zeroPhi=true,
353  bool forceHomogeneous = false);
354 
355  ///use if you want final residual
356  virtual void solveNoInitResid(Vector<T*>& a_phi, Vector<T*>& a_finalResid,
357  const Vector<T*>& a_rhs,
358  int l_max, int l_base, bool a_zeroPhi=true,
359  bool forceHomogeneous = false);
360 
361  void relaxOnlyHomogeneous(Vector<T*>& a_phi, const Vector<T*>& a_rhs,
362  int l_max, int l_base);
363 
364  virtual void AMRVCycle(Vector<T*>& a_correction,
365  Vector<T*>& a_residual,
366  int l, int l_max, int l_base);
367 
368  void setMGCycle(int a_numMG);
369 
370  void init(const Vector<T*>& a_phi, const Vector<T*>& a_rhs,
371  int l_max, int l_base);
372  //init messes with multigrid depth. this puts it back
373  void revert(const Vector<T*>& a_phi, const Vector<T*>& a_rhs,
374  int l_max, int l_base);
375 
376 
381  /// max no. of coarsenings -- -1 (default) means coarsen as far as possible
382  /** If using a value besides the default, need to set it _before_
383  define function is called */
385  AMRLevelOp<T>& levelOp(int level);
386  // default m_convergenceMetric = 0.: initial residual will be set to
387  // result of computeAMRResidual.
388  // if m_convergenceMetric > 0., then initial residual will be set to
389  // m_convergenceMetric.
391 
392  // used to give an additional cushion in the EPS used for bottom solves
394 
395  ///
396  /**
397  resid = L(phi) - rhs
398  */
400  Vector<T*>& a_phi,
401  const Vector<T*>& a_rhs,
402  int l_max,
403  int l_base,
404  bool a_homogeneousBC=false,
405  bool a_computeNorm=true);
406 
407 
408  /** just return the normed value of computeAMRResidual. used for benchmarking */
410  const Vector<T*>& a_rhs,
411  int l_max,
412  int l_min);
413 
414  ///
415  /**
416  lph = L(phi)
417  */
418  void computeAMROperator(Vector<T*>& a_lph,
419  Vector<T*>& a_phi,
420  int l_max,
421  int l_base,
422  bool a_homogeneousBC=false);
423 
424  ///
425  /**
426  For changing coefficients. Use at thy own peril.
427  */
432  {
433  return m_op;
434  }
435 
436  ///
437  /**
438  Set parameters of the solve.
439  a_pre is the number of smoothings before averaging.
440  a_post is the number of smoothings after averaging.
441  a_bottom is the number of smoothings at the bottom level.
442  a_numMG = 1 for vcycle, =2 for wcycle (use 1).
443  a_itermax is the max number of v cycles.
444  a_hang is the minimum amount of change per vcycle.
445  a_eps is the solution tolerance.
446  a_normThresh is how close to zero eps*resid is allowed to get.
447  */
448  void setSolverParameters(const int& a_pre,
449  const int& a_post,
450  const int& a_bottom,
451  const int& a_numMG,
452  const int& a_iterMax,
453  const Real& a_eps,
454  const Real& a_hang,
455  const Real& a_normThresh);
456 
457 
458  /// set up bottom solver for internal MG solver
459  /**
460  This function is normally called by the solve(...) function.
461  However, it must be called if solve will not be called (in particular,
462  if only the V-cycle is being used)
463  */
464  void setBottomSolver(int l_max, int l_base);
465 
466  void setBottomSolverEpsCushion(Real a_bottomSolverEpsCushion);
467 
468 protected:
469 
470  void relax(T& phi, T& R, int depth, int nRelax = 2);
471 
472  void computeAMRResidualLevel(Vector<T*>& a_resid,
473  Vector<T*>& a_phi,
474  const Vector<T*>& a_rhs,
475  int l_max, int l_base, int ilev,
476  bool a_homogeneousBC);
477 
485 
487 
489 
491 
492  void clear();
493 
494 private:
495 
496  // A list of inspectors maintained by this instance.
498 
499  // Forbidden copiers.
502 };
503 
504 //*******************************************************
505 // AMRMultigrid Implementation
506 //*******************************************************
507 
508 //===================================================================
509 
510 template <class T>
511 void
512 AMRMultiGrid<T>::outputAMR(Vector<T*>& a_data, string& a_name, int a_lmax, int a_lbase)
513 {
514  Vector<T*> outputData;
515  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
516  {
517  outputData.push_back(a_data[ilev]);
518  }
519  m_op[a_lbase]->outputAMR(outputData, a_name);
520 }
521 
522 template <class T>
523 void
525  const int& a_post,
526  const int& a_bottom,
527  const int& a_numMG,
528  const int& a_iterMax,
529  const Real& a_eps,
530  const Real& a_hang,
531  const Real& a_normThresh)
532 {
533  m_solverParamsSet = true;
534  m_pre = a_pre;
535  m_post = a_post;
536  m_bottom = a_bottom;
537  m_eps = a_eps;
538  m_hang = a_hang;
539  m_normThresh = a_normThresh;
540  m_iterMax = a_iterMax;
541  for (int img = 0; img < m_mg.size(); img++)
542  {
543  m_mg[img]->m_pre = a_pre;
544  m_mg[img]->m_post = a_post;
545  m_mg[img]->m_bottom = a_bottom;
546  }
547  setMGCycle(a_numMG);
548  m_bottomSolverEpsCushion = 1.0;
549 }
550 template <class T>
553 {
554  Vector< MGLevelOp<T>* > retval;
555  for (int iop = 0; iop < m_op.size(); iop++)
556  {
557  MGLevelOp<T>* operPtr = (MGLevelOp<T>*) m_op[iop];
558  retval.push_back(operPtr);
559  }
560 
561  for (int img = 0; img < m_mg.size(); img++)
562  {
563  Vector< MGLevelOp<T>* > mgOps = m_mg[img]->getAllOperators();
564  retval.append(mgOps);
565  }
566  return retval;
567 }
568 
569 template <class T>
572 {
573  Vector< MGLevelOp<T>* > retval;
574  for (int iop = 0; iop < m_op.size(); iop++)
575  {
576  MGLevelOp<T>* operPtr = (MGLevelOp<T>*) m_op[iop];
577  retval.push_back(operPtr);
578  }
579  return retval;
580 }
581 
582 template <class T>
585 {
586  Vector< Vector< MGLevelOp<T>* > > retval(m_mg.size());
587 
588  for (int img = 0; img < m_mg.size(); img++)
589  {
590  retval[img] = m_mg[img]->getAllOperators();
591  }
592  return retval;
593 }
594 
595 template <class T>
597  :m_eps(1E-6),
598  m_hang(1E-15),
599  m_normThresh(1E-30),
600  m_imin(5),
601  m_iterMax(20),
602  m_verbosity(3),
603  m_pre(2),
604  m_post(2),
605  m_bottom(2),
606  m_numMG(1),
607  m_maxDepth(-1),
608  m_convergenceMetric(0.),
609  m_bottomSolverEpsCushion(1.0),
610  m_bottomSolver(NULL),
611  m_inspectors()
612 {
613  m_solverParamsSet = false;
614 
615  //m_maxDepth = 4;
616  //pout() << " AMRMultiGrid<T>::AMRMultiGrid() = " << m_maxDepth << endl;
617 }
618 template <class T>
620 {
621  for (int ilev = 0; ilev < m_op.size(); ilev++)
622  {
623  m_mg[ilev]->m_numMG = a_numMG;
624  m_mg[ilev]->m_cycle = a_numMG;
625  }
626  m_numMG = a_numMG;
627 }
628 template <class T>
629 void AMRMultiGrid<T>::relax(T& a_correction, T& a_residual, int depth, int a_numSmooth)
630 {
631  CH_TIME("AMRMultiGrid::relax");
632 
633 
634  if (m_op[depth]->refToCoarser() > 2)
635  {
636  // intermediate multigrid levels exist between this level and the
637  // next coarser AMR level, so do a mini v-cyle to smooth on those
638 
639  int intermediateDepth = 0;
640  int r = m_op[depth]->refToCoarser();
641  while (r>2)
642  {
643  r/=2;
644  intermediateDepth++;
645  }
646  int tmp = m_mg[depth]->m_depth;
647  m_mg[depth]->m_depth = intermediateDepth;
648  //ok, there is an intermediate multigrid level that is
649  // not an AMR level. We use regular MultiGrid smoothing algorithm for this.
650  // note that MultiGrid::cycle handles relaxation on this level
651  // as well, so there's no need to call the LinearOp's relax function.
652  m_mg[depth]->cycle(0, a_correction, a_residual);
653  m_mg[depth]->m_depth = tmp;
654  }
655  else
656  {
657  // no intermediate multigrid levels between AMR levels, so
658  // just call relax on this level and be done with it.
659  m_op[depth]->relax(a_correction, a_residual, a_numSmooth); //numSmoothDown
660  }
661 
662 }
663 /************/
664 template <class T>
666 {
667  CH_TIME("~AMRMultiGrid");
668  clear();
669 }
670 
671 template <class T>
673 {
674  return *(m_op[level]);
675 }
676 
677 /************/
678 template <class T>
680  Vector<T*>& a_phi,
681  const Vector<T*>& a_rhs,
682  int l_max,
683  int l_base,
684  bool a_homogeneousBC,
685  bool a_computeNorm)
686 {
687  CH_TIME("AMRMultiGrid::computeAMRResidual");
688 
689  Real rnorm = 0;
690  Real localNorm = 0;
691  for (int ilev = l_base; ilev <= l_max; ilev++)
692  {
693  //always used at top level where bcs are inhomogeneous
694  computeAMRResidualLevel(a_resid,
695  a_phi,
696  a_rhs,
697  l_max, l_base, ilev, a_homogeneousBC);
698  if (a_computeNorm)
699  {
700  if (ilev == l_max)
701  {
702  localNorm =m_op[ilev]->localMaxNorm(*a_resid[ilev]);
703  }
704  else
705  {
706  m_op[ilev]->zeroCovered(*a_resid[ilev], *m_resC[ilev+1], m_resCopier[ilev+1]);
707  localNorm = m_op[ilev]->localMaxNorm(*a_resid[ilev]);
708  }
709  const int normType = 0;
710  if (normType==2)
711  localNorm = m_op[ilev]->norm(*a_resid[ilev],normType);
712  rnorm = Max(localNorm, rnorm);
713  }
714  }
715 #ifdef CH_MPI
716  if (a_computeNorm)
717  {
718  CH_TIME("MPI_Allreduce");
719  Real recv;
720  int result = MPI_Allreduce(&rnorm, &recv, 1, MPI_CH_REAL,
721  MPI_MAX, Chombo_MPI::comm);
722  if (result != MPI_SUCCESS)
723  {
724  //bark!!!
725  MayDay::Error("sorry, but I had a communcation error on norm");
726  }
727  rnorm = recv;
728  }
729 #endif
730 
731  return rnorm; // if a_computeNorm is false, then this just returns zero.
732 }
733 
734 template <class T>
736  const Vector<T*>& a_rhs,
737  int l_max,
738  int l_min)
739 {
740  return computeAMRResidual(m_residual,
741  a_phi,
742  a_rhs,
743  l_max,
744  l_min,
745  false,
746  true);
747 
748 
749 }
750 /************/
751 template <class T>
753  Vector<T*>& a_phi,
754  int l_max,
755  int l_base,
756  bool a_homogeneousBC)
757 {
758  for (int ilev = l_base; ilev <= l_max; ilev++)
759  {
760  m_op[ilev]->create(*m_residual[ilev], *a_lph[ilev]);
761  m_op[ilev]->setToZero(*m_residual[ilev]);
762  }
763  computeAMRResidual(a_lph, a_phi, m_residual, l_max, l_base, a_homogeneousBC, false);
764  // fixing the bug of not negating the result --- Qinghai Zhang 12/10/2009
765  // In loving memory of the two weeks I lost over this bug!!!
766  for (int ilev = l_base; ilev <= l_max; ilev++)
767  {
768  m_op[ilev]->scale(*a_lph[ilev], -1.);
769  }
770 }
771 
772 
773 /************/
774 template <class T>
776  Vector<T*>& a_phi,
777  const Vector<T*>& a_rhs,
778  int l_max, int l_base, int ilev,
779  bool a_homogeneousBC)
780 {
781  CH_TIME("AMRMultiGrid<T>::computeAMRResidualLevel");
782  CH_assert(m_hasInitBeenCalled[ilev]=='t');
783 
784  //m_op[ilev]->setToZero(*(a_resid[l_max]));
785  if (l_max != l_base)
786  {
787  if (ilev == l_max)
788  {
789  m_op[l_max]->AMRResidualNF(*(a_resid[l_max]), *(a_phi[l_max]),
790  *(a_phi[l_max-1]), *(a_rhs[l_max]),
791  a_homogeneousBC);
792  }
793  else if (ilev == l_base)
794  {
795  if (l_base == 0)
796  {
797  m_op[l_base]->AMRResidualNC(*(a_resid[l_base]), *(a_phi[l_base+1]),
798  *(a_phi[l_base]), *(a_rhs[l_base]),
799  a_homogeneousBC, m_op[l_base+1]);
800  }
801  else
802  {
803  m_op[l_base]->AMRResidual(*a_resid[l_base], *a_phi[l_base+1], *a_phi[l_base],
804  *a_phi[l_base-1], *a_rhs[l_base],
805  a_homogeneousBC, m_op[l_base+1]);
806  }
807  }
808  else
809  {
810  m_op[ilev]->AMRResidual(*a_resid[ilev], *a_phi[ilev+1], *a_phi[ilev],
811  *a_phi[ilev-1], *a_rhs[ilev],
812  a_homogeneousBC, m_op[ilev+1]);
813  }
814  }
815  else
816  {
817  CH_assert(ilev == l_base);
818  if (l_base == 0)
819  {
820  m_op[l_max]->residual(*a_resid[l_max], *a_phi[l_max], *a_rhs[l_max],a_homogeneousBC);
821  }
822  else
823  {
824  m_op[l_max]->AMRResidualNF(*(a_resid[l_max]), *(a_phi[l_max]),
825  *(a_phi[l_max-1]), *(a_rhs[l_max]),
826  a_homogeneousBC);
827  }
828  }
829 
830 }
831 
832 template<class T>
833 void AMRMultiGrid<T>::solve(Vector<T*>& a_phi, const Vector<T*>& a_rhs,
834  int l_max, int l_base, bool a_zeroPhi,
835  bool a_forceHomogeneous)
836 {
837  CH_TIME("AMRMultiGrid::solve");
838  init(a_phi, a_rhs, l_max, l_base);
839  solveNoInit(a_phi, a_rhs, l_max, l_base, a_zeroPhi, a_forceHomogeneous);
840  //puts multigrid depth back where it started so the solver can be reused
841  revert(a_phi, a_rhs, l_max, l_base);
842 }
843 
844 template<class T>
846  const Vector<T*>& a_rhs,
847  int l_max, int l_base, bool a_zeroPhi,
848  bool a_forceHomogeneous)
849 {
850  Vector<T*> uberResidual(a_rhs.size());
851  int lowlim = l_base;
852  if (l_base > 0) // we need an ubercorrection one level lower than l_base
853  lowlim--;
854  // for (int ilev = l_base; ilev <= l_max; ilev++)
855  for (int ilev = lowlim; ilev <= l_max; ilev++)
856  {
857  uberResidual[ilev] = new T();
858  if (ilev >= l_base)
859  {
860  m_op[ilev]->create(*uberResidual[ilev], *a_rhs[ilev]);
861  }
862  }
863  solveNoInitResid(a_phi, uberResidual, a_rhs, l_max, l_base, a_zeroPhi, a_forceHomogeneous);
864  for (int i = lowlim; i <= l_max; i++)
865  {
866  m_op[i]->clear(*uberResidual[i]);
867  delete uberResidual[i];
868  }
869 }
870 
871 template<class T>
873  const Vector<T*>& a_rhs,
874  int l_max, int l_base, bool a_zeroPhi,
875  bool a_forceHomogeneous)
876 {
877  CH_TIMERS("AMRMultiGrid::solveNoInit");
878  CH_TIMER("AMRMultiGrid::AMRVcycle", vtimer);
879  setBottomSolver(l_max, l_base);
880 
881  CH_assert(l_base <= l_max);
882  CH_assert(a_rhs.size() == a_phi.size());
883 
884  //these correspond to the residual and correction
885  //that live in AMRSolver
886  Vector<T*> uberCorrection(a_rhs.size());
887 
888  int lowlim = l_base;
889  if (l_base > 0) // we need an ubercorrection one level lower than l_base
890  lowlim--;
891  // for (int ilev = l_base; ilev <= l_max; ilev++)
892  bool outputIntermediates = false;
893 
894  for (int ilev = lowlim; ilev <= l_max; ilev++)
895  {
896  uberCorrection[ilev] = new T();
897  m_op[ilev]->create(*uberCorrection[ilev], *a_phi[ilev]);
898  if (ilev >= l_base)
899  {
900  m_op[ilev]->create(*uberResidual[ilev], *a_rhs[ilev]);
901  m_op[ilev]->setToZero(*(uberResidual[ilev]));
902  }
903  m_op[ilev]->setToZero(*(uberCorrection[ilev]));
904  }
905  // m_op[0]->dumpStuff(uberResidual, string("initialRes.hdf5"));
906  if (a_zeroPhi)
907  for (int ilev = l_base; ilev <=l_max; ++ilev)
908  {
909  m_op[ilev]->setToZero(*(a_phi[ilev]));
910  }
911  //compute initial residual and initialize internal residual to it
912 
913  Real initial_rnorm = 0;
914  {
915  CH_TIME("Initial AMR Residual");
916  initial_rnorm=computeAMRResidual(uberResidual, a_phi, a_rhs, l_max, l_base, a_forceHomogeneous, true);
917  }
918 
919  if (m_convergenceMetric != 0.)
920  {
921  initial_rnorm = m_convergenceMetric;
922  }
923 
924  Real rnorm = initial_rnorm;
925  Real norm_last = 2*initial_rnorm;
926 
927  /// set bottom solver convergence norm and solver tolerance
928  m_bottomSolver->setConvergenceMetrics(initial_rnorm, m_bottomSolverEpsCushion*m_eps);
929 
930  int iter=0;
931  if (m_verbosity >= 2) ////
932  {
933  pout() << " AMRMultiGrid:: iteration = " << iter << ", residual norm = " << rnorm << std::endl;
934  }
935 
936  bool goNorm = rnorm > m_normThresh; //iterate if norm is not small enough
937  bool goRedu = rnorm > m_eps*initial_rnorm; //iterate if initial norm is not reduced enough
938  bool goIter = iter < m_iterMax; //iterate if iter < max iteration count
939  bool goHang = iter < m_imin || rnorm <(1-m_hang)*norm_last;//iterate if we didn't hang
940  while (goIter && goRedu && goHang && goNorm)
941  {
942 
943  // Report the residuals to the inspectors.
944  for (int i = 0; i < m_inspectors.size(); ++i)
945  m_inspectors[i]->recordResiduals(uberResidual, l_base, l_max, iter);
946 
947  if (outputIntermediates)
948  {
949  char strresname[100];
950  sprintf(strresname, "amrmg.res.iter.%03d", iter);
951  string nameres(strresname);
952  outputAMR(uberResidual, nameres, l_max, l_base);
953  }
954 
955  norm_last = rnorm;
956 
957  //this generates a correction from the current residual
958  CH_START(vtimer);
959  AMRVCycle(uberCorrection, uberResidual, l_max, l_max, l_base);
960  CH_STOP(vtimer);
961 
962 
963  char charname[100];
964  sprintf(charname, "resid_iter%03d.%dd.hdf5", iter, SpaceDim);
965  string sname(charname);
966  // m_op[0]->dumpAMR(uberResidual, sname);
967  // Report the corrections to the inspectors.
968  for (int i = 0; i < m_inspectors.size(); ++i)
969  m_inspectors[i]->recordCorrections(uberCorrection, l_base, l_max, iter);
970 
971  //increment phi by correction and reset correction to zero
972  for (int ilev = l_base; ilev <= l_max; ilev++)
973  {
974  if (outputIntermediates)
975  {
976  char strcorname[100];
977  sprintf(strcorname, "amrmg.phi.iter.%03d", iter);
978  string namecor(strcorname);
979  outputAMR(a_phi, namecor, l_max, l_base);
980  }
981 
982  m_op[ilev]->incr(*(a_phi[ilev]), *(uberCorrection[ilev]), 1.0);
983 
984  if (outputIntermediates)
985  {
986  char strcorname[100];
987  sprintf(strcorname, "amrmg.cor.iter.%03d", iter);
988  string namecor(strcorname);
989  outputAMR(uberCorrection, namecor, l_max, l_base);
990  }
991 
992  m_op[ilev]->setToZero(*(uberCorrection[ilev]));
993  // clean const out
994  }
995 
996 
997  // For solvers with accuracy higher than 2nd order
998  // consistency between levels has to be explicitly enforced.
999  // Qinghai Zhang
1000  if (m_op[0]->orderOfAccuracy()>2)
1001  {
1002  for (int ilev=l_max; ilev>l_base; ilev--)
1003  {
1004  m_op[ilev]->enforceCFConsistency(*a_phi[ilev-1], *a_phi[ilev]);
1005  }
1006  }
1007  //------end enforcing consistency.
1008 
1009  //recompute residual
1010  rnorm = computeAMRResidual(uberResidual, a_phi, a_rhs, l_max, l_base, a_forceHomogeneous, true);
1011  iter++;
1012  if (m_verbosity >= 3) ////
1013  {
1014  pout() << " AMRMultiGrid:: iteration = " << iter << ", residual norm = " << rnorm;
1015  if (rnorm > 0.0)
1016  {
1017  pout() << ", rate = " << norm_last/rnorm;
1018  }
1019  pout() << std::endl;
1020  }
1021 
1022  goNorm = rnorm > m_normThresh; //keep iterating if norm is not small enough
1023  goRedu = rnorm > m_eps*initial_rnorm; //keep iterating if initial norm is not reduced enough
1024  goIter = iter < m_iterMax; //keep iterating if iter < max iteration count
1025  goHang = iter < m_imin || rnorm <(1-m_hang)*norm_last;//keep iterating if we didn't hang
1026  }
1027  if ((rnorm > 10.*initial_rnorm) && (rnorm > 10.*m_eps))
1028  {
1029  pout() << "solver seems to have blown up" << endl;
1030  MayDay::Error("kaboom");
1031  }
1032  m_exitStatus = int(!goRedu) + int(!goIter)*2 + int(!goHang)*4 + int(!goNorm)*8;
1033  if (m_verbosity >= 2)
1034  {
1035  pout() << " AMRMultiGrid:: iteration = " << iter << ", residual norm = " << rnorm << std::endl;
1036  }
1037  if (m_verbosity > 1)
1038  {
1039  if (!goIter && goRedu && goNorm) // goRedu=T, goIter=F, goHang=?, goNorm=T
1040  { // m_exitStatus == 0 + 2 + 0|4 + 0 = 2|6
1041  pout() << " AMRMultiGrid:: WARNING: Exit because max iteration count exceeded" << std::endl;
1042  }
1043  if (!goHang && goRedu && goNorm) // goRedu=T, goIter=?, goHang=F, goNorm=T
1044  { // m_exitStatus == 0 + 0|2 + 4 + 0 = 4|6
1045  pout() << " AMRMultiGrid:: WARNING: Exit because of solver hang" << std::endl;
1046  }
1047  if (m_verbosity > 4)
1048  {
1049  pout() << " AMRMultiGrid:: exitStatus = " << m_exitStatus << std::endl;
1050  }
1051  }
1052  for (int i = lowlim; i <= l_max; i++)
1053  {
1054  m_op[i]->clear(*uberCorrection[i]);
1055  delete uberCorrection[i];
1056  }
1057 }
1058 
1059 
1060 template<class T>
1062  int l_max, int l_base)
1063 {
1064  CH_TIME("AMRMultiGrid::relaxOnly");
1065  CH_assert(l_max == l_base);
1066  CH_assert(l_max == 0);
1067  init(a_phi, a_rhs, l_max, l_base);
1068 
1069  CH_assert(a_rhs.size() == a_phi.size());
1070 
1071  Vector<T*> uberResidual(a_rhs.size());
1072 
1073  for (int ilev = 0; ilev < m_op.size(); ilev++)
1074  {
1075  uberResidual[ilev] = new T();
1076  m_op[ilev]->create(*uberResidual[ilev], *a_rhs[ilev]);
1077  }
1078 
1079  //compute initial residual and initialize internal residual to it
1080  m_op[0]->residual(*uberResidual[0], *a_phi[0], *a_rhs[0], true);
1081  Real initial_rnorm = m_op[0]->norm(*uberResidual[0], 0);
1082  Real rnorm = initial_rnorm;
1083  Real norm_last = 2*initial_rnorm;
1084 
1085 
1086  int iter=0;
1087  if (m_verbosity >= 3)
1088  {
1089  pout() << " AMRMultiGrid::relaxOnly iteration = " << iter << ", residual norm = " << rnorm << std::endl;
1090  }
1091  bool goNorm = rnorm > m_normThresh; //iterate if norm is not small enough
1092  bool goRedu = rnorm > m_eps*initial_rnorm; //iterate if initial norm is not reduced enough
1093  bool goIter = iter < m_iterMax; //iterate if iter < max iteration count
1094  bool goHang = iter < m_imin || rnorm <(1-m_hang)*norm_last;//iterate if we didn't hang
1095  while (goIter && goRedu && goHang && goNorm)
1096  {
1097  norm_last = rnorm;
1098  m_op[0]->relax(*a_phi[0], *a_rhs[0], 1);
1099 
1100  iter++;
1101  //recompute residual
1102  m_op[0]->residual(*uberResidual[0], *a_phi[0], *a_rhs[0], true);
1103  rnorm = m_op[0]->norm(*uberResidual[0], 2);
1104  if (m_verbosity >= 4)
1105  {
1106  pout() << " AMRMultiGrid::relaxOnly iteration = " << iter << ", residual norm = " << rnorm
1107  << ", rate = " << norm_last/rnorm << std::endl;
1108  }
1109  goNorm = rnorm > m_normThresh; //keep iterating if norm is not small enough
1110  goRedu = rnorm > m_eps*initial_rnorm; //keep iterating if initial norm is not reduced enough
1111  goIter = iter < m_iterMax; //keep iterating if iter < max iteration count
1112  goHang = iter < m_imin || rnorm <(1-m_hang)*norm_last;//keep iterating if we didn't hang
1113  }
1114  m_exitStatus = 0;
1115  for (int i = 0; i < m_op.size(); i++)
1116  {
1117  m_op[i]->clear(*uberResidual[i]);
1118  delete uberResidual[i];
1119  }
1120 }
1121 
1122 template <class T>
1124 {
1125  for (int i = 0; i < m_op.size(); i++)
1126  {
1127  m_op[i]->clear(*m_correction[i]);
1128  m_op[i]->clear(*m_residual[i]);
1129  m_op[i]->clear(*m_resC[i]);
1130 
1131  delete m_correction[i];
1132  delete m_residual[i];
1133  delete m_resC[i];
1134  delete m_op[i];
1135  delete m_mg[i];
1136 
1137  m_correction[i] = NULL;
1138  m_residual[i] = NULL;
1139  m_resC[i] = NULL;
1140  m_op[i] = NULL;
1141  m_mg[i] = NULL;
1142  }
1143  m_hasInitBeenCalled.resize(0);
1144 }
1145 
1146 template <class T>
1147 void AMRMultiGrid<T>::init(const Vector<T*>& a_phi, const Vector<T*>& a_rhs,
1148  int l_max, int l_base)
1149 {
1150  //if (m_hasInitBeenCalled[l_max]=='t' && m_hasInitBeenCalled[l_base]=='t') return;
1151  CH_TIME("AMRMultiGrid::init");
1152 // CH_assert(a_phi.size()>=m_op.size());
1153 // CH_assert(a_rhs.size()>=m_op.size());
1154  for (int i = l_base; i <= l_max; i++)
1155  {
1156  m_hasInitBeenCalled[i] = 't';
1157  AMRLevelOp<T>& op = *(m_op[i]);
1158 
1159  op.create(*m_correction[i], *a_phi[i]);
1160  op.create(*m_residual[i], *a_rhs[i]);
1161  int r = op.refToCoarser();
1162  if (i!= l_base)
1163  {
1164  int intermediateDepth = 0;
1165 
1166  while (r>2)
1167  {
1168  r/=2;
1169  intermediateDepth++;
1170  }
1171  m_mg[i]->m_depth = intermediateDepth;
1172  }
1173  m_mg[i]->init(*a_phi[i], *a_rhs[i]);
1174 
1175  if (i != l_base)
1176  {
1177  r = op.refToCoarser();
1178  op.createCoarsened(*m_resC[i], *a_rhs[i], r);
1179  op.buildCopier(m_resCopier[i], *a_rhs[i-1], *m_resC[i]);
1180  m_reverseCopier[i] = m_resCopier[i];
1181  m_reverseCopier[i].reverse();
1182  }
1183  }
1184 }
1185 template <class T>
1186 void AMRMultiGrid<T>::revert(const Vector<T*>& a_phi, const Vector<T*>& a_rhs,
1187  int l_max, int l_base)
1188 {
1189  CH_TIME("AMRMultiGrid::revert");
1190  for (int i = l_base; i <= l_max; i++)
1191  {
1192  if (i!= l_base)
1193  {
1194  m_mg[i]->m_depth = m_mg[i]->m_defaultDepth;
1195  }
1196  }
1197 }
1198 
1199 template <class T>
1200 void AMRMultiGrid<T>::setBottomSolver(int l_max, int l_base)
1201 {
1202  // for (int ilev = 0; ilev <= m_mg.size(); ilev++)
1203  for (int ilev = l_base; ilev <= l_max; ilev++)
1204  {
1205  // only do this if we haven't already set these
1206  if (!m_solverParamsSet)
1207  {
1208  m_mg[ilev]->m_pre = m_pre;
1209  m_mg[ilev]->m_post = m_post;
1210  m_mg[ilev]->m_bottom = m_bottom;
1211  }
1212  m_mg[ilev]->m_bottomSolver = &m_nosolve;
1213  }
1214 
1215  m_mg[l_base]->setBottomSolver(m_bottomSolver);
1216 }
1217 
1218 template <class T>
1219 void AMRMultiGrid<T>::setBottomSolverEpsCushion(Real a_bottomSolverEpsCushion)
1220 {
1221  CH_assert((a_bottomSolverEpsCushion > 0.) &&
1222  (a_bottomSolverEpsCushion <= 1.));
1223  m_bottomSolverEpsCushion = a_bottomSolverEpsCushion;
1224 }
1225 
1226 template <class T>
1227 void AMRMultiGrid<T>::define(const ProblemDomain& a_coarseDomain,
1228  AMRLevelOpFactory<T>& a_factory,
1229  LinearSolver<T>* a_bottomSolver,
1230  int a_maxAMRLevels)
1231 {
1232  CH_TIME("AMRMultiGrid::define");
1233  this->clear();
1234  m_op.resize( a_maxAMRLevels, NULL);
1235  m_mg.resize( a_maxAMRLevels, NULL);
1236  m_hasInitBeenCalled.resize(a_maxAMRLevels, 'f');
1237 
1238  m_correction.resize( a_maxAMRLevels, NULL);
1239  m_residual. resize( a_maxAMRLevels, NULL);
1240  m_resC. resize( a_maxAMRLevels, NULL);
1241  m_resCopier. resize( a_maxAMRLevels);
1242  m_reverseCopier.resize( a_maxAMRLevels);
1243  m_bottomSolver = a_bottomSolver;
1244 
1245  ProblemDomain current = a_coarseDomain;
1246  for (int i = 0; i < a_maxAMRLevels; i++)
1247  {
1248 
1249  m_correction[i] = new T();
1250  m_residual[i] = new T();
1251  m_resC[i] = new T();
1252  m_mg[i]= new MultiGrid<T>();
1253  m_op[i] = a_factory.AMRnewOp(current);
1254  m_mg[i]->define(a_factory, &m_nosolve, current, m_maxDepth, m_op[i]);
1255 
1256  // Only do this if it will be used (avoiding a reference to invalid
1257  // and/or unavailable refinement ratios)
1258  if (i < a_maxAMRLevels-1)
1259  {
1260  current.refine(a_factory.refToFiner(current));
1261  }
1262  }
1263 }
1264 
1265 template<class T>
1267  Vector<T*>& a_uberResidual,
1268  int ilev, int l_max, int l_base)
1269 {
1270  if (ilev == l_max)
1271  {
1272  for (int level = l_base; level <= l_max; level++)
1273  {
1274  m_op[level]->assignLocal(*m_residual[level], *a_uberResidual[level]);
1275  m_op[level]->setToZero(*m_correction[level]);
1276  }
1277  }
1278 
1279  if (l_max == l_base)
1280  {
1281  CH_assert(ilev == l_base);
1282  m_mg[l_base]->oneCycle(*(a_uberCorrection[ilev]), *(a_uberResidual[ilev]));
1283  }
1284  else if (ilev == l_base)
1285  {
1286  m_mg[l_base]->oneCycle(*(m_correction[ilev]), *(m_residual[ilev]));
1287 
1288  m_op[ilev]->incr(*(a_uberCorrection[ilev]), *(m_correction[ilev]), 1.0);
1289  }
1290  else
1291  {
1292  //============= Downsweep ========================
1293 
1294  this->relax(*(m_correction[ilev]), *(m_residual[ilev]), ilev, m_pre);
1295  m_op[ilev]->incr(*(a_uberCorrection[ilev]), *(m_correction[ilev]), 1.0);
1296 
1297  // Set next coarser level correction to zero
1298  m_op[ilev-1]->setToZero(*(m_correction[ilev-1]));
1299 
1300  // Recompute residual on next coarser level
1301  // for the valid region NOT covered by this level.
1302  computeAMRResidualLevel(m_residual,
1303  a_uberCorrection,
1304  a_uberResidual,
1305  l_max, l_base, ilev-1,
1306  true);
1307 
1308  // Compute the restriction of the residual to the coarser level resC.
1309  m_op[ilev]->AMRRestrictS(*(m_resC[ilev]),
1310  *(m_residual[ilev]),
1311  *(m_correction[ilev]),
1312  *(m_correction[ilev-1]),
1313  *(a_uberCorrection[ilev]));
1314 
1315  // Overwrite residual on the valid region of the next coarser level
1316  // with coarsened residual from this level
1317  m_op[ilev-1]->assignCopier(*m_residual[ilev-1], *(m_resC[ilev]), m_resCopier[ilev]);
1318 
1319  //============finish Compute residual for the next coarser level======
1320 
1321  for (int img = 0; img < m_numMG; img++)
1322  {
1323  AMRVCycle(a_uberCorrection, a_uberResidual, ilev-1, l_max, l_base);
1324  }
1325 
1326  //================= Upsweep ======================
1327  //increment the correction with coarser version
1328  //m_op[ilev]->AMRProlong(*(m_correction[ilev]), *(m_correction[ilev-1]));
1329  m_op[ilev]->AMRProlongS(*(m_correction[ilev]), *(m_correction[ilev-1]),
1330  *m_resC[ilev], m_reverseCopier[ilev]);
1331  //recompute residual
1332  m_op[ilev]->AMRUpdateResidual(*(m_residual[ilev]), *(m_correction[ilev]), *(m_correction[ilev-1]));
1333 
1334 
1335  //compute correction to the correction
1336  T& dCorr = *(a_uberCorrection[ilev]); // user uberCorrection as holder for correction to correction
1337  m_op[ilev]->setToZero(dCorr);
1338  this->relax(dCorr, *(m_residual[ilev]), ilev, m_post);
1339 
1340  //correct the correction with the correction to the correction
1341  m_op[ilev]->incr(*(m_correction[ilev]), dCorr, 1.0);
1342 
1343 
1344  m_op[ilev]->assignLocal(*(a_uberCorrection[ilev]), *(m_correction[ilev]));
1345 
1346  }
1347 }
1348 
1349 #include "NamespaceFooter.H"
1350 #endif
std::ostream & pout()
Use this in place of std::cout for program output.
int m_bottom
Definition: AMRMultiGrid.H:380
virtual AMRLevelOp< T > * AMRnewOp(const ProblemDomain &a_indexSpace)=0
virtual int refToFiner(const ProblemDomain &a_indexSpace) const =0
#define CH_TIMERS(name)
Definition: CH_Timer.H:70
virtual Real AMRNorm(const T &a_coarResid, const T &a_fineResid, const int &a_refRat, const int &a_ord)
Definition: AMRMultiGrid.H:63
void init(const Vector< T * > &a_phi, const Vector< T * > &a_rhs, int l_max, int l_base)
Definition: AMRMultiGrid.H:1147
bool m_solverParamsSet
Definition: AMRMultiGrid.H:378
A reference-counting handle class.
Definition: RefCountedPtr.H:66
virtual void AMRUpdateResidual(T &a_residual, const T &a_correction, const T &a_coarseCorrection)=0
virtual void AMROperator(T &a_LofPhi, const T &a_phiFine, const T &a_phi, const T &a_phiCoarse, bool a_homogeneousDomBC, AMRLevelOp< T > *a_finerOp)=0
virtual void AMRResidualNC(T &a_residual, const T &a_phiFine, const T &a_phi, const T &a_rhs, bool a_homogeneousBC, AMRLevelOp< T > *a_finerOp)=0
#define CH_assert(cond)
Definition: CHArray.H:37
virtual void buildCopier(Copier &a_copier, const T &a_lhs, const T &a_rhs)
Definition: AMRMultiGrid.H:184
A class to facilitate interaction with physical boundary conditions.
Definition: ProblemDomain.H:130
virtual void define(const ProblemDomain &a_coarseDomain, AMRLevelOpFactory< T > &a_factory, LinearSolver< T > *a_bottomSolver, int a_numLevels)
Definition: AMRMultiGrid.H:1227
virtual ~AMRMultiGridInspector()
Destructor.
Definition: AMRMultiGrid.H:267
int m_iterMax
Definition: AMRMultiGrid.H:379
virtual void AMROperatorNF(T &a_LofPhi, const T &a_phi, const T &a_phiCoarse, bool a_homogeneousBC)=0
LinearSolver< T > * m_bottomSolver
Definition: AMRMultiGrid.H:488
AMRLevelOp< T > & levelOp(int level)
Definition: AMRMultiGrid.H:672
A strange but true thing to make copying from one boxlayoutdata to another fast.
Definition: Copier.H:137
int m_verbosity
Definition: AMRMultiGrid.H:379
void revert(const Vector< T * > &a_phi, const Vector< T * > &a_rhs, int l_max, int l_base)
Definition: AMRMultiGrid.H:1186
virtual void AMRRestrictS(T &a_resCoarse, const T &a_residual, const T &a_correction, const T &a_coarseCorrection, T &scratch)
Definition: AMRMultiGrid.H:209
#define CH_START(tpointer)
Definition: CH_Timer.H:78
virtual void solve(Vector< T * > &a_phi, const Vector< T * > &a_rhs, int l_max, int l_base, bool a_zeroPhi=true, bool forceHomogeneous=false)
Definition: AMRMultiGrid.H:833
void computeAMROperator(Vector< T * > &a_lph, Vector< T * > &a_phi, int l_max, int l_base, bool a_homogeneousBC=false)
Definition: AMRMultiGrid.H:752
Vector< T * > m_residual
Definition: AMRMultiGrid.H:481
AMRMultiGridInspector()
Base class constructor. This must be called by all subclasses.
Definition: AMRMultiGrid.H:262
virtual Real norm(const T &a_rhs, int a_ord)=0
int m_pre
Definition: AMRMultiGrid.H:380
void addInspector(RefCountedPtr< AMRMultiGridInspector< T > > &a_inspector)
Definition: AMRMultiGrid.H:332
Definition: NoOpSolver.H:20
virtual void dumpStuff(Vector< T * > data, string filename)
Definition: AMRMultiGrid.H:54
Vector< Copier > m_reverseCopier
Definition: AMRMultiGrid.H:484
virtual void AMRProlongS(T &a_correction, const T &a_coarseCorrection, T &a_temp, const Copier &a_copier)
Definition: AMRMultiGrid.H:204
virtual void createCoarsened(T &a_lhs, const T &a_rhs, const int &a_refRat)=0
Real m_eps
Definition: AMRMultiGrid.H:377
Definition: AMRMultiGrid.H:306
void clear()
Definition: AMRMultiGrid.H:1123
virtual void outputAMR(Vector< T * > &a_rhs, string &a_name)
Definition: AMRMultiGrid.H:77
Vector< MultiGrid< T > * > m_mg
Definition: AMRMultiGrid.H:479
Real m_convergenceMetric
Definition: AMRMultiGrid.H:390
virtual void AMRResidualNF(T &a_residual, const T &a_phi, const T &a_phiCoarse, const T &a_rhs, bool a_homogeneousBC)=0
const int SpaceDim
Definition: SPACE.H:39
void append(const Vector< T > &invec)
Definition: Vector.H:307
Definition: AMRMultiGrid.H:35
#define MPI_CH_REAL
Definition: REAL.H:34
virtual void assignCopier(T &a_lhs, const T &a_rhs, const Copier &a_copier)
Definition: AMRMultiGrid.H:188
virtual void setToZero(T &a_lhs)=0
#define CH_TIMER(name, tpointer)
Definition: CH_Timer.H:55
virtual void AMRResidual(T &a_residual, const T &a_phiFine, const T &a_phi, const T &a_phiCoarse, const T &a_rhs, bool a_homogeneousDomBC, AMRLevelOp< T > *a_finerOp)=0
int m_numMG
Definition: AMRMultiGrid.H:380
void outputAMR(Vector< T * > &a_data, string &a_name, int a_lmax, int a_lbase)
Definition: AMRMultiGrid.H:512
void push_back(const T &in)
Definition: Vector.H:283
virtual int refToCoarser()=0
#define CH_TIME(name)
Definition: CH_Timer.H:59
virtual void solveNoInitResid(Vector< T * > &a_phi, Vector< T * > &a_finalResid, const Vector< T * > &a_rhs, int l_max, int l_base, bool a_zeroPhi=true, bool forceHomogeneous=false)
use if you want final residual
Definition: AMRMultiGrid.H:872
Definition: MultiGrid.H:259
virtual void solveNoInit(Vector< T * > &a_phi, const Vector< T * > &a_rhs, int l_max, int l_base, bool a_zeroPhi=true, bool forceHomogeneous=false)
Definition: AMRMultiGrid.H:845
const char * name(const FArrayBox &a_dummySpecializationArg)
Definition: CH_HDF5.H:741
Real m_bottomSolverEpsCushion
Definition: AMRMultiGrid.H:393
int m_maxDepth
max no. of coarsenings – -1 (default) means coarsen as far as possible
Definition: AMRMultiGrid.H:384
double Real
Definition: REAL.H:33
Definition: MultiGrid.H:30
virtual void dumpAMR(Vector< T * > &a_data, string name)
Definition: AMRMultiGrid.H:40
AMRMultiGrid()
Definition: AMRMultiGrid.H:596
Real m_hang
Definition: AMRMultiGrid.H:377
Vector< RefCountedPtr< AMRMultiGridInspector< T > > > m_inspectors
Definition: AMRMultiGrid.H:497
virtual void AMRRestrict(T &a_resCoarse, const T &a_residual, const T &a_correction, const T &a_coarseCorrection)=0
virtual unsigned int orderOfAccuracy(void) const
Definition: AMRMultiGrid.H:215
Vector< T * > m_correction
Definition: AMRMultiGrid.H:480
Vector< T * > m_resC
Definition: AMRMultiGrid.H:482
Real computeAMRResidual(Vector< T * > &a_resid, Vector< T * > &a_phi, const Vector< T * > &a_rhs, int l_max, int l_base, bool a_homogeneousBC=false, bool a_computeNorm=true)
Definition: AMRMultiGrid.H:679
AMRLevelOp()
Constructor.
Definition: AMRMultiGrid.H:49
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.
AMRMultiGridInspector & operator=(const AMRMultiGridInspector &)
virtual void enforceCFConsistency(T &a_coarseCorrection, const T &a_correction)
This routine is for operators with orderOfAccuracy()>2.
Definition: AMRMultiGrid.H:221
Definition: MultiGrid.H:296
Definition: AMRMultiGrid.H:257
void relaxOnlyHomogeneous(Vector< T * > &a_phi, const Vector< T * > &a_rhs, int l_max, int l_base)
Definition: AMRMultiGrid.H:1061
virtual void outputLevel(T &a_rhs, string &a_name)
Definition: AMRMultiGrid.H:72
int m_exitStatus
Definition: AMRMultiGrid.H:379
ProblemDomain & refine(int a_refinement_ratio)
Refine this problem domain.
virtual void AMROperatorNC(T &a_LofPhi, const T &a_phiFine, const T &a_phi, bool a_homogeneousBC, AMRLevelOp< T > *a_finerOp)=0
virtual void create(T &a_lhs, const T &a_rhs)=0
virtual void AMRProlong(T &a_correction, const T &a_coarseCorrection)=0
Definition: LinearSolver.H:158
void computeAMRResidualLevel(Vector< T * > &a_resid, Vector< T * > &a_phi, const Vector< T * > &a_rhs, int l_max, int l_base, int ilev, bool a_homogeneousBC)
Definition: AMRMultiGrid.H:775
void relax(T &phi, T &R, int depth, int nRelax=2)
Definition: AMRMultiGrid.H:629
int m_post
Definition: AMRMultiGrid.H:380
int m_imin
Definition: AMRMultiGrid.H:379
virtual void recordCorrections(const Vector< T * > &a_corrections, int a_minLevel, int a_maxLevel, int a_iter)=0
Vector< MGLevelOp< T > * > getOperatorsOp()
Definition: AMRMultiGrid.H:571
#define CH_STOP(tpointer)
Definition: CH_Timer.H:80
Vector< char > m_hasInitBeenCalled
Definition: AMRMultiGrid.H:490
void setSolverParameters(const int &a_pre, const int &a_post, const int &a_bottom, const int &a_numMG, const int &a_iterMax, const Real &a_eps, const Real &a_hang, const Real &a_normThresh)
Definition: AMRMultiGrid.H:524
AMRMultiGrid & operator=(const AMRMultiGrid< T > &)
size_t size() const
Definition: Vector.H:177
Vector< Vector< MGLevelOp< T > * > > getOperatorsMG()
Definition: AMRMultiGrid.H:584
T Max(const T &a_a, const T &a_b)
Definition: Misc.H:39
virtual void assign(T &a_lhs, const T &a_rhs)=0
Vector< Copier > m_resCopier
Definition: AMRMultiGrid.H:483
virtual Real localMaxNorm(const T &a_phi)
Definition: AMRMultiGrid.H:198
virtual void AMRVCycle(Vector< T * > &a_correction, Vector< T * > &a_residual, int l, int l_max, int l_base)
Definition: AMRMultiGrid.H:1266
virtual void dumpLevel(T &a_data, string name)
Definition: AMRMultiGrid.H:44
Definition: AMRMultiGrid.H:231
NoOpSolver< T > m_nosolve
Definition: AMRMultiGrid.H:486
virtual void recordResiduals(const Vector< T * > &a_residuals, int a_minLevel, int a_maxLevel, int a_iter)=0
Vector< MGLevelOp< T > * > getAllOperators()
Definition: AMRMultiGrid.H:552
Real m_normThresh
Definition: AMRMultiGrid.H:377
virtual ~AMRLevelOpFactory()
Definition: AMRMultiGrid.H:234
Vector< AMRLevelOp< T > * > & getAMROperators()
Definition: AMRMultiGrid.H:431
virtual ~AMRLevelOp()
Destructor.
Definition: AMRMultiGrid.H:59
virtual void zeroCovered(T &a_lhs, T &a_rhs, const Copier &a_copier)
Definition: AMRMultiGrid.H:193
void setMGCycle(int a_numMG)
Definition: AMRMultiGrid.H:619
void setBottomSolverEpsCushion(Real a_bottomSolverEpsCushion)
Definition: AMRMultiGrid.H:1219
Vector< AMRLevelOp< T > * > m_op
Definition: AMRMultiGrid.H:478
virtual ~AMRMultiGrid()
Definition: AMRMultiGrid.H:665
void setBottomSolver(int l_max, int l_base)
set up bottom solver for internal MG solver
Definition: AMRMultiGrid.H:1200