Chombo + EB  3.2
FASMultiGrid.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 
12 #ifndef _FASMULTIGRID_H_
13 #define _FASMULTIGRID_H_
14 
15 #include "AMRMultiGrid.H" // So we can cast MGLevelOp's as AMRLevelOp's
16 #include "MultiGrid.H"
17 
18 
19 #include "NamespaceHeader.H"
20 
21 template <class T>
22 class FASMultiGrid : public MultiGrid<T>
23 {
24 public:
25  FASMultiGrid();
26  virtual ~FASMultiGrid();
27 
28 
29  /// Function to define a FASMultiGrid object.
30  /* Author: Jamie Parkinson, February 2017.
31  An extension of AMRFASMultiGrid to even coarser levels*/
32  /**
33  a_factory is the factory for generating operators.
34  a_bottomSolver is called at the bottom of v-cycle.
35  a_domain is the problem domain at the top of the vcycle.
36  maxDepth defines the location of the bottom of the v-cycle.
37  The vycle will terminate (hit bottom) when the factory returns NULL for a paticular
38  depth if maxdepth = -1. Otherwise the vycle terminates at maxdepth.
39  */
40  virtual void define(MGLevelOpFactory<T>& a_factory,
41  LinearSolver<T>* a_bottomSolver,
42  const ProblemDomain& a_domain,
43  int a_maxDepth = -1,
44  MGLevelOp<T>* a_finestLevelOp = NULL);
45 
46  ///
47  /**
48  solve L(a_phi) = a_rhs. Tolerance is how much you want the norm of the error reduced.
49  verbosity is how chatty you want the function to be. maxIterations is the maximum number
50  of v-cycles. This does the whole residual correction switcharoo and calls oneCycle up to
51  maxIterations times, evaluating the residual as it goes.
52  */
53  virtual void solve(T& a_phi, const T& a_rhs, Real a_tolerance, int a_maxIterations, int verbosity= 0);
54 
55  /// Execute ONE v-cycle of multigrid.
56  /**
57  If you want the solution to converge, you need to iterate this.
58  See solve() or AMRFASMultiGrid::solve for a more automatic solve() function.
59  */
60  virtual void oneCycle(T& a_e, const T& a_res);
61  void oneCycle(T& a_e, const T& a_res, T* a_phiCoarse);
62 
63  //internal function
64  virtual void cycle(int a_depth, T& a_correction, const T& a_residual);
65 
66 protected:
67 
69 
70 
71 
72 };
73 
74 
75 //Constructor
76 template <class T>
78 : MultiGrid<T>()
79  {
80  }
81 
82 
83 // Destructor
84 template <class T>
86 {
87  CH_TIME("~FASMultiGrid");
88 }
89 
90 template <class T>
92  LinearSolver<T>* a_bottomSolver,
93  const ProblemDomain& a_domain,
94  int a_maxDepth,
95  MGLevelOp<T>* a_finestOp)
96  {
97  MultiGrid<T>::define(a_factory,a_bottomSolver,a_domain, a_maxDepth, a_finestOp);
98 
99  // Should never be homogeneous in the FAS case, as we're solving the full problem
100  // so need domain BCs
101  this->m_homogeneous = false;
102 
103 }
104 
105 
106 template <class T>
107 void FASMultiGrid<T>::solve(T& a_phi, const T& a_rhs, Real a_tolerance, int a_maxIterations, int verbosity)
108 {
109  CH_TIME("FASMultiGrid::solve");
110  this->init(a_phi, a_rhs);
111 
112  T correction, residual;
113  this->m_op[0]->create(correction, a_phi);
114  this->m_op[0]->create(residual, a_rhs);
115  this->m_op[0]->setToZero(a_phi);
116  this->m_op[0]->residual(residual, a_phi, a_rhs, false);
117 
118  Real errorno = this->m_op[0]->norm(residual, 0);
119  if (verbosity > 2)
120  {
121  pout() << "FASmultigrid::solve initial residual = " << errorno << std::endl;
122  }
123  Real compval = a_tolerance*errorno;
124  Real epsilon = 1.0e-16;
125  compval = Max(compval, epsilon);
126  Real error = errorno;
127  int iter = 0;
128  while ((error > compval) && (error > a_tolerance*errorno) && (iter < a_maxIterations))
129  {
130  this->m_op[0]->setToZero(correction);
131  this->m_op[0]->residual(residual, a_phi, a_rhs, false);
132  error = this->m_op[0]->norm(residual, 0);
133  if (verbosity > 3)
134  {
135  pout() << "FASMultigrid::solve iter = " << iter << ", residual = " << error << std::endl;
136  }
137 
138  this->cycle(0, correction, residual);
139  this->m_op[0]->incr(a_phi, correction, 1.0);
140 
141  iter++;
142  }
143  if (verbosity > 2)
144  {
145  pout() << "FASMultigrid::solve final residual = " << error << std::endl;
146  }
147 
148  this->m_op[0]->clear(correction);
149  this->m_op[0]->clear(residual);
150 }
151 
152 // Do one cycle
153 // in FAS we solve the full problem, so instead of having the
154 // error and residual, we have phi and the rhs
155 
156 template <class T>
157 void FASMultiGrid<T>::oneCycle(T& a_phi, const T& a_rhs)
158 {
159  // We should never do this - always need a coarse phi
160  CH_assert(0);
161 }
162 
163 template <class T>
164 void FASMultiGrid<T>::oneCycle(T& a_phi, const T& a_rhs, T* a_phiCoarse)
165 {
166  CH_TIME("FASMultigrid::oneCycle");
167 
168  // First need to create coarser versions of phiCoarse at all MG depths
169  m_phiCoarse.resize(this->m_depth);
170  m_phiCoarse[0] = a_phiCoarse;
171 
172  for (int depth = 1; depth < this->m_depth; depth++)
173  {
174  if (m_phiCoarse[depth-1] == NULL)
175  {
176  m_phiCoarse[depth] = NULL;
177  }
178  else
179  {
180  const T& phiCoarseGridFiner = *m_phiCoarse[depth-1];
181 
182  m_phiCoarse[depth] = new T();
183 
184  T scratchCrse;
185  this->m_op[depth]->create(scratchCrse, phiCoarseGridFiner);
186  this->m_op[depth]->setToZero(scratchCrse);
187 
188  this->m_op[depth]->createCoarser(*(m_phiCoarse[depth]), phiCoarseGridFiner, true);
189  this->m_op[depth]->restrictR(*(m_phiCoarse[depth]), phiCoarseGridFiner);
190 
191  // If the coarse level is so coarse that it only has one cell then we can't do
192  // CF interpolation with it, so delete it
193  const DisjointBoxLayout& dbl = m_phiCoarse[depth]->disjointBoxLayout();
194  const Box& coarseDomBox = dbl.physDomain().domainBox();
195  if (coarseDomBox.smallEnd() == coarseDomBox.bigEnd())
196  {
197  delete m_phiCoarse[depth];
198  m_phiCoarse[depth] = NULL;
199  }
200  }
201  }
202 
203  // We pass in a_phi and a_rhs here
204  this->cycle(0, a_phi, a_rhs);
205 
206  //Cleanup
207  // Note that m_phiCoarse[0] = a_phiCoarse, so we start the loop from i=1
208  for (int i=1; i<m_phiCoarse.size(); i++)
209  {
210  if (m_phiCoarse[i]!=NULL)
211  {
212  delete m_phiCoarse[i];
213  m_phiCoarse[i] = NULL;
214  }
215  }
216 }
217 //-----------------------------------------------------------------------
218 
219 //-----------------------------------------------------------------------
220 template <class T>
221 void FASMultiGrid<T>::cycle(int depth, T& a_phi, const T& a_rhs)
222 {
223  CH_TIME("FASMultigrid::cycle");
224 
225  // currently I can't drop a timer in this function because it is recursive. doh
226  if (depth == this->m_depth-1)
227  {
228  CH_TIME("FASMultigrid::cycle:bottom-solve");
229 
230  // Store the original phi so we can compute the 'correction' later
231  this->m_op[0]->setToZero(*this->m_correction[depth]);
232  this->m_op[0]->incr(*this->m_correction[depth], a_phi, -1.0);
233 
234  // On the coarsest level, just do loads of relaxing as the bottom solve
235  this->m_op[depth ]->relaxNF(a_phi, m_phiCoarse[depth], a_rhs, this->m_bottom);
236 
237  // Subtract off new phi to get actual correction from this level
238  this->m_op[0]->incr(*this->m_correction[depth], a_phi, 1.0);
239 
240 
241  }
242  else
243  {
244  int cycles = this->m_cycle;
245  if ( cycles < 0 )
246  {
247  MayDay::Error("FASMultigrid::cycle - non-V cycle options not implemented");
248  }
249  else
250  {
251  // Create some temporary data structures
252  T phi_coarse, op_coarse;
253 
254  this->m_op[depth+1]->create(phi_coarse, *(this->m_correction[depth+1]));
255  this->m_op[depth+1]->setToZero(phi_coarse);
256 
257 
258  this->m_op[depth+1]->create(op_coarse, *(this->m_correction[depth+1]));
259  this->m_op[depth+1]->setToZero(op_coarse);
260 
261  // Store initial guess at phi on this level
262  // Note this is negative, as the eventual correction = u_new - u_old (and this is u_old)
263  this->m_op[depth]->setToZero( *(this->m_correction[depth]));
264  this->m_op[depth]->incr(*(this->m_correction[depth]), a_phi, -1.0);
265 
266  // Relax phi towards actual phi
267  this->m_op[depth ]->relaxNF( a_phi, m_phiCoarse[depth], a_rhs, this->m_pre );
268 
269 
270  // restrict this new phi to a coarser grid and store it in 'm_correction'
271 
272  T scratch;
273  this->m_op[depth]->create(scratch, a_phi);
274  this->m_op[depth]->setToZero(scratch);
275 
276  this->m_op[depth]->createCoarser(phi_coarse, a_phi, true); // create coarser grids for phi
277 
278  this->m_op[depth ]->restrictR(phi_coarse, a_phi);
279 
280  // Construct modified RHS
281 
282  // Calculate Op(Restrict(phi_fine))
283  this->m_op[depth+1]->applyOpMg(op_coarse, phi_coarse, m_phiCoarse[depth+1], false); // false - not homogeneous
284 
285  // This sets rhs = Restrict(residual fine)
286  this->m_op[depth ]->restrictResidual(*(this->m_residual[depth+1]), a_phi, m_phiCoarse[depth], a_rhs, false); // false - not homogeneous
287 
288  // This sets rhs = Restrict(residual fine) + Op(Restrict(phi fine))
289  this->m_op[depth+1]->incr(*(this->m_residual[depth+1]), op_coarse, 1.0);
290 
291  for (int img = 0; img < cycles; img++)
292  {
293  cycle(depth+1, phi_coarse, *(this->m_residual[depth+1]));
294  }
295 
296  // Interpolate the correction from the coarser level and add to phi on this level
297  this->m_op[depth ]->prolongIncrement(a_phi, *(this->m_correction[depth+1]));
298 
299  // Do some more smoothing of phi on this level
300  this->m_op[depth ]->relaxNF(a_phi, m_phiCoarse[depth], a_rhs, this->m_post);
301 
302  // Finally, compute the 'correction' (initial phi on this level - current phi on this level)
303  // so we can interpolate it when we get back to finer levels
304  this->m_op[depth]->incr(*(this->m_correction[depth]), a_phi, 1.0);
305 
306  // Clear temporary data objects
307  this->m_op[depth+1]->clear(phi_coarse);
308  this->m_op[depth+1]->clear(op_coarse);
309  }
310  }
311 }
312 
313 
314 #include "NamespaceFooter.H"
315 #endif /*_FASMULTIGRID_H_*/
std::ostream & pout()
Use this in place of std::cout for program output.
#define CH_assert(cond)
Definition: CHArray.H:37
A class to facilitate interaction with physical boundary conditions.
Definition: ProblemDomain.H:141
virtual void solve(T &a_phi, const T &a_rhs, Real a_tolerance, int a_maxIterations, int verbosity=0)
Definition: FASMultiGrid.H:107
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
const ProblemDomain & physDomain() const
FASMultiGrid()
Definition: FASMultiGrid.H:77
virtual void oneCycle(T &a_e, const T &a_res)
Execute ONE v-cycle of multigrid.
Definition: FASMultiGrid.H:157
#define CH_TIME(name)
Definition: CH_Timer.H:82
Definition: MultiGrid.H:294
double Real
Definition: REAL.H:33
Definition: MultiGrid.H:30
const IntVect & bigEnd() const
Definition: Box.H:1768
virtual ~FASMultiGrid()
Definition: FASMultiGrid.H:85
Definition: FASMultiGrid.H:22
A BoxLayout that has a concept of disjointedness.
Definition: DisjointBoxLayout.H:30
const Box & domainBox() const
Returns the logical computational domain.
Definition: ProblemDomain.H:887
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.
Definition: MultiGrid.H:331
virtual void cycle(int a_depth, T &a_correction, const T &a_residual)
Definition: FASMultiGrid.H:221
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:465
Definition: LinearSolver.H:156
T Max(const T &a_a, const T &a_b)
Definition: Misc.H:39
const IntVect & smallEnd() const
{ Accessors}
Definition: Box.H:1754
Vector< T * > m_phiCoarse
Definition: FASMultiGrid.H:68
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 FASMultiGrid object.
Definition: FASMultiGrid.H:91