Chombo + EB + MF  3.2
MultiGrid.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 _MULTIGRID_H_
16 #define _MULTIGRID_H_
17 
18 #include "CH_Timer.H"
19 #include "LinearSolver.H"
20 #include "ProblemDomain.H"
21 #include "REAL.H"
22 #include "Box.H"
23 #include "parstream.H"
24 #include "RefCountedPtr.H"
25 #include <vector>
26 #include "NamespaceHeader.H"
27 
28 // Forward declaration of MGLevelOp.
29 template <typename T>
30 class MGLevelOp;
31 
32 //! \class MGLevelOpObserver
33 //! This observer class allows objects to be notified when coefficient
34 //! data for an operator changes.
35 template <typename T>
37 {
38  public:
39 
40  //! Base level Constructor. Called by all subclasses.
42  :m_op(NULL)
43  {
44  }
45 
46  //! Destructor.
47  virtual ~MGLevelOpObserver();
48 
49  //! Use this to implement the response of the observer to changes
50  //! in the observee.
51  //! \param a_operator The operator whose state has changed. Passed to
52  //! the observer by the operator itself.
53  virtual void operatorChanged(const MGLevelOp<T>& a_operator)
54  {
55  }
56 
57  //! This gets called by the observee when this observer is added to it.
58  //! DO NOT CALL THIS METHOD DIRECTLY.
59  //! \param a_observee The object being observed by this observer.
60  void setObservee(MGLevelOp<T>* a_observee)
61  {
62  CH_assert(m_op == NULL);
63  m_op = a_observee;
64  }
65 
66  //! This gets called by the observee when it is destroyed.
67  //! DO NOT CALL THIS METHOD DIRECTLY.
69  {
70  m_op = NULL;
71  }
72 
73  private:
74 
75  // A pointer to the observed operator.
77 
78  // Forbidden ops.
81 };
82 
83 //! \class MGLevelOp
84 //! This class handles the additional tasks of coordinating operations
85 //! between this level and the level one coarser than this 'level'.
86 //! This class is also an MGLevelOpObserver, so objects of this type can
87 //! observe other MGLevelOp objects.
88 template <typename T>
89 class MGLevelOp : public LinearOp<T>, public MGLevelOpObserver<T>
90 {
91  public:
92 
93  //! Constructor.
95  :LinearOp<T>(),
96  MGLevelOpObserver<T>(),
97  m_observers(),
99  {
100  }
101 
102  //! Destructor.
103  virtual ~MGLevelOp()
104  {
105  // Go through the list of observers and clear their observees.
106  for (size_t i = 0; i < m_observers.size(); ++i)
108  }
109 
110  ///
111  /**
112  Create a coarsened (by two) version of the input data. This does not include averaging
113  the data. So if a_fine is over a Box of (0, 0, 0) (63, 63, 63), a_fine should
114  be over a Box (0, 0, 0) (31, 31, 31).
115  */
116  virtual void createCoarser(T& a_coarse, const T& a_fine, bool ghosted) = 0;
117 
118  ///
119  /**
120  Use your relaxation operator to remove the high frequency wave numbers from
121  the correction so that it may be averaged to a coarser refinement.
122  A point relaxtion scheme, for example takes the form
123  a_correction -= lambda*(L(a_correction) - a_residual).
124  */
125  virtual void relax(T& a_correction, const T& a_residual, int a_iterations) = 0 ;
126 
127  /// specialized no-fine relax function, useful for full-multigrid schemes, defaults to regular relax
128  /**
129  Necessary to implement FAS correctly
130  */
131  virtual void relaxNF(T& a_phi, const T* a_phiCoarse, const T& a_rhs, int a_iterations)
132  {
133  this->relax(a_phi, a_rhs, a_iterations);
134  }
135 
136  ///
137  /**
138  calculate restricted residual
139  a_resCoarse[2h] = I[h->2h] (rhsFine[h] - L[h](phiFine[h])
140  */
141  virtual void restrictResidual(T& a_resCoarse, T& a_phiFine, const T& a_rhsFine) = 0;
142 
143  /// full-multigrid version of restrictResidual, useful for FAS-type schemes. Defaults to standard restriction
144  virtual void restrictResidual(T& a_resCoarse, T& a_phiFine, const T* a_phiCoarse, const T& a_rhsFine, bool homogeneous)
145  {
146  restrictResidual(a_resCoarse, a_phiFine, a_rhsFine);
147  }
148 
149  /// simple restriction function
150  virtual void restrictR(T& a_phiCoarse, const T& a_phiFine)
151  {
152 
153  }
154 
155  ///
156  /**
157  correct the fine solution based on coarse correction
158  a_phiThisLevel += I[2h->h](a_correctCoarse)
159  */
160  virtual void prolongIncrement(T& a_phiThisLevel, const T& a_correctCoarse) = 0;
161 
162  //! This adds a new observer to this operator. Note that this operator does not
163  //! own the resources for the observer, so you must be careful to ensure that
164  //! the observer does not go out of scope while the operator lives. If the observer
165  //! is already observing the operator, this operation has no effect.
166  //! \param a_observer An observer to be notified when the operator changes.
168  {
169  if (std::find(m_observers.begin(), m_observers.end(), a_observer) == m_observers.end())
170  {
171  m_observers.push_back(a_observer);
172  m_coarseningFactors.push_back(-1);
173  a_observer->setObservee(this);
174  }
175  }
176 
177  /// Apply an operator
178  /**
179  useful for full-multigrid/FAS types of algorithms
180  */
181  virtual void applyOpMg(T& a_lhs, T& a_phi, T* a_phiCoarse, bool a_homogeneous)
182  {
183 
184  }
185 
186  // For use by full-multigrid schemes if we have a coarser level. Just do standard residual by default.
187  virtual void residualNF(T& a_lhs, T& a_phi, const T* a_phiCoarse, const T& a_rhs, bool a_homogeneous = false)
188  {
189  this->residual(a_lhs, a_phi, a_rhs, a_homogeneous);
190  }
191 
192  //! Removes the given observer from the operator's list of observers. This
193  //! method has no effect if the observer is not found in the list.
194  //! \param a_observer An observer to be removed from the operator's observer list.
196  {
197  typename std::vector<MGLevelOpObserver<T>*>::iterator iter =
198  std::find(m_observers.begin(), m_observers.end(), a_observer);
199  if (iter != m_observers.end())
200  {
201  (*iter)->clearObservee();
202  m_observers.erase(iter);
203  // Remove the coarsening factor entry, too.
204  size_t index = std::distance(m_observers.begin(), iter);
205  m_coarseningFactors.erase(m_coarseningFactors.begin() + index);
206  }
207  }
208 
209  //! Adds a MGLevelOp observer to this operator that is notified when the
210  //! operator changes so that it can coarsen its corresponding data. This
211  //! is used to implement dependencies between multigrid operators at different
212  //! grid levels.
213  //! \param a_operator The multigrid operator that will observe this operator.
214  //! No operator may observe itself. If \a a_operator is
215  //! already observing this operator, this method has no
216  //! effect.
217  //! \param a_coarseningFactor The refinement factor between the grid
218  //! level for \a a_operator and that for this one.
220  int a_coarseningFactor)
221  {
222  CH_assert(a_operator != NULL);
223  CH_assert(a_operator != this);
224  // Add the operator to the list of observers.
225  if (std::find(m_observers.begin(), m_observers.end(), a_operator) == m_observers.end())
226  {
227  // Add the operator to our list of observers.
228  addObserver(a_operator);
229 
230  // Jot down this operator's multigrid coarsening factor.
231  m_coarseningFactors.back() = a_coarseningFactor;
232  }
233  }
234 
235  //! This should be called whenever the operator's data is updated.
237  {
238  //! Notify our observers that the operator's data has changed,
239  //! taking special care to pass multigrid operators their
240  //! level.
241  for (int i = 0; i < m_observers.size(); ++i)
242  {
243  if (m_coarseningFactors[i] != -1)
244  {
245  MGLevelOp<T>* op = dynamic_cast<MGLevelOp<T>*>(m_observers[i]);
246  CH_assert(op != NULL);
247  // This observer is a multigrid operator. Pass it its grid level.
249  }
250  else
251  {
252  // Simply tell it that its operator has changed.
253  m_observers[i]->operatorChanged(*this);
254  }
255  }
256  }
257 
258  //! Overload this method to implement the response of a multigrid
259  //! operator to a change in a finer operator that it is observing. This
260  //! allows the coarser observer to coarsen any data that has been changed
261  //! in the finer operator.
262  //! \param a_operator The operator whose state has changed.
263  //! \param a_coarseningFactor The factor by which this operator's data is
264  //! coarsened with respect to that of \a a_operator.
265  virtual void finerOperatorChanged(const MGLevelOp<T>& a_operator,
266  int a_coarseningFactor)
267  {
268  }
269 
270  //! Returns the number of objects observing this operator.
271  int numObservers() const
272  {
273  return m_observers.size();
274  }
275 
276  private:
277 
278  //! List of observers for the operator.
279  std::vector<MGLevelOpObserver<T>*> m_observers;
280 
281  //! Multigrid coarsening factors.
282  std::vector<int> m_coarseningFactors;
283 
284  // Forbidden operations.
285  MGLevelOp(const MGLevelOp&);
286  MGLevelOp& operator=(const MGLevelOp&);
287 };
288 
289 ///
290 /**
291  Factory class for generating MGLevelOps
292  */
293 template <class T>
295 {
296  public:
297 
298  //! Base class constructor.
300  {
301  }
302 
303  //! Destructor.
305  {
306  }
307 
308  /**
309  Create an operator at an index space = coarsen(a_fineIndexSpace, 2^a_depth)
310  Return NULL if no such Multigrid level can be created at this a_depth.
311  If a_homoOnly = true, then only homogeneous boundary conditions will be needed.
312  */
313  virtual MGLevelOp<T>* MGnewOp(const ProblemDomain& a_FineindexSpace,
314  int a_depth,
315  bool a_homoOnly = true) = 0;
316 
317  private:
318 
319  // Forbidden operations.
322 };
323 
324 ///
325 /**
326  Class to execute geometric multigrid. This class is not meant to be particularly
327  user-friendly and a good option for people who want something a tad less raw is to
328  use AMRMultiGrid instead.
329  */
330 template <class T>
332 {
333 public:
334  MultiGrid();
335  virtual ~MultiGrid();
336 
337  /// Function to define a MultiGrid object.
338  /**
339  a_factory is the factory for generating operators.
340  a_bottomSolver is called at the bottom of v-cycle.
341  a_domain is the problem domain at the top of the vcycle.
342  maxDepth defines the location of the bottom of the v-cycle.
343  The vycle will terminate (hit bottom) when the factory returns NULL for a paticular
344  depth if maxdepth = -1. Otherwise the vycle terminates at maxdepth.
345  */
346  virtual void define(MGLevelOpFactory<T>& a_factory,
347  LinearSolver<T>* a_bottomSolver,
348  const ProblemDomain& a_domain,
349  int a_maxDepth = -1,
350  MGLevelOp<T>* a_finestLevelOp = NULL);
351 
352  ///
353  /**
354  solve L(a_phi) = a_rhs. Tolerance is how much you want the norm of the error reduced.
355  verbosity is how chatty you want the function to be. maxIterations is the maximum number
356  of v-cycles. This does the whole residual correction switcharoo and calls oneCycle up to
357  maxIterations times, evaluating the residual as it goes.
358  */
359  virtual void solve(T& a_phi, const T& a_rhs, Real a_tolerance, int a_maxIterations, int verbosity= 0);
360 
361  /// Execute ONE v-cycle of multigrid.
362  /**
363  If you want the solution to converge, you need to iterate this.
364  See solve() or AMRMultiGrid::solve for a more automatic solve() function.
365  This operates residual-correction form of solution
366  so all boundary conditions are assumed to be homogeneous.
367  L(a_e) = a_res
368  */
369  virtual void oneCycle(T& a_e, const T& a_res);
370 
371  // used by AMRMultiGrid--not a part of the public API
372  void init(const T& a_correction, const T& a_residual);
373 
374  //internal function
375  void cycle(int a_depth, T& a_correction, const T& a_residual);
376 
377  //internal function
378  void clear();
379 
380  // used by AMRMultiGrid--not a part of the public API
381  void setBottomSolver(LinearSolver<T>* a_bottomSolver);
382 
383  ///
384  /**
385  Public solver parameters. m_pre and m_post are the ones
386  that usually get set and are the number of relaxations performed
387  before and after multigrid recursion. See AMRMultiGrid for a
388  more user-friendly interface.
389  */
393 
394  ///
395  /**
396  for changing coefficients --- not for the faint of heart.
397  */
399 
401 protected:
402  bool m_defined;
403 
405 
409  std::vector<bool> m_ownOp;
410 
411 private:
412  MultiGrid(const MultiGrid<T>& a_opin)
413  {
414  MayDay::Error("invalid operator");
415  }
416 
417  void operator=(const MultiGrid<T>& a_opin)
418  {
419  MayDay::Error("invalid operator");
420  }
421 };
422 
423 // *******************************************************
424 // MultiGrid Implementation
425 // *******************************************************
426 
427 //-----------------------------------------------------------------------
428 template <class T>
431 {
432  Vector<MGLevelOp<T>*> retval = m_op;
433  return retval;
434 }
435 //-----------------------------------------------------------------------
436 
437 //-----------------------------------------------------------------------
438 template <class T>
440  :m_pre(3),
441  m_post(3),
442  m_bottom(0),
443  m_cycle(1),
444  m_numMG(1),
445  m_homogeneous(true),
446  m_defined(false)
447 {
448 }
449 //-----------------------------------------------------------------------
450 
451 //-----------------------------------------------------------------------
452 template <class T>
454 {
455  clear();
456 }
457 //-----------------------------------------------------------------------
458 
459 //-----------------------------------------------------------------------
460 template <class T>
462  LinearSolver<T>* a_bottomSolver,
463  const ProblemDomain& a_domain,
464  int a_maxDepth,
465  MGLevelOp<T>* a_finestOp)
466 {
467 
468  CH_TIME("MultiGrid::define");
469  clear();
470  m_depth = 0;
471  m_bottomSolver = a_bottomSolver;
472  m_residual.resize(m_depth); // the zero'th entry is left null constructed
473  m_correction.resize(m_depth); // the zero'th entry is left null constructed
474  m_op.resize(m_depth);
475  m_ownOp.resize(0);
476  m_topLevelDomain = a_domain;
477 
478  m_bottomCells = a_domain.domainBox().numPts();
479 
480  MGLevelOp<T>* nextOp;
481  bool ownOp;
482  if (a_finestOp == NULL)
483  {
484  nextOp = a_factory.MGnewOp(a_domain, 0, true);
485  ownOp = true;
486  }
487  else
488  {
489  nextOp = a_finestOp;
490  ownOp = false;
491  }
492  while (nextOp != NULL)
493  {
494  m_residual.push_back(new T());
495  m_correction.push_back(new T());
496  m_op.push_back(nextOp);
497  m_ownOp.push_back(ownOp);
498  m_depth++;
499  if ((m_depth < a_maxDepth) || (a_maxDepth < 0))
500  {
501  // Create a coarser multigrid operator.
502  nextOp = a_factory.MGnewOp(a_domain , m_depth, true);
503  ownOp = true;
504  if (nextOp != NULL)
505  {
506  // This multigrid operator should observe the next finer operator. -JNJ
507  // NOTE: the coarsening ratio is hardwired to 2 here.
508  m_op.back()->addCoarserObserver(nextOp, 2);
509  for (int i=0; i<CH_SPACEDIM; ++i) m_bottomCells /= 2;
510  }
511  }
512  else
513  {
514  nextOp = NULL;
515  }
516  }
517  m_defaultDepth = m_depth;
518  m_bottomSolver->define(m_op[m_depth-1], true);
519  m_defined = true;
520 }
521 //-----------------------------------------------------------------------
522 
523 //-----------------------------------------------------------------------
524 template <class T>
526 {
527  m_bottomSolver = a_bottomSolver;
528  m_bottomSolver->define(m_op[m_depth-1], true);
529 }
530 //-----------------------------------------------------------------------
531 
532 //-----------------------------------------------------------------------
533 template <class T>
534 void MultiGrid<T>::init(const T& a_e, const T& a_residual)
535 {
536  CH_TIME("MutliGrid::init");
537  if (m_depth > 1)
538  {
539  m_op[0]->createCoarser(*(m_residual[1]), a_residual, false);
540  m_op[0]->createCoarser(*(m_correction[1]), a_e, true);
541  }
542 
543  for (int i = 2; i < m_depth; i++)
544  {
545  m_op[i-1]->createCoarser(*(m_residual[i]), *(m_residual[i-1]), false);
546  m_op[i-1]->createCoarser(*(m_correction[i]), *(m_correction[i-1]), true);
547  }
548 }
549 //-----------------------------------------------------------------------
550 
551 //-----------------------------------------------------------------------
552 template <class T>
553 void MultiGrid<T>::solve(T& a_phi, const T& a_rhs, Real a_tolerance, int a_maxIterations, int verbosity)
554 {
555  CH_TIME("MultiGrid::solve");
556  this->init(a_phi, a_rhs);
557 
558  T correction, residual;
559  m_op[0]->create(correction, a_phi);
560  m_op[0]->create(residual, a_rhs);
561  m_op[0]->setToZero(a_phi);
562  m_op[0]->residual(residual, a_phi, a_rhs, false);
563 
564  Real errorno = m_op[0]->norm(residual, 0);
565  if (verbosity > 2)
566  {
567  pout() << "multigrid::solve initial residual = " << errorno << std::endl;
568  }
569  Real compval = a_tolerance*errorno;
570  Real epsilon = 1.0e-16;
571  compval = Max(compval, epsilon);
572  Real error = errorno;
573  int iter = 0;
574  while ((error > compval) && (error > a_tolerance*errorno) && (iter < a_maxIterations))
575  {
576  m_op[0]->setToZero(correction);
577  m_op[0]->residual(residual, a_phi, a_rhs, false);
578  error = m_op[0]->norm(residual, 0);
579  if (verbosity > 3)
580  {
581  pout() << "multigrid::solve iter = " << iter << ", residual = " << error << std::endl;
582  }
583 
584  this->cycle(0, correction, residual);
585  m_op[0]->incr(a_phi, correction, 1.0);
586 
587  iter++;
588  }
589  if (verbosity > 2)
590  {
591  pout() << "multigrid::solve final residual = " << error << std::endl;
592  }
593 
594  m_op[0]->clear(correction);
595  m_op[0]->clear(residual);
596 }
597 //-----------------------------------------------------------------------
598 
599 //-----------------------------------------------------------------------
600 template <class T>
601 void MultiGrid<T>::oneCycle(T& a_e, const T& a_residual)
602 {
603  CH_TIME("Multigrid::oneCycle");
604 
605  //this->init(a_e, a_residual); oneCycle is only used by AMRMultigrid
606  if (m_homogeneous)
607  {
608  cycle(0, a_e, a_residual);
609  }
610  else
611  {
612  T correction, residual;
613  m_op[0]->create(correction, a_e);
614  m_op[0]->create(residual, a_residual);
615 
616  m_op[0]->residual(residual, a_e, a_residual);
617 
618  m_op[0]->setToZero(correction);
619  this->cycle(0, correction, residual);
620 
621  m_op[0]->incr(a_e, correction, 1.0);
622  m_op[0]->clear(correction);
623  m_op[0]->clear(residual);
624  }
625 }
626 //-----------------------------------------------------------------------
627 
628 //-----------------------------------------------------------------------
629 template <class T>
630 void MultiGrid<T>::cycle(int depth, T& correction, const T& residual)
631 {
632  CH_TIME("Multigrid::cycle");
633 
634  // currently I can't drop a timer in this function because it is recursive. doh
635  if (depth == m_depth-1)
636  {
637  CH_TIME("Multigrid::cycle:bottom-solve");
638  if (m_bottomCells == 1)
639  {
640  m_op[depth ]->relax(correction, residual, 1);
641  }
642  else
643  {
644  m_op[depth ]->relax(correction, residual, m_bottom);
645  m_bottomSolver->solve(correction, residual);
646  }
647  //m_op[depth ]->relax(correction, residual, 1);
648  }
649  else
650  {
651  int cycles = m_cycle;
652  if ( cycles < 0 )
653  {
654  cycles = -cycles;
655  // F cycle multigrid
656  m_op[depth ]->restrictResidual(*(m_residual[depth+1]), correction, residual);
657  m_op[depth ]->setToZero(*(m_correction[depth+1]));
658 
659  // recursive call
660  cycle(depth+1, *(m_correction[depth+1]), *(m_residual[depth+1]));
661 
662  m_op[depth ]->prolongIncrement(correction, *(m_correction[depth+1]));
663 
664  m_op[depth ]->relax(correction, residual, m_pre);
665 
666  for (int img = 0; img < cycles; img++)
667  {
668  m_op[depth ]->restrictResidual(*(m_residual[depth+1]), correction, residual);
669  m_op[depth ]->setToZero(*(m_correction[depth+1]));
670 
671  m_cycle = 1; // hack to get a V-cycle
672  cycle(depth+1, *(m_correction[depth+1]), *(m_residual[depth+1]));
673  m_cycle = -cycles;
674  //
675  m_op[depth ]->prolongIncrement(correction, *(m_correction[depth+1]));
676  }
677 
678  m_op[depth ]->relax(correction, residual, m_post);
679  }
680  else
681  {
682  m_op[depth ]->relax( correction, residual, m_pre );
683 
684  m_op[depth ]->restrictResidual(*(m_residual[depth+1]), correction, residual);
685  m_op[depth ]->setToZero(*(m_correction[depth+1]));
686 
687  for (int img = 0; img < cycles; img++)
688  {
689  cycle(depth+1, *(m_correction[depth+1]), *(m_residual[depth+1]));
690  }
691  m_op[depth ]->prolongIncrement(correction, *(m_correction[depth+1]));
692 
693  m_op[depth ]->relax(correction, residual, m_post);
694  }
695  }
696 }
697 //-----------------------------------------------------------------------
698 
699 //-----------------------------------------------------------------------
700 template <class T>
702 {
703  std::vector<bool>::reverse_iterator it = m_ownOp.rbegin();
704  for (int i=m_op.size()-1; i>=0; --i, ++it)
705  {
706  if (*it) delete m_op[i];
707  delete m_correction[i];
708  delete m_residual[i];
709  }
710  m_op.resize(0);
711  m_residual.resize(0);
712  m_correction.resize(0);
713  m_defined = false;
714 }
715 //-----------------------------------------------------------------------
716 
717 //-----------------------------------------------------------------------
718 // This destructor implementation had to appear below the MGLevelOp class
719 // definition.
720 template <typename T>
722 {
723  if (m_op != NULL)
724  {
725  // Remove this observer from its observee.
726  // For now, we skip this, since we haven't established rules
727  // for scoping between multigrid objects and AMRMultigrid objects
728  // (which own the AMR operators).
729  m_op->removeObserver(this);
730  }
731 }
732 //-----------------------------------------------------------------------
733 
734 #include "NamespaceFooter.H"
735 #endif /*_MULTIGRID_H_*/
std::ostream & pout()
Use this in place of std::cout for program output.
void addObserver(MGLevelOpObserver< T > *a_observer)
Definition: MultiGrid.H:167
MGLevelOp()
Constructor.
Definition: MultiGrid.H:94
Definition: LinearSolver.H:28
void clear()
Definition: MultiGrid.H:701
int m_cycle
Definition: MultiGrid.H:390
#define CH_SPACEDIM
Definition: SPACE.H:51
#define CH_assert(cond)
Definition: CHArray.H:37
A class to facilitate interaction with physical boundary conditions.
Definition: ProblemDomain.H:141
Vector< T * > m_residual
Definition: MultiGrid.H:407
int m_defaultDepth
Definition: MultiGrid.H:390
MGLevelOpObserver & operator=(const MGLevelOpObserver &)
MGLevelOpObserver()
Base level Constructor. Called by all subclasses.
Definition: MultiGrid.H:41
one dimensional dynamic array
Definition: Vector.H:53
void removeObserver(MGLevelOpObserver< T > *a_observer)
Definition: MultiGrid.H:195
void setObservee(MGLevelOp< T > *a_observee)
Definition: MultiGrid.H:60
virtual ~MGLevelOpFactory()
Destructor.
Definition: MultiGrid.H:304
void init(const T &a_correction, const T &a_residual)
Definition: MultiGrid.H:534
virtual void define(MGLevelOpFactory< T > &a_factory, LinearSolver< T > *a_bottomSolver, const ProblemDomain &a_domain, int a_maxDepth=-1, MGLevelOp< T > *a_finestLevelOp=NULL)
Function to define a MultiGrid object.
Definition: MultiGrid.H:461
Vector< MGLevelOp< T > * > m_op
Definition: MultiGrid.H:406
virtual void restrictResidual(T &a_resCoarse, T &a_phiFine, const T *a_phiCoarse, const T &a_rhsFine, bool homogeneous)
full-multigrid version of restrictResidual, useful for FAS-type schemes. Defaults to standard restric...
Definition: MultiGrid.H:144
virtual void oneCycle(T &a_e, const T &a_res)
Execute ONE v-cycle of multigrid.
Definition: MultiGrid.H:601
bool m_homogeneous
Definition: MultiGrid.H:391
int m_bottomCells
Definition: MultiGrid.H:404
Definition: MultiGrid.H:36
std::vector< MGLevelOpObserver< T > * > m_observers
List of observers for the operator.
Definition: MultiGrid.H:279
MultiGrid()
Definition: MultiGrid.H:439
void cycle(int a_depth, T &a_correction, const T &a_residual)
Definition: MultiGrid.H:630
int m_post
Definition: MultiGrid.H:390
std::vector< bool > m_ownOp
Definition: MultiGrid.H:409
void setBottomSolver(LinearSolver< T > *a_bottomSolver)
Definition: MultiGrid.H:525
ProblemDomain m_topLevelDomain
Definition: MultiGrid.H:400
MGLevelOp & operator=(const MGLevelOp &)
int m_depth
Definition: MultiGrid.H:390
virtual void restrictR(T &a_phiCoarse, const T &a_phiFine)
simple restriction function
Definition: MultiGrid.H:150
int numObservers() const
Returns the number of objects observing this operator.
Definition: MultiGrid.H:271
#define CH_TIME(name)
Definition: CH_Timer.H:82
Definition: MultiGrid.H:294
bool m_defined
Definition: MultiGrid.H:402
virtual void residualNF(T &a_lhs, T &a_phi, const T *a_phiCoarse, const T &a_rhs, bool a_homogeneous=false)
Definition: MultiGrid.H:187
double Real
Definition: REAL.H:33
Definition: MultiGrid.H:30
virtual void prolongIncrement(T &a_phiThisLevel, const T &a_correctCoarse)=0
int m_numMG
Definition: MultiGrid.H:390
const Box & domainBox() const
Returns the logical computational domain.
Definition: ProblemDomain.H:887
int m_bottom
Definition: MultiGrid.H:390
MultiGrid(const MultiGrid< T > &a_opin)
Definition: MultiGrid.H:412
virtual void finerOperatorChanged(const MGLevelOp< T > &a_operator, int a_coarseningFactor)
Definition: MultiGrid.H:265
void clearObservee()
Definition: MultiGrid.H:68
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.
MGLevelOpFactory & operator=(const MGLevelOpFactory &)
void operator=(const MultiGrid< T > &a_opin)
Definition: MultiGrid.H:417
Definition: MultiGrid.H:331
virtual void define(LinearOp< T > *a_operator, bool a_homogeneous=false)=0
Vector< T * > m_correction
Definition: MultiGrid.H:408
size_t numPts() const
Definition: LinearSolver.H:156
virtual MGLevelOp< T > * MGnewOp(const ProblemDomain &a_FineindexSpace, int a_depth, bool a_homoOnly=true)=0
virtual void relax(T &a_correction, const T &a_residual, int a_iterations)=0
void addCoarserObserver(MGLevelOp< T > *a_operator, int a_coarseningFactor)
Definition: MultiGrid.H:219
virtual void createCoarser(T &a_coarse, const T &a_fine, bool ghosted)=0
virtual void applyOpMg(T &a_lhs, T &a_phi, T *a_phiCoarse, bool a_homogeneous)
Apply an operator.
Definition: MultiGrid.H:181
virtual void restrictResidual(T &a_resCoarse, T &a_phiFine, const T &a_rhsFine)=0
Vector< MGLevelOp< T > * > getAllOperators()
Definition: MultiGrid.H:430
virtual ~MultiGrid()
Definition: MultiGrid.H:453
MGLevelOpFactory()
Base class constructor.
Definition: MultiGrid.H:299
void notifyObserversOfChange()
This should be called whenever the operator's data is updated.
Definition: MultiGrid.H:236
T Max(const T &a_a, const T &a_b)
Definition: Misc.H:39
std::vector< int > m_coarseningFactors
Multigrid coarsening factors.
Definition: MultiGrid.H:282
virtual ~MGLevelOp()
Destructor.
Definition: MultiGrid.H:103
MGLevelOp< T > * m_op
Definition: MultiGrid.H:76
virtual ~MGLevelOpObserver()
Destructor.
Definition: MultiGrid.H:721
virtual void operatorChanged(const MGLevelOp< T > &a_operator)
Definition: MultiGrid.H:53
int m_pre
Definition: MultiGrid.H:390
virtual void relaxNF(T &a_phi, const T *a_phiCoarse, const T &a_rhs, int a_iterations)
specialized no-fine relax function, useful for full-multigrid schemes, defaults to regular relax ...
Definition: MultiGrid.H:131
virtual void solve(T &a_phi, const T &a_rhs, Real a_tolerance, int a_maxIterations, int verbosity=0)
Definition: MultiGrid.H:553
virtual void residual(T &a_lhs, const T &a_phi, const T &a_rhs, bool a_homogeneous=false)=0
LinearSolver< T > * m_bottomSolver
Definition: MultiGrid.H:392