Chombo + EB  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  // used to give an additional cushion in the EPS used for bottom solves
395 
396  ///
397  /**
398  resid = L(phi) - rhs
399  */
401  Vector<T*>& a_phi,
402  const Vector<T*>& a_rhs,
403  int l_max,
404  int l_base,
405  bool a_homogeneousBC=false,
406  bool a_computeNorm=true);
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  /// set up bottom solver for internal MG solver
458  /**
459  This function is normally called by the solve(...) function.
460  However, it must be called if solve will not be called (in particular,
461  if only the V-cycle is being used)
462  */
463  void setBottomSolver(int l_max, int l_base);
464 
465  void setBottomSolverEpsCushion(Real a_bottomSolverEpsCushion);
466 
467  /// Dump out the state of the solver. AMR levels, MG levels, bottom solver, box counts, etc
468  void getInfo() const;
469 protected:
470 
471  void relax(T& phi, T& R, int depth, int nRelax = 2);
472 
473  void computeAMRResidualLevel(Vector<T*>& a_resid,
474  Vector<T*>& a_phi,
475  const Vector<T*>& a_rhs,
476  int l_max, int l_base, int ilev,
477  bool a_homogeneousBC);
478 
486 
488 
490 
492 
493  void clear();
494 
495 private:
496 
497  // A list of inspectors maintained by this instance.
499 
500  // Forbidden copiers.
503 };
504 
505 //*******************************************************
506 // AMRMultigrid Implementation
507 //*******************************************************
508 
509 //===================================================================
510 
511 template <class T>
512 void
513 AMRMultiGrid<T>::outputAMR(Vector<T*>& a_data, string& a_name, int a_lmax, int a_lbase)
514 {
515  Vector<T*> outputData;
516  for (int ilev = a_lbase; ilev <= a_lmax; ilev++)
517  {
518  outputData.push_back(a_data[ilev]);
519  }
520  m_op[a_lbase]->outputAMR(outputData, a_name);
521 }
522 
523 template <class T>
524 void
526  const int& a_post,
527  const int& a_bottom,
528  const int& a_numMG,
529  const int& a_iterMax,
530  const Real& a_eps,
531  const Real& a_hang,
532  const Real& a_normThresh)
533 {
534  m_solverParamsSet = true;
535  m_pre = a_pre;
536  m_post = a_post;
537  m_bottom = a_bottom;
538  m_eps = a_eps;
539  m_hang = a_hang;
540  m_normThresh = a_normThresh;
541  m_iterMax = a_iterMax;
542  for (int img = 0; img < m_mg.size(); img++)
543  {
544  m_mg[img]->m_pre = a_pre;
545  m_mg[img]->m_post = a_post;
546  m_mg[img]->m_bottom = a_bottom;
547  }
548  setMGCycle(a_numMG);
549  m_bottomSolverEpsCushion = 1.0;
550 }
551 template <class T>
554 {
555  Vector< MGLevelOp<T>* > retval;
556  for (int iop = 0; iop < m_op.size(); iop++)
557  {
558  MGLevelOp<T>* operPtr = (MGLevelOp<T>*) m_op[iop];
559  retval.push_back(operPtr);
560  }
561 
562  for (int img = 0; img < m_mg.size(); img++)
563  {
564  Vector< MGLevelOp<T>* > mgOps = m_mg[img]->getAllOperators();
565  retval.append(mgOps);
566  }
567  return retval;
568 }
569 
570 template <class T>
573 {
574  Vector< MGLevelOp<T>* > retval;
575  for (int iop = 0; iop < m_op.size(); iop++)
576  {
577  MGLevelOp<T>* operPtr = (MGLevelOp<T>*) m_op[iop];
578  retval.push_back(operPtr);
579  }
580  return retval;
581 }
582 
583 template <class T>
586 {
587  Vector< Vector< MGLevelOp<T>* > > retval(m_mg.size());
588 
589  for (int img = 0; img < m_mg.size(); img++)
590  {
591  retval[img] = m_mg[img]->getAllOperators();
592  }
593  return retval;
594 }
595 
596 template <class T>
598  :m_eps(1E-6),
599  m_hang(1E-15),
600  m_normThresh(1E-30),
601  m_imin(5),
602  m_iterMax(20),
603  m_iterMin(-1),
604  m_verbosity(3),
605  m_pre(2),
606  m_post(2),
607  m_bottom(2),
608  m_numMG(1),
609  m_maxDepth(-1),
610  m_convergenceMetric(0.),
611  m_bottomSolverEpsCushion(1.0),
612  m_bottomSolver(NULL),
613  m_inspectors()
614 {
615  m_solverParamsSet = false;
616 
617  //m_maxDepth = 4;
618  //pout() << " AMRMultiGrid<T>::AMRMultiGrid() = " << m_maxDepth << endl;
619 }
620 template <class T>
622 {
623  for (int ilev = 0; ilev < m_op.size(); ilev++)
624  {
625  m_mg[ilev]->m_numMG = a_numMG;
626  m_mg[ilev]->m_cycle = a_numMG;
627  }
628  m_numMG = a_numMG;
629 }
630 template <class T>
631 void AMRMultiGrid<T>::relax(T& a_correction, T& a_residual, int depth, int a_numSmooth)
632 {
633  CH_TIME("AMRMultiGrid::relax");
634 
635  if (m_op[depth]->refToCoarser() > 2)
636  {
637  // intermediate multigrid levels exist between this level and the
638  // next coarser AMR level, so do a mini v-cyle to smooth on those
639 
640  int intermediateDepth = 0;
641  int r = m_op[depth]->refToCoarser();
642  while (r>2)
643  {
644  r/=2;
645  intermediateDepth++;
646  }
647  int tmp = m_mg[depth]->m_depth;
648  m_mg[depth]->m_depth = intermediateDepth;
649  //ok, there is an intermediate multigrid level that is
650  // not an AMR level. We use regular MultiGrid smoothing algorithm for this.
651  // note that MultiGrid::cycle handles relaxation on this level
652  // as well, so there's no need to call the LinearOp's relax function.
653  m_mg[depth]->cycle(0, a_correction, a_residual);
654  m_mg[depth]->m_depth = tmp;
655  }
656  else
657  {
658  // no intermediate multigrid levels between AMR levels, so
659  // just call relax on this level and be done with it.
660  m_op[depth]->relax(a_correction, a_residual, a_numSmooth); //numSmoothDown
661  }
662 
663 }
664 /************/
665 template <class T>
667 {
668  CH_TIME("~AMRMultiGrid");
669  clear();
670 }
671 
672 template <class T>
674 {
675  return *(m_op[level]);
676 }
677 
678 /************/
679 template <class T>
681  Vector<T*>& a_phi,
682  const Vector<T*>& a_rhs,
683  int l_max,
684  int l_base,
685  bool a_homogeneousBC,
686  bool a_computeNorm)
687 {
688  CH_TIME("AMRMultiGrid::computeAMRResidual");
689 
690  Real rnorm = 0;
691  Real localNorm = 0;
692  for (int ilev = l_base; ilev <= l_max; ilev++)
693  {
694  //always used at top level where bcs are inhomogeneous
695  computeAMRResidualLevel(a_resid,
696  a_phi,
697  a_rhs,
698  l_max, l_base, ilev, a_homogeneousBC);
699  if (a_computeNorm)
700  {
701  if (ilev == l_max)
702  {
703  localNorm =m_op[ilev]->localMaxNorm(*a_resid[ilev]);
704  }
705  else
706  {
707  m_op[ilev]->zeroCovered(*a_resid[ilev], *m_resC[ilev+1], m_resCopier[ilev+1]);
708  localNorm = m_op[ilev]->localMaxNorm(*a_resid[ilev]);
709  }
710  const int normType = 0;
711  if (normType==2)
712  localNorm = m_op[ilev]->norm(*a_resid[ilev],normType);
713  rnorm = Max(localNorm, rnorm);
714  }
715  }
716 #ifdef CH_MPI
717  if (a_computeNorm)
718  {
719  CH_TIME("MPI_Allreduce");
720  Real recv;
721  int result = MPI_Allreduce(&rnorm, &recv, 1, MPI_CH_REAL,
722  MPI_MAX, Chombo_MPI::comm);
723  if (result != MPI_SUCCESS)
724  {
725  //bark!!!
726  MayDay::Error("sorry, but I had a communcation error on norm");
727  }
728  rnorm = recv;
729  }
730 #endif
731 
732  return rnorm; // if a_computeNorm is false, then this just returns zero.
733 }
734 
735 template <class T>
737  const Vector<T*>& a_rhs,
738  int l_max,
739  int l_min)
740 {
741  return computeAMRResidual(m_residual,
742  a_phi,
743  a_rhs,
744  l_max,
745  l_min,
746  false,
747  true);
748 }
749 
750 /************/
751 // old inefficient way of doing this
752 /*
753 template <class T>
754 void AMRMultiGrid<T>::computeAMROperator(Vector<T*>& a_lph,
755  Vector<T*>& a_phi,
756  int l_max,
757  int l_base,
758  bool a_homogeneousBC)
759 {
760 
761  for (int ilev = l_base; ilev <= l_max; ilev++)
762  {
763  m_op[ilev]->create(*m_residual[ilev], *a_lph[ilev]);
764  m_op[ilev]->setToZero(*m_residual[ilev]);
765  }
766  computeAMRResidual(a_lph, a_phi, m_residual, l_max, l_base, a_homogeneousBC, false);
767  // fixing the bug of not negating the result --- Qinghai Zhang 12/10/2009
768  // In loving memory of the two weeks I lost over this bug!!!
769  for (int ilev = l_base; ilev <= l_max; ilev++)
770  {
771  m_op[ilev]->scale(*a_lph[ilev], -1.);
772  }
773 }
774 */
775 template <class T>
777  Vector<T*>& a_phi,
778  int l_max,
779  int l_base,
780  bool a_homogeneousBC)
781 {
782  CH_TIME("AMRMultiGrid<T>::computeAMROperator");
783 
784  for(int ilev=l_base; ilev<=l_max; ilev++)
785  {
786  if (l_max != l_base)
787  {
788  if (ilev == l_max)
789  {
790  m_op[l_max]->AMROperatorNF(*(a_lphi[l_max]),*(a_phi[l_max]),
791  *(a_phi[l_max-1]), a_homogeneousBC);
792  }
793  else if (ilev == l_base)
794  {
795  if (l_base == 0)
796  {
797  m_op[l_base]->AMROperatorNC(*(a_lphi[l_base]), *(a_phi[l_base+1]),
798  *(a_phi[l_base]),
799  a_homogeneousBC, m_op[l_base+1]);
800  }
801  else
802  {
803  m_op[l_base]->AMROperator(*a_lphi[l_base], *a_phi[l_base+1], *a_phi[l_base],
804  *a_phi[l_base-1],
805  a_homogeneousBC, m_op[l_base+1]);
806  }
807  }
808  else
809  {
810  m_op[ilev]->AMROperator(*a_lphi[ilev], *a_phi[ilev+1], *a_phi[ilev],
811  *a_phi[ilev-1],
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]->applyOp(*a_lphi[l_max], *a_phi[l_max], a_homogeneousBC);
821  }
822  else
823  {
824  m_op[l_max]->AMROperatorNF(*(a_lphi[l_max]), *(a_phi[l_max]),
825  *(a_phi[l_max-1]),
826  a_homogeneousBC);
827  }
828  }
829  }
830 }
831 
832 /************/
833 template <class T>
835  Vector<T*>& a_phi,
836  const Vector<T*>& a_rhs,
837  int l_max, int l_base, int ilev,
838  bool a_homogeneousBC)
839 {
840  CH_TIME("AMRMultiGrid<T>::computeAMRResidualLevel");
841  CH_assert(m_hasInitBeenCalled[ilev]=='t');
842 
843  //m_op[ilev]->setToZero(*(a_resid[l_max]));
844  if (l_max != l_base)
845  {
846  if (ilev == l_max)
847  {
848  m_op[l_max]->AMRResidualNF(*(a_resid[l_max]), *(a_phi[l_max]),
849  *(a_phi[l_max-1]), *(a_rhs[l_max]),
850  a_homogeneousBC);
851  }
852  else if (ilev == l_base)
853  {
854  if (l_base == 0)
855  {
856  m_op[l_base]->AMRResidualNC(*(a_resid[l_base]), *(a_phi[l_base+1]),
857  *(a_phi[l_base]), *(a_rhs[l_base]),
858  a_homogeneousBC, m_op[l_base+1]);
859  }
860  else
861  {
862  m_op[l_base]->AMRResidual(*a_resid[l_base], *a_phi[l_base+1], *a_phi[l_base],
863  *a_phi[l_base-1], *a_rhs[l_base],
864  a_homogeneousBC, m_op[l_base+1]);
865  }
866  }
867  else
868  {
869  m_op[ilev]->AMRResidual(*a_resid[ilev], *a_phi[ilev+1], *a_phi[ilev],
870  *a_phi[ilev-1], *a_rhs[ilev],
871  a_homogeneousBC, m_op[ilev+1]);
872  }
873  }
874  else
875  {
876  CH_assert(ilev == l_base);
877  if (l_base == 0)
878  {
879  m_op[l_max]->residual(*a_resid[l_max], *a_phi[l_max], *a_rhs[l_max],a_homogeneousBC);
880  }
881  else
882  {
883  m_op[l_max]->AMRResidualNF(*(a_resid[l_max]), *(a_phi[l_max]),
884  *(a_phi[l_max-1]), *(a_rhs[l_max]),
885  a_homogeneousBC);
886  }
887  }
888 
889 }
890 
891 template<class T>
892 void AMRMultiGrid<T>::solve(Vector<T*>& a_phi, const Vector<T*>& a_rhs,
893  int l_max, int l_base, bool a_zeroPhi,
894  bool a_forceHomogeneous)
895 {
896  CH_TIME("AMRMultiGrid::solve");
897  init(a_phi, a_rhs, l_max, l_base);
898  solveNoInit(a_phi, a_rhs, l_max, l_base, a_zeroPhi, a_forceHomogeneous);
899  //puts multigrid depth back where it started so the solver can be reused
900  revert(a_phi, a_rhs, l_max, l_base);
901 }
902 
903 template<class T>
905  const Vector<T*>& a_rhs,
906  int l_max, int l_base, bool a_zeroPhi,
907  bool a_forceHomogeneous)
908 {
909  CH_TIME("AMRMultiGrid::solveNoInit");
910  Vector<T*> uberResidual(a_rhs.size());
911  int lowlim = l_base;
912  if (l_base > 0) // we need an ubercorrection one level lower than l_base
913  lowlim--;
914  // for (int ilev = l_base; ilev <= l_max; ilev++)
915  for (int ilev = lowlim; ilev <= l_max; ilev++)
916  {
917  uberResidual[ilev] = new T();
918  if (ilev >= l_base)
919  {
920  m_op[ilev]->create(*uberResidual[ilev], *a_rhs[ilev]);
921  }
922  }
923  solveNoInitResid(a_phi, uberResidual, a_rhs, l_max, l_base, a_zeroPhi, a_forceHomogeneous);
924  for (int i = lowlim; i <= l_max; i++)
925  {
926  m_op[i]->clear(*uberResidual[i]);
927  delete uberResidual[i];
928  }
929 }
930 
931 template<class T>
933  const Vector<T*>& a_rhs,
934  int l_max, int l_base, bool a_zeroPhi,
935  bool a_forceHomogeneous)
936 {
937  CH_TIMERS("AMRMultiGrid::solveNo-InitResid");
938  CH_TIMER("AMRMultiGrid::AMRVcycle", vtimer);
939  setBottomSolver(l_max, l_base);
940 
941  CH_assert(l_base <= l_max);
942  CH_assert(a_rhs.size() == a_phi.size());
943 
944  //these correspond to the residual and correction
945  //that live in AMRSolver
946  Vector<T*> uberCorrection(a_rhs.size());
947 
948  int lowlim = l_base;
949  if (l_base > 0) // we need an ubercorrection one level lower than l_base
950  lowlim--;
951  // for (int ilev = l_base; ilev <= l_max; ilev++)
952  bool outputIntermediates = false;
953 
954  for (int ilev = lowlim; ilev <= l_max; ilev++)
955  { CH_TIME("uberCorrection");
956  uberCorrection[ilev] = new T();
957  m_op[ilev]->create(*uberCorrection[ilev], *a_phi[ilev]);
958  if (ilev >= l_base)
959  {
960  m_op[ilev]->create(*uberResidual[ilev], *a_rhs[ilev]);
961  m_op[ilev]->setToZero(*(uberResidual[ilev]));
962  }
963  m_op[ilev]->setToZero(*(uberCorrection[ilev]));
964  }
965  // m_op[0]->dumpStuff(uberResidual, string("initialRes.hdf5"));
966  if (a_zeroPhi)
967  for (int ilev = l_base; ilev <=l_max; ++ilev)
968  {
969  m_op[ilev]->setToZero(*(a_phi[ilev]));
970  }
971  //compute initial residual and initialize internal residual to it
972 
973  Real initial_rnorm = 0;
974  {
975  CH_TIME("Initial AMR Residual");
976  initial_rnorm=computeAMRResidual(uberResidual, a_phi, a_rhs, l_max, l_base, a_forceHomogeneous, true);
977  }
978 
979  if (m_convergenceMetric != 0.)
980  {
981  initial_rnorm = m_convergenceMetric;
982  }
983 
984  Real rnorm = initial_rnorm;
985  Real norm_last = 2*initial_rnorm;
986 
987  /// set bottom solver convergence norm and solver tolerance
988  m_bottomSolver->setConvergenceMetrics(initial_rnorm, m_bottomSolverEpsCushion*m_eps);
989 
990  int iter=0;
991  if (m_verbosity >= 2) ////
992  {
993  pout() << setw(12)
994  << setprecision(6)
995  << setiosflags(ios::showpoint)
996  << setiosflags(ios::scientific);
997  pout() << " AMRMultiGrid:: iteration = " << iter << ", residual norm = " << rnorm << std::endl;
998  }
999 
1000  bool goNorm = rnorm > m_normThresh; //iterate if norm is not small enough
1001  bool goRedu = rnorm > m_eps*initial_rnorm; //iterate if initial norm is not reduced enough
1002  bool goIter = iter < m_iterMax; //iterate if iter < max iteration count
1003  bool goHang = iter < m_imin || rnorm <(1-m_hang)*norm_last;//iterate if we didn't hang
1004  bool goMin = iter < m_iterMin ; // iterate if iter < min
1005  while (goMin || (goIter && goRedu && goHang && goNorm))
1006  {
1007 
1008  // Report the residuals to the inspectors.
1009  for (int i = 0; i < m_inspectors.size(); ++i)
1010  m_inspectors[i]->recordResiduals(uberResidual, l_base, l_max, iter);
1011 
1012  if (outputIntermediates)
1013  {
1014  char strresname[100];
1015  sprintf(strresname, "amrmg.res.iter.%03d", iter);
1016  string nameres(strresname);
1017  outputAMR(uberResidual, nameres, l_max, l_base);
1018  }
1019 
1020  norm_last = rnorm;
1021 
1022  //this generates a correction from the current residual
1023  CH_START(vtimer);
1024  AMRVCycle(uberCorrection, uberResidual, l_max, l_max, l_base);
1025  CH_STOP(vtimer);
1026 
1027  char charname[100];
1028  sprintf(charname, "resid_iter%03d.%dd.hdf5", iter, SpaceDim);
1029  string sname(charname);
1030  // m_op[0]->dumpAMR(uberResidual, sname);
1031  // Report the corrections to the inspectors.
1032  for (int i = 0; i < m_inspectors.size(); ++i)
1033  m_inspectors[i]->recordCorrections(uberCorrection, l_base, l_max, iter);
1034 
1035  //increment phi by correction and reset correction to zero
1036  for (int ilev = l_base; ilev <= l_max; ilev++)
1037  {
1038  if (outputIntermediates)
1039  {
1040  char strcorname[100];
1041  sprintf(strcorname, "amrmg.phi.iter.%03d", iter);
1042  string namecor(strcorname);
1043  outputAMR(a_phi, namecor, l_max, l_base);
1044  }
1045 
1046  m_op[ilev]->incr(*(a_phi[ilev]), *(uberCorrection[ilev]), 1.0);
1047 
1048  if (outputIntermediates)
1049  {
1050  char strcorname[100];
1051  sprintf(strcorname, "amrmg.cor.iter.%03d", iter);
1052  string namecor(strcorname);
1053  outputAMR(uberCorrection, namecor, l_max, l_base);
1054  }
1055 
1056  m_op[ilev]->setToZero(*(uberCorrection[ilev]));
1057  // clean const out
1058  }
1059 
1060  // For solvers with accuracy higher than 2nd order
1061  // consistency between levels has to be explicitly enforced.
1062  // Qinghai Zhang
1063  if (m_op[0]->orderOfAccuracy()>2)
1064  {
1065  for (int ilev=l_max; ilev>l_base; ilev--)
1066  {
1067  m_op[ilev]->enforceCFConsistency(*a_phi[ilev-1], *a_phi[ilev]);
1068  }
1069  }
1070  //------end enforcing consistency.
1071 
1072  //recompute residual
1073  rnorm = computeAMRResidual(uberResidual, a_phi, a_rhs, l_max, l_base, a_forceHomogeneous, true);
1074  iter++;
1075  if (m_verbosity >= 3) ////
1076  {
1077  pout() << " AMRMultiGrid:: iteration = " << iter << ", residual norm = " << rnorm;
1078  if (rnorm > 0.0)
1079  {
1080  pout() << ", rate = " << norm_last/rnorm;
1081  }
1082  pout() << std::endl;
1083  }
1084 
1085  goNorm = rnorm > m_normThresh; //keep iterating if norm is not small enough
1086  goRedu = rnorm > m_eps*initial_rnorm; //keep iterating if initial norm is not reduced enough
1087  goIter = iter < m_iterMax; //keep iterating if iter < max iteration count
1088  goHang = iter < m_imin || rnorm <(1-m_hang)*norm_last;//keep iterating if we didn't hang
1089  goMin = iter < m_iterMin ; // keep iterating if iter < min
1090  }
1091  // if ((rnorm > 10.*initial_rnorm) && (rnorm > 10.*m_eps))
1092  // {
1093  // pout() << "solver seems to have blown up" << endl;
1094  // MayDay::Error("kaboom");
1095  // }
1096  m_exitStatus = int(!goRedu) + int(!goIter)*2 + int(!goHang)*4 + int(!goNorm)*8;
1097  if (m_verbosity >= 2)
1098  {
1099  pout() << " AMRMultiGrid:: iteration = " << iter << ", residual norm = " << rnorm << std::endl;
1100  }
1101  if (m_verbosity > 1)
1102  {
1103  if (!goIter && goRedu && goNorm) // goRedu=T, goIter=F, goHang=?, goNorm=T
1104  { // m_exitStatus == 0 + 2 + 0|4 + 0 = 2|6
1105  pout() << " AMRMultiGrid:: WARNING: Exit because max iteration count exceeded" << std::endl;
1106  }
1107  if (!goHang && goRedu && goNorm) // goRedu=T, goIter=?, goHang=F, goNorm=T
1108  { // m_exitStatus == 0 + 0|2 + 4 + 0 = 4|6
1109  pout() << " AMRMultiGrid:: WARNING: Exit because of solver hang" << std::endl;
1110  }
1111  if (m_verbosity > 4)
1112  {
1113  pout() << " AMRMultiGrid:: exitStatus = " << m_exitStatus << std::endl;
1114  }
1115  }
1116  for (int i = lowlim; i <= l_max; i++)
1117  {
1118  m_op[i]->clear(*uberCorrection[i]);
1119  delete uberCorrection[i];
1120  }
1121 }
1122 
1123 template<class T>
1125  int l_max, int l_base)
1126 {
1127  CH_TIME("AMRMultiGrid::relaxOnly");
1128  CH_assert(l_max == l_base);
1129  CH_assert(l_max == 0);
1130  init(a_phi, a_rhs, l_max, l_base);
1131 
1132  CH_assert(a_rhs.size() == a_phi.size());
1133 
1134  Vector<T*> uberResidual(a_rhs.size());
1135 
1136  for (int ilev = 0; ilev < m_op.size(); ilev++)
1137  {
1138  uberResidual[ilev] = new T();
1139  m_op[ilev]->create(*uberResidual[ilev], *a_rhs[ilev]);
1140  }
1141 
1142  //compute initial residual and initialize internal residual to it
1143  m_op[0]->residual(*uberResidual[0], *a_phi[0], *a_rhs[0], true);
1144  Real initial_rnorm = m_op[0]->norm(*uberResidual[0], 0);
1145  Real rnorm = initial_rnorm;
1146  Real norm_last = 2*initial_rnorm;
1147 
1148  int iter=0;
1149  if (m_verbosity >= 3)
1150  {
1151  pout() << " AMRMultiGrid::relaxOnly iteration = " << iter << ", residual norm = " << rnorm << std::endl;
1152  }
1153  bool goNorm = rnorm > m_normThresh; //iterate if norm is not small enough
1154  bool goRedu = rnorm > m_eps*initial_rnorm; //iterate if initial norm is not reduced enough
1155  bool goIter = iter < m_iterMax; //iterate if iter < max iteration count
1156  bool goHang = iter < m_imin || rnorm <(1-m_hang)*norm_last;//iterate if we didn't hang
1157  while (goIter && goRedu && goHang && goNorm)
1158  {
1159  norm_last = rnorm;
1160  m_op[0]->relax(*a_phi[0], *a_rhs[0], 1);
1161 
1162  iter++;
1163  //recompute residual
1164  m_op[0]->residual(*uberResidual[0], *a_phi[0], *a_rhs[0], true);
1165  rnorm = m_op[0]->norm(*uberResidual[0], 2);
1166  if (m_verbosity >= 4)
1167  {
1168  pout() << " AMRMultiGrid::relaxOnly iteration = " << iter << ", residual norm = " << rnorm
1169  << ", rate = " << norm_last/rnorm << std::endl;
1170  }
1171  goNorm = rnorm > m_normThresh; //keep iterating if norm is not small enough
1172  goRedu = rnorm > m_eps*initial_rnorm; //keep iterating if initial norm is not reduced enough
1173  goIter = iter < m_iterMax; //keep iterating if iter < max iteration count
1174  goHang = iter < m_imin || rnorm <(1-m_hang)*norm_last;//keep iterating if we didn't hang
1175  }
1176  m_exitStatus = 0;
1177  for (int i = 0; i < m_op.size(); i++)
1178  {
1179  m_op[i]->clear(*uberResidual[i]);
1180  delete uberResidual[i];
1181  }
1182 }
1183 
1184 template <class T>
1186 {
1187  for (int i = 0; i < m_op.size(); i++)
1188  {
1189  m_op[i]->clear(*m_correction[i]);
1190  m_op[i]->clear(*m_residual[i]);
1191  m_op[i]->clear(*m_resC[i]);
1192 
1193  delete m_correction[i];
1194  delete m_residual[i];
1195  delete m_resC[i];
1196  delete m_op[i];
1197  delete m_mg[i];
1198 
1199  m_correction[i] = NULL;
1200  m_residual[i] = NULL;
1201  m_resC[i] = NULL;
1202  m_op[i] = NULL;
1203  m_mg[i] = NULL;
1204  }
1205  m_hasInitBeenCalled.resize(0);
1206 }
1207 
1208 template <class T>
1209 void AMRMultiGrid<T>::init(const Vector<T*>& a_phi, const Vector<T*>& a_rhs,
1210  int l_max, int l_base)
1211 {
1212  //if (m_hasInitBeenCalled[l_max]=='t' && m_hasInitBeenCalled[l_base]=='t') return;
1213  CH_TIME("AMRMultiGrid::init");
1214 // CH_assert(a_phi.size()>=m_op.size());
1215 // CH_assert(a_rhs.size()>=m_op.size());
1216  for (int i = l_base; i <= l_max; i++)
1217  {
1218  m_hasInitBeenCalled[i] = 't';
1219  AMRLevelOp<T>& op = *(m_op[i]);
1220 
1221  op.create(*m_correction[i], *a_phi[i]);
1222  op.create(*m_residual[i], *a_rhs[i]);
1223  int r = op.refToCoarser();
1224  if (i!= l_base)
1225  {
1226  int intermediateDepth = 0;
1227 
1228  while (r>2)
1229  {
1230  r/=2;
1231  intermediateDepth++;
1232  }
1233  m_mg[i]->m_depth = intermediateDepth;
1234  }
1235  m_mg[i]->init(*a_phi[i], *a_rhs[i]);
1236 
1237  if (i != l_base)
1238  {
1239  r = op.refToCoarser();
1240  op.createCoarsened(*m_resC[i], *a_rhs[i], r);
1241  op.buildCopier(m_resCopier[i], *a_rhs[i-1], *m_resC[i]);
1242  m_reverseCopier[i] = m_resCopier[i];
1243  m_reverseCopier[i].reverse();
1244  }
1245  }
1246 }
1247 
1248 template <class T>
1249 void AMRMultiGrid<T>::revert(const Vector<T*>& a_phi, const Vector<T*>& a_rhs,
1250  int l_max, int l_base)
1251 {
1252  CH_TIME("AMRMultiGrid::revert");
1253  for (int i = l_base; i <= l_max; i++)
1254  {
1255  if (i!= l_base)
1256  {
1257  m_mg[i]->m_depth = m_mg[i]->m_defaultDepth;
1258  }
1259  }
1260 }
1261 
1262 template <class T>
1263 void AMRMultiGrid<T>::setBottomSolver(int l_max, int l_base)
1264 {
1265  // for (int ilev = 0; ilev <= m_mg.size(); ilev++)
1266  for (int ilev = l_base; ilev <= l_max; ilev++)
1267  {
1268  // only do this if we haven't already set these
1269  if (!m_solverParamsSet)
1270  {
1271  m_mg[ilev]->m_pre = m_pre;
1272  m_mg[ilev]->m_post = m_post;
1273  m_mg[ilev]->m_bottom = m_bottom;
1274  }
1275  m_mg[ilev]->m_bottomSolver = &m_nosolve;
1276  }
1277 
1278  m_mg[l_base]->setBottomSolver(m_bottomSolver);
1279 }
1280 
1281 template <class T>
1282 void AMRMultiGrid<T>::setBottomSolverEpsCushion(Real a_bottomSolverEpsCushion)
1283 {
1284  CH_assert((a_bottomSolverEpsCushion > 0.) &&
1285  (a_bottomSolverEpsCushion <= 1.));
1286  m_bottomSolverEpsCushion = a_bottomSolverEpsCushion;
1287 }
1288 
1289 template <class T>
1290 void AMRMultiGrid<T>::define(const ProblemDomain& a_coarseDomain,
1291  AMRLevelOpFactory<T>& a_factory,
1292  LinearSolver<T>* a_bottomSolver,
1293  int a_maxAMRLevels)
1294 {
1295  CH_TIME("AMRMultiGrid::define");
1296  this->clear();
1297  m_op.resize( a_maxAMRLevels, NULL);
1298  m_mg.resize( a_maxAMRLevels, NULL);
1299  m_hasInitBeenCalled.resize(a_maxAMRLevels, 'f');
1300 
1301  m_correction.resize( a_maxAMRLevels, NULL);
1302  m_residual. resize( a_maxAMRLevels, NULL);
1303  m_resC. resize( a_maxAMRLevels, NULL);
1304  m_resCopier. resize( a_maxAMRLevels);
1305  m_reverseCopier.resize( a_maxAMRLevels);
1306  m_bottomSolver = a_bottomSolver;
1307 
1308  ProblemDomain current = a_coarseDomain;
1309  for (int i = 0; i < a_maxAMRLevels; i++)
1310  {
1311  m_correction[i] = new T();
1312  m_residual[i] = new T();
1313  m_resC[i] = new T();
1314  m_mg[i]= new MultiGrid<T>();
1315  m_op[i] = a_factory.AMRnewOp(current);
1316  m_mg[i]->define(a_factory, &m_nosolve, current, m_maxDepth, m_op[i]);
1317 
1318  // Only do this if it will be used (avoiding a reference to invalid
1319  // and/or unavailable refinement ratios)
1320  if (i < a_maxAMRLevels-1)
1321  {
1322  current.refine(a_factory.refToFiner(current));
1323  }
1324  }
1325 }
1326 
1327 template<class T>
1329 {
1330  pout()<<"AMR Levels "<<m_op.size()<<"\n";
1331  for(unsigned int i=0; i<m_mg.size(); ++i)
1332  {
1333  pout()<<"MG Level "<<i<< " depth "<<m_mg[i]->m_depth
1334  <<" cycles "<<m_mg[i]->m_cycle<< "domain "<<m_mg[i]->m_topLevelDomain<<"\n";
1335  }
1336 }
1337 
1338 template<class T>
1340  Vector<T*>& a_uberResidual,
1341  int ilev, int l_max, int l_base)
1342 {
1343  if (ilev == l_max)
1344  {
1345  for (int level = l_base; level <= l_max; level++)
1346  {
1347  m_op[level]->assignLocal(*m_residual[level], *a_uberResidual[level]);
1348  m_op[level]->setToZero(*m_correction[level]);
1349  }
1350  }
1351 
1352  if (l_max == l_base)
1353  {
1354  CH_assert(ilev == l_base);
1355  m_mg[l_base]->oneCycle(*(a_uberCorrection[ilev]), *(a_uberResidual[ilev]));
1356  }
1357  else if (ilev == l_base)
1358  {
1359  m_mg[l_base]->oneCycle(*(m_correction[ilev]), *(m_residual[ilev]));
1360  m_op[ilev]->incr(*(a_uberCorrection[ilev]), *(m_correction[ilev]), 1.0);
1361  }
1362  else
1363  {
1364  //============= Downsweep ========================
1365 
1366  this->relax(*(m_correction[ilev]), *(m_residual[ilev]), ilev, m_pre);
1367  m_op[ilev]->incr(*(a_uberCorrection[ilev]), *(m_correction[ilev]), 1.0);
1368 
1369  // Set next coarser level correction to zero
1370  m_op[ilev-1]->setToZero(*(m_correction[ilev-1]));
1371 
1372  // Recompute residual on next coarser level
1373  // for the valid region NOT covered by this level.
1374  computeAMRResidualLevel(m_residual,
1375  a_uberCorrection,
1376  a_uberResidual,
1377  l_max, l_base, ilev-1,
1378  true);
1379 
1380  // Compute the restriction of the residual to the coarser level resC.
1381  m_op[ilev]->AMRRestrictS(*(m_resC[ilev]),
1382  *(m_residual[ilev]),
1383  *(m_correction[ilev]),
1384  *(m_correction[ilev-1]),
1385  *(a_uberCorrection[ilev]));
1386 
1387  // Overwrite residual on the valid region of the next coarser level
1388  // with coarsened residual from this level
1389  m_op[ilev-1]->assignCopier(*m_residual[ilev-1], *(m_resC[ilev]), m_resCopier[ilev]);
1390 
1391  //============finish Compute residual for the next coarser level======
1392 
1393  for (int img = 0; img < m_numMG; img++)
1394  {
1395  AMRVCycle(a_uberCorrection, a_uberResidual, ilev-1, l_max, l_base);
1396  }
1397 
1398  //================= Upsweep ======================
1399  //increment the correction with coarser version
1400  //m_op[ilev]->AMRProlong(*(m_correction[ilev]), *(m_correction[ilev-1]));
1401  m_op[ilev]->AMRProlongS(*(m_correction[ilev]), *(m_correction[ilev-1]),
1402  *m_resC[ilev], m_reverseCopier[ilev]);
1403  //recompute residual
1404  m_op[ilev]->AMRUpdateResidual(*(m_residual[ilev]), *(m_correction[ilev]), *(m_correction[ilev-1]));
1405 
1406  //compute correction to the correction
1407  T& dCorr = *(a_uberCorrection[ilev]); // user uberCorrection as holder for correction to correction
1408  m_op[ilev]->setToZero(dCorr);
1409  this->relax(dCorr, *(m_residual[ilev]), ilev, m_post);
1410 
1411  //correct the correction with the correction to the correction
1412  m_op[ilev]->incr(*(m_correction[ilev]), dCorr, 1.0);
1413 
1414  m_op[ilev]->assignLocal(*(a_uberCorrection[ilev]), *(m_correction[ilev]));
1415  }
1416 }
1417 
1418 #include "NamespaceFooter.H"
1419 
1420 #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:1209
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:1290
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:489
AMRLevelOp< T > & levelOp(int level)
Definition: AMRMultiGrid.H:673
A strange but true thing to make copying from one boxlayoutdata to another fast.
Definition: Copier.H:145
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:1249
#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:892
void computeAMROperator(Vector< T * > &a_lph, Vector< T * > &a_phi, int l_max, int l_base, bool a_homogeneousBC=false)
Definition: AMRMultiGrid.H:776
Vector< T * > m_residual
Definition: AMRMultiGrid.H:482
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:485
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:1185
virtual void outputAMR(Vector< T * > &a_rhs, string &a_name)
Definition: AMRMultiGrid.H:81
Vector< MultiGrid< T > * > m_mg
Definition: AMRMultiGrid.H:480
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:513
void getInfo() const
Dump out the state of the solver. AMR levels, MG levels, bottom solver, box counts, etc.
Definition: AMRMultiGrid.H:1328
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:932
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:904
const char * name(const FArrayBox &a_dummySpecializationArg)
Definition: CH_HDF5.H:864
Real m_bottomSolverEpsCushion
Definition: AMRMultiGrid.H:394
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:597
Real m_hang
Definition: AMRMultiGrid.H:378
int m_iterMin
Definition: AMRMultiGrid.H:380
Vector< RefCountedPtr< AMRMultiGridInspector< T > > > m_inspectors
Definition: AMRMultiGrid.H:498
virtual unsigned int orderOfAccuracy(void) const
Definition: AMRMultiGrid.H:217
Vector< T * > m_correction
Definition: AMRMultiGrid.H:481
Vector< T * > m_resC
Definition: AMRMultiGrid.H:483
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:680
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:1124
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:834
void relax(T &phi, T &R, int depth, int nRelax=2)
Definition: AMRMultiGrid.H:631
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:572
#define CH_STOP(tpointer)
Definition: CH_Timer.H:150
Vector< char > m_hasInitBeenCalled
Definition: AMRMultiGrid.H:491
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:525
AMRMultiGrid & operator=(const AMRMultiGrid< T > &)
size_t size() const
Definition: Vector.H:192
Vector< Vector< MGLevelOp< T > * > > getOperatorsMG()
Definition: AMRMultiGrid.H:585
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:484
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:1339
virtual void dumpLevel(T &a_data, string name)
Definition: AMRMultiGrid.H:48
Definition: AMRMultiGrid.H:233
NoOpSolver< T > m_nosolve
Definition: AMRMultiGrid.H:487
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:553
Real m_normThresh
Definition: AMRMultiGrid.H:378
virtual ~AMRLevelOpFactory()
Definition: AMRMultiGrid.H:236
Vector< AMRLevelOp< T > * > & getAMROperators()
Definition: AMRMultiGrid.H:431
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:621
void setBottomSolverEpsCushion(Real a_bottomSolverEpsCushion)
Definition: AMRMultiGrid.H:1282
Vector< AMRLevelOp< T > * > m_op
Definition: AMRMultiGrid.H:479
virtual ~AMRMultiGrid()
Definition: AMRMultiGrid.H:666
void setBottomSolver(int l_max, int l_base)
set up bottom solver for internal MG solver
Definition: AMRMultiGrid.H:1263