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