Chombo + EB + MF  3.2
AMRFASMultiGrid.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 _AMRFASMULTIGRID_H_
16 #define _AMRFASMULTIGRID_H_
17 
18 #include "AMRMultiGrid.H"
19 #include "FASMultiGrid.H"
20 
21 #include "NamespaceHeader.H"
22 
23 
25 
26 ///
27 /**
28  Class to solve elliptic equations using the FAS multigrid
29  */
30 template <class T>
31 class AMRFASMultiGrid : public AMRMultiGrid<T>
32 {
33 public:
35  virtual ~AMRFASMultiGrid();
36 
37  virtual void define(const ProblemDomain& a_coarseDomain,
38  AMRLevelOpFactory<T>& a_factory,
39  LinearSolver<T>* a_bottomSolver,
40  int a_numLevels);
41 
42  virtual void solveNoInit(Vector<T*>& a_phi, const Vector<T*>& a_rhs,
43  int l_max, int l_base, bool a_zeroPhi=true,
44  bool forceHomogeneous = false);
45 
47  {
48  m_type = a_type;
49  }
50  void setAvoidNorms( bool b = true )
51  {
52  m_avoid_norms = b;
53  }
54  void setNumVcycles( int n )
55  {
56  m_numVcycles = n;
57  }
58 
59  // Number of multigrid iterations taken
60  int m_iter;
61 
62  // initial rnorm
64 
65  // current rnorm
67 
68 private:
69  virtual void FMG(Vector<T*>& a_phi,
70  const Vector<T*>& a_rhs,
71  int l, int l_max, int l_base);
72 
73  virtual void VCycle(Vector<T*>& a_phi,
74  const Vector<T*>& a_rhs,
75  int l, int l_max, int l_base);
76 
77  virtual void init(const Vector<T*>& a_phi, const Vector<T*>& a_rhs,
78  int l_max, int l_base);
79 
80  void clear_private();
81 
82  // data
84  bool m_avoid_norms; // flag to avoid norms and residuals (for convergence checking)
88 
90  ProblemDomain m_coarseDomain; // need to cache this
91 };
92 
93 //*******************************************************
94 // AMRFASMultigrid Implementation
95 //*******************************************************
96 
97 ////////////////////////////////////////////////////////////////////////
98 // AMRFASMultiGrid<T>::AMRFASMultiGrid
99 ////////////////////////////////////////////////////////////////////////
100 template <class T>
102  : AMRMultiGrid<T>()
103 {
104  m_numVcycles = 1;
105  m_type = VCYCLE;
106  m_avoid_norms = false;
107 }
108 
109 ////////////////////////////////////////////////////////////////////////
110 // AMRFASMultiGrid<T>::~AMRFASMultiGrid
111 ////////////////////////////////////////////////////////////////////////
112 template <class T>
114 {
115  CH_TIME("~AMRFASMultiGrid");
116  clear_private();
117 }
118 
119 ////////////////////////////////////////////////////////////////////////
120 // AMRFASMultiGrid<T>::define
121 ////////////////////////////////////////////////////////////////////////
122 template <class T>
123 void AMRFASMultiGrid<T>::define( const ProblemDomain& a_coarseDomain,
124  AMRLevelOpFactory<T>& a_factory,
125  LinearSolver<T>* a_bottomSolver,
126  int a_maxAMRLevels )
127 {
128  CH_TIME("AMRFASMultiGrid::define");
129 
130  AMRMultiGrid<T>::define( a_coarseDomain,a_factory,a_bottomSolver,a_maxAMRLevels);
131 
132  m_HOCopier.resize( a_maxAMRLevels );
133  m_revHOCopier.resize( a_maxAMRLevels );
134  m_HOCornerCopier.resize( a_maxAMRLevels );
135  m_coarseDomain = a_coarseDomain;
136 
137  // We actually want the MG solver to be FASMultigrid, not the default MultiGrid
138  // This is so we can solve down to the very coarsest grid in a nonlinear way
139  ProblemDomain current = a_coarseDomain;
140  for (int i = 0; i < a_maxAMRLevels; i++)
141  {
142 
143  // first delete initially defined standard MG object
144  if (this->m_mg[i] != NULL)
145  {
146  delete this->m_mg[i];
147  this->m_mg[i] = NULL;
148  }
149 
150  // Create FASMultigrid object instead
151  this->m_mg[i]= new FASMultiGrid<T>();
152  this->m_mg[i]->define(a_factory, &this->m_nosolve, current, this->m_maxDepth, this->m_op[i]);
153 
154  // Only do this if it will be used (avoiding a reference to invalid
155  // and/or unavailable refinement ratios)
156  if (i < a_maxAMRLevels-1)
157  {
158  current.refine(a_factory.refToFiner(current));
159  }
160  }
161 }
162 
163 ////////////////////////////////////////////////////////////////////////
164 // AMRFASMultiGrid<T>::init
165 ////////////////////////////////////////////////////////////////////////
166 template <class T>
167 void AMRFASMultiGrid<T>::init( const Vector<T*>& a_phi, const Vector<T*>& a_rhs,
168  int l_max, int l_base )
169 {
170  CH_TIME("AMRFASMultiGrid::init");
171 
172  AMRMultiGrid<T>::init( a_phi, a_rhs, l_max,l_base );
173 
174  // set new copiers for high order interp
175  // why were we doing this?
176 
178  for (int i = l_base+1; i <= l_max; i++)
179  {
180  AMRLevelOp<T>& op = *(this->m_op[i]);
181  int r = op.refToCoarser();
182  m_HOCopier[i].define( a_phi[i-1]->disjointBoxLayout(), this->m_resC[i]->disjointBoxLayout(),
183  a_phi[i-1]->ghostVect() );
184 
185  m_revHOCopier[i] = m_HOCopier[i];
186  m_revHOCopier[i].reverse();
187 
188 
189  m_HOCornerCopier[i].define( this->m_resC[i]->disjointBoxLayout(), this->m_resC[i]->disjointBoxLayout(), dom, this->m_resC[i]->ghostVect(), true ); //-- smart corner copier -- needed for AMR
190  //m_HOCornerCopier[i].define( this->m_resC[i]->disjointBoxLayout(), this->m_resC[i]->disjointBoxLayout(), this->m_resC[i]->ghostVect() ); -- dumb copier - not needed w/o AMR
191 
192  dom = dom.refine(r);
193  }
194 }
195 
196 ////////////////////////////////////////////////////////////////////////
197 // AMRFASMultiGrid<T>::clear_private
198 ////////////////////////////////////////////////////////////////////////
199 template <class T>
201 {
202  CH_TIME("AMRFASMultiGrid::clear_private");
203 
204  for (int i = 0; i < this->m_op.size(); i++)
205  {
206  }
207 }
208 
209 ////////////////////////////////////////////////////////////////////////
210 // AMRFASMultiGrid<T>::solveNoInit
211 ////////////////////////////////////////////////////////////////////////
212 template<class T>
214  const Vector<T*>& a_rhs,
215  int l_max, int l_base, bool a_zeroPhi,
216  bool a_forceHomogeneous)
217 {
218  CH_TIMERS("AMRFASMultiGrid::solveNoInit");
219  CH_TIMER("AMRFASMultiGrid::cycle", vtimer);
220 // CH_assert(!a_forceHomogeneous); // this is now allowed!
221  bool outputIntermediates = false;
222 
223  this->setBottomSolver( l_max, l_base );
224 
225  CH_assert(l_base <= l_max);
226  CH_assert(a_rhs.size() == a_phi.size());
227 
228  if (a_zeroPhi)
229  for (int ilev = l_base; ilev <=l_max; ++ilev)
230  {
231  this->m_op[ilev]->setToZero(*(a_phi[ilev]));
232  }
233 
234  m_initial_rnorm = 0;
235  if ( !m_avoid_norms )
236  {
237  CH_TIME("Initial AMR Residual");
238  m_initial_rnorm = this->computeAMRResidual( a_phi, a_rhs, l_max, l_base );
239  }
240 
241  if (this->m_convergenceMetric != 0.)
242  {
244  }
245 
246  // set bottom solver convergence norm and solver tolerance
247  this->m_bottomSolver->setConvergenceMetrics(m_initial_rnorm, this->m_bottomSolverEpsCushion*this->m_eps);
248 
250  Real norm_last = 2*m_initial_rnorm;
251 
252  int iter=0;
253  if (this->m_verbosity >= 2 && !m_avoid_norms )
254  {
255  // indentation in the output is to ensure that the :: lines up with the output from AMRMultigrid
256  pout() << " AMRFASMultiGrid:: iteration = " << iter << ", residual norm = " << m_rnorm << std::endl;
257  }
258 
259  bool goNorm = m_rnorm > this->m_normThresh || m_avoid_norms; //iterate if norm is not small enough
260  bool goRedu = m_rnorm > this->m_eps*m_initial_rnorm || m_avoid_norms;//iterate if initial norm is not reduced enough
261  bool goIter = iter < this->m_iterMax; //iterate if iter < max iteration count
262  bool goHang = iter < this->m_imin || m_rnorm <(1-this->m_hang)*norm_last; //iterate if we didn't hang
263  bool goMin = iter < this->m_iterMin ; // iterate if iter < min
264  while (goMin || (goIter && goRedu && goHang && goNorm))
265  {
266  norm_last = m_rnorm;
267 
268  //this generates a correction from the current residual
269  CH_START(vtimer);
270  if ( m_type == FULL )
271  {
272  FMG( a_phi, a_rhs, l_max, l_max, l_base);
273  }
274  else if ( m_type == VCYCLE )
275  {
276  // Set m_residual[l_max] = a_rhs[l_max]
277  this->m_op[l_max]->assignLocal( *(this->m_residual[l_max]), *(a_rhs[l_max]) );
278  VCycle( a_phi, a_rhs, l_max, l_max, l_base);
279  }
280  else
281  {
282  MayDay::Error("unknown FAS type");
283  }
284  CH_STOP(vtimer);
285  //increment phi by correction and reset correction to zero
286  for (int ilev = l_base; ilev <= l_max; ilev++)
287  {
288  if (outputIntermediates)
289  {
290  char strcorname[100];
291  sprintf(strcorname, "amrFASmg.phi.iter.%03d", iter);
292  string namecor(strcorname);
293  this->outputAMR(a_phi, namecor, l_max, l_base);
294  }
295  }
296 
297  // For solvers with accuracy higher than 2nd order
298  // consistency between levels has to be explicitly enforced.
299  if (this->m_op[0]->orderOfAccuracy()>2)
300  {
301  for (int ilev=l_max; ilev>l_base; ilev--)
302  {
303  this->m_op[ilev]->enforceCFConsistency(*a_phi[ilev-1], *a_phi[ilev]);
304  }
305  }
306 
307  // recompute residual
308  iter++;
309  if ( !m_avoid_norms )
310  {
311  m_rnorm = this->computeAMRResidual( a_phi, a_rhs, l_max, l_base );
312 
313  if (this->m_verbosity >= 2)
314  {
315  pout() << " AMRFASMultiGrid:: iteration = " << iter << ", residual norm = " << m_rnorm;
316  if (m_rnorm > 0.0)
317  {
318  pout() << ", rate = " << norm_last/m_rnorm;
319  }
320  pout() << std::endl;
321  }
322 
323  goNorm = m_rnorm > this->m_normThresh; //keep iterating if norm is not small enough
324  goRedu = m_rnorm > this->m_eps*m_initial_rnorm; //keep iterating if initial norm is not reduced enough
325  goHang = iter < this->m_imin || m_rnorm <(1-this->m_hang)*norm_last;//keep iterating if we didn't hang
326  }
327  goIter = iter < this->m_iterMax; //keep iterating if iter < max iteration count
328  goMin = iter < this->m_iterMin ; // keep iterating if iter < min
329  }
330 
331  this->m_exitStatus = int(!goRedu) + int(!goIter)*2 + int(!goHang)*4 + int(!goNorm)*8;
332  if (this->m_verbosity >= 2 && !m_avoid_norms)
333  {
334  pout() << " AMRFASMultiGrid:: iteration = " << iter << ", residual norm = " << m_rnorm << std::endl;
335  }
336  if (this->m_verbosity > 1)
337  {
338  if (!goIter && goRedu && goNorm) // goRedu=T, goIter=F, goHang=?, goNorm=T
339  { // m_exitStatus == 0 + 2 + 0|4 + 0 = 2|6
340  pout() << " AMRFASMultiGrid:: WARNING: Exit because max iteration count exceeded" << std::endl;
341  }
342  if (!goHang && goRedu && goNorm) // goRedu=T, goIter=?, goHang=F, goNorm=T
343  { // m_exitStatus == 0 + 0|2 + 4 + 0 = 4|6
344  pout() << " AMRFASMultiGrid:: WARNING: Exit because of solver hang" << std::endl;
345  }
346  if (this->m_verbosity > 4)
347  {
348  pout() << " AMRFASMultiGrid:: exitStatus = " << this->m_exitStatus << std::endl;
349  }
350  }
351 
352  // Store number or iterations taken so we can access this information
353  // from other objects
354  this->m_iter = iter;
355 }
356 
357 ////////////////////////////////////////////////////////////////////////
358 // AMRFASMultiGrid<T>::FMG
359 ////////////////////////////////////////////////////////////////////////
360 template<class T>
362  const Vector<T*>& a_rhs,
363  int i, int l_max, int l_base)
364 {
365  // coarse grid
366  this->m_op[l_base]->assignLocal( *(this->m_residual[l_base]), *(a_rhs[l_base]) );
367  this->m_mg[l_base]->oneCycle( *(a_phi[l_base]), *(this->m_residual[l_base]) );
368 
369  for ( int ilev = l_base+1 ; ilev <= l_max ; ilev++ )
370  {
371  // FMG interpolation -- high order
372  this->m_op[ilev]->AMRProlongS_2( *(a_phi[ilev]), *(a_phi[ilev-1]),
373  *(this->m_resC[ilev]), this->m_revHOCopier[ilev],
374  this->m_HOCornerCopier[ilev], this->m_op[ilev-1] );
375  // V-cycle
376  this->m_op[ilev]->assignLocal( *(this->m_residual[ilev]), *(a_rhs[ilev]) );
377  for (int i=0;i<m_numVcycles;i++)
378  {
379  VCycle( a_phi, a_rhs, ilev, l_max, l_base );
380  }
381  }
382 }
383 
384 ////////////////////////////////////////////////////////////////////////
385 // AMRFASMultiGrid<T>::VCycle
386 // 'm_residual' is really the full RHS ala FAS.
387 // 'm_correction' is just used as a temp
388 ////////////////////////////////////////////////////////////////////////
389 //#define FAS_PRINT
390 template<class T>
392  const Vector<T*>& a_rhs_dummy,
393  int ilev, int l_max, int l_base )
394 {
395  if (ilev == l_base)
396  {
397  CH_TIME("AMRFASMultiGrid::VCycle:coarse_grid_solver");
398  // exact solver
399 
400 #ifdef FAS_PRINT
401  this->m_op[ilev]->write(a_phi[ilev],"zh_coarsest_phi_0.hdf5");
402  this->m_op[ilev]->write(this->m_residual[ilev],"zh_coarsest_rhs.hdf5");
403 #endif
404 
405  // Keep doing FAS solves on the coarser levels
406  // Remember that m_residual is just the full rhs with FAS
407  // Cast as FASMultiGrid so we can use a different onyCycle
408  T* phiCoarsePtr = NULL;
409  if (l_base > 0)
410  {
411  phiCoarsePtr = a_phi[l_base-1];
412  }
413  FASMultiGrid<T>* fasMg = dynamic_cast<FASMultiGrid<T>*> (this->m_mg[l_base]);
414 
415  // This gives us back the full solution (NOT a correction)
416  // m_residual[l_base] is the full FAS right hand side (with tau correction)
417  fasMg->oneCycle( *(a_phi[l_base]), *(this->m_residual[l_base]), phiCoarsePtr);
418 
419 #ifdef FAS_PRINT
420 this->m_op[ilev]->write(a_phi[ilev],"zh_coarsest_phi_1.hdf5");
421 #endif
422  }
423  else
424  {
425  //Note: this else clause is currently untested
426  //============= Downsweep ========================
427 
428  T* rhs = NULL;
429 
430  // on the finest level, the rhs is the full rhs
431  // on coarser levels, it is whatever we calculate and put in m_residual
432  // need to be careful about this, as the fine level m_residual get's overwritten
433  // whenever we compute the AMR residual (e.g. at the end of each V-cycle)
434  if (ilev == l_max)
435  {
436  rhs = a_rhs_dummy[ilev];
437  }
438  else
439  {
440  rhs = this->m_residual[ilev];
441  }
442 
443 
444 #ifdef FAS_PRINT
445 this->m_op[ilev]->write(this->m_residual[ilev],"za_fine_res_0.hdf5");
446 this->m_op[ilev]->write(a_phi[ilev],"zb_fine_phi_0.hdf5");
447 #endif
448  // Replacing relax with relaxNF
449  // this->m_op[ilev]->relax( *(a_phi[ilev]), *(this->m_residual[ilev]), this->m_pre );
450  this->m_op[ilev]->relaxNF( *(a_phi[ilev]), (a_phi[ilev-1]), *rhs, this->m_pre );
451 
452 #ifdef FAS_PRINT
453 this->m_op[ilev]->write(this->m_residual[ilev],"za_fine_res_1.hdf5");
454 this->m_op[ilev]->write(a_phi[ilev],"zb_fine_phi_1.hdf5");
455 this->m_op[ilev]->residual( *(this->m_correction[ilev]), *a_phi[ilev], *(this->m_residual[ilev]), false);
456 this->m_op[ilev]->write(this->m_correction[ilev],"zz_fine_res_0.hdf5");
457 #endif
458 
459 
460  // Compute residual on next coarser level
461  // for the valid region NOT covered by this level.
462  // includes reflux corrections
463  // do this before we fill a_phi[ilev-1] with junk
465  a_phi,
466  a_rhs_dummy,
467  l_max, l_base, ilev-1,
468  false);
469 
470 
471  // Compute the restriction of the solution to the coarser level phiC -- R(u_f)
472  // Need this so that later we can compute L_c(R(u_f))
473  // This just restricts a_phi[ilev] to m_resC[ilev] (due to the skip residuals flag)
474  this->m_op[ilev]->AMRRestrictS(*(this->m_resC[ilev]), // output
475  *(a_phi[ilev]), // input
476  *(a_phi[ilev]), // dummy
477  *(a_phi[ilev-1]), // dummy
478  *(this->m_correction[ilev]),// scratch
479  true ); // skip residuals
480 
481 #ifdef FAS_PRINT
482 this->m_op[ilev]->write(a_phi[ilev],"zb_fine_phi_2.hdf5");
483 this->m_op[ilev-1]->write(this->m_resC[ilev],"zb_coarse_phi_1.hdf5");
484 #endif
485  // Overwrite R(u_f) on the valid region of the next coarser level a_phi[ilev-1]
486  //this->m_op[ilev-1]->assignCopier( *(a_phi[ilev-1]), *(m_phiC[ilev]), m_phiCopier[ilev] );
487  this->m_resC[ilev]->copyTo(this->m_resC[ilev]->interval(), *(a_phi[ilev-1]), a_phi[ilev-1]->interval(), this->m_resCopier[ilev]);
488 
489 
490  // Compute the restriction of the _residual_ to the coarser level
491  // for the valid region covered by this level
492  // resC -- R(f - L_f(u_f))
493  this->m_op[ilev]->AMRRestrictS(*(this->m_resC[ilev]), // output
494  *rhs, // const input
495  *(a_phi[ilev]), // const but C-F interp modified
496  *(a_phi[ilev-1]), // coarse phi, for C-F interp.
497  *(this->m_correction[ilev])); // scratch
498 
499 
500 #ifdef FAS_PRINT
501 this->m_op[ilev]->write(this->m_resC[ilev],"zd_resC.hdf5");
502 this->m_op[ilev]->write(a_phi[ilev-1],"ze_coarse_phi_2.hdf5");
503 this->m_op[ilev]->write(this->m_residual[ilev],"za_fine_res_2.hdf5");
504 #endif
505 
506  // Overwrite residual on the valid region of the next coarser level
507  // with coarsened residual from this level
508  //this->m_op[ilev-1]->assignCopier( *(this->m_residual[ilev-1]), *(this->m_resC[ilev]), this->m_resCopier[ilev]);
509  this->m_resC[ilev]->copyTo(this->m_resC[ilev]->interval(), *(this->m_residual[ilev-1]), this->m_residual[ilev-1]->interval(), this->m_resCopier[ilev] );
510 
511 #ifdef FAS_PRINT
512 this->m_op[ilev-1]->write(this->m_residual[ilev-1],"zf_coarse_res_1.hdf5");
513 #endif
514 
515  // Compute the tau correction on the covered region, L_c (R(u_f)), and put it in m_correction[ilev-1]
516  // Do not pass in fine level here - we've already done refluxing
517  // this->m_op[ilev-1]->AMROperatorNC( *(this->m_correction[ilev-1]), *(a_phi[ilev]), *(a_phi[ilev-1]), true, this->m_op[ilev]);
518  T undefined;
519  this->m_op[ilev-1]->AMROperatorNC( *(this->m_correction[ilev-1]), undefined, *(a_phi[ilev-1]), false, this->m_op[ilev]);
520 
521 #ifdef FAS_PRINT
522 this->m_op[ilev-1]->write(this->m_correction[ilev-1],"zf_coarse_Lu.hdf5");
523 this->m_op[ilev]->write(a_phi[ilev],"zb_fine_phi_3.hdf5");
524 #endif
525 
526  // compute the full RHS, R( f - L_f(u_f) ) + L_c( R(u_f) )
527  this->m_op[ilev-1]->axby( *(this->m_residual[ilev-1]), *(this->m_residual[ilev-1]), *(this->m_correction[ilev-1]), 1.0, 1.0);
528 
529 #ifdef FAS_PRINT
530 this->m_op[ilev-1]->write(this->m_residual[ilev-1],"zf_coarse_res_3.hdf5");
531 #endif
532 
533  {
534  // store correction in R_u_f -- R(u_f)
535  T R_u_f;
536  this->m_op[ilev-1]->create( R_u_f, *a_phi[ilev-1]);
537  this->m_op[ilev-1]->assignLocal( R_u_f, *(a_phi[ilev-1]) );
538  //============finish Compute residual for the next coarser level======
539  for (int img = 0; img < this->m_numMG; img++)
540  {
541  VCycle( a_phi, a_rhs_dummy, ilev-1, l_max, l_base );
542  }
543 
544 #ifdef FAS_PRINT
545 this->m_op[ilev-1]->write(a_phi[ilev-1],"zi_coarse_phi_2.hdf5");
546 #endif
547 
548  // subtract off initial solution to get an increment (m_correction[ilev-1] is an increment)
549  this->m_op[ilev-1]->axby( *this->m_correction[ilev-1], *a_phi[ilev-1], R_u_f, 1.0, -1.0 );
550  this->m_op[ilev-1]->clear( R_u_f );
551  }
552 
553 #ifdef FAS_PRINT
554 this->m_op[ilev-1]->residual( *(this->m_correction[ilev-1]), *a_phi[ilev-1], *(this->m_residual[ilev-1]), true);
555 this->m_op[ilev-1]->write(this->m_correction[ilev-1],"zz_coarse_res.hdf5");
556 #endif
557 
558  //================= Upsweep ======================
559  // increment the correction with coarser version
560 
561 // There's a bug with the copiers and AMRProlongS_2, so doing this the old way for now
562 // this->m_op[ilev]->AMRProlongS_2( *(a_phi[ilev]), *(a_phi[ilev-1]),
563 // *(this->m_resC[ilev]), this->m_HOCopier[ilev],
564 // this->m_HOCornerCopier[ilev], this->m_op[ilev-1] );
565 
566 
567  this->m_op[ilev]->AMRProlongS( *(a_phi[ilev]), *(this->m_correction[ilev-1]),
568  *(this->m_resC[ilev]),
569  this->m_reverseCopier[ilev] );
570 
571 
572 #ifdef FAS_PRINT
573 this->m_op[ilev]->write(a_phi[ilev],"zl_fine_phi_1.hdf5");
574 this->m_op[ilev]->residual( *(this->m_correction[ilev]), *a_phi[ilev], *(this->m_residual[ilev]), false);
575 this->m_op[ilev]->write(this->m_correction[ilev],"zz_fine_res_1.hdf5");
576 #endif
577  // Replacing relax with relaxNF
578 // this->m_op[ilev]->relax( *(a_phi[ilev]), *(this->m_residual[ilev]), this->m_post );
579  this->m_op[ilev]->relaxNF( *(a_phi[ilev]), (a_phi[ilev-1]), *rhs, this->m_post );
580 
581 #ifdef FAS_PRINT
582 this->m_op[ilev]->write(a_phi[ilev],"zl_fine_phi_2.hdf5");
583 this->m_op[ilev]->write(this->m_residual[ilev],"za_fine_res_3.hdf5");
584 this->m_op[ilev]->write(a_rhs_dummy[ilev],"za_rhs.hdf5");
585 this->m_op[ilev]->residual( *(this->m_correction[ilev]), *a_phi[ilev], *(this->m_residual[ilev]), false);
586 this->m_op[ilev]->write(this->m_correction[ilev],"zz_fine_res_2.hdf5");
587 #endif
588  }
589 }
590 
591 #include "NamespaceFooter.H"
592 
593 #endif
std::ostream & pout()
Use this in place of std::cout for program output.
virtual int refToFiner(const ProblemDomain &a_indexSpace) const =0
#define CH_TIMERS(name)
Definition: CH_Timer.H:133
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: AMRFASMultiGrid.H:213
#define CH_assert(cond)
Definition: CHArray.H:37
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 void VCycle(Vector< T *> &a_phi, const Vector< T *> &a_rhs, int l, int l_max, int l_base)
Definition: AMRFASMultiGrid.H:391
int m_iterMax
Definition: AMRMultiGrid.H:380
LinearSolver< T > * m_bottomSolver
Definition: AMRMultiGrid.H:492
int m_verbosity
Definition: AMRMultiGrid.H:380
OLD_FASMG_type
Definition: AMRFASMultiGrid.H:24
#define CH_START(tpointer)
Definition: CH_Timer.H:145
Vector< T * > m_residual
Definition: AMRMultiGrid.H:485
Definition: AMRFASMultiGrid.H:24
Definition: AMRFASMultiGrid.H:24
int m_pre
Definition: AMRMultiGrid.H:381
Vector< Copier > m_reverseCopier
Definition: AMRMultiGrid.H:488
virtual void define(const ProblemDomain &a_coarseDomain, AMRLevelOpFactory< T > &a_factory, LinearSolver< T > *a_bottomSolver, int a_numLevels)
Definition: AMRFASMultiGrid.H:123
void setCycleType(OLD_FASMG_type a_type)
Definition: AMRFASMultiGrid.H:46
Real m_eps
Definition: AMRMultiGrid.H:378
Definition: AMRMultiGrid.H:308
ProblemDomain m_coarseDomain
Definition: AMRFASMultiGrid.H:90
OLD_FASMG_type m_type
Definition: AMRFASMultiGrid.H:83
Vector< MultiGrid< T > * > m_mg
Definition: AMRMultiGrid.H:483
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
Vector< Copier > m_revHOCopier
Definition: AMRFASMultiGrid.H:87
Real m_convergenceMetric
Definition: AMRMultiGrid.H:391
Definition: AMRMultiGrid.H:39
#define CH_TIMER(name, tpointer)
Definition: CH_Timer.H:63
Real m_initial_rnorm
Definition: AMRFASMultiGrid.H:63
int m_numVcycles
Definition: AMRFASMultiGrid.H:85
void resize(unsigned int isize)
Definition: Vector.H:346
virtual void oneCycle(T &a_e, const T &a_res)
Execute ONE v-cycle of multigrid.
Definition: FASMultiGrid.H:157
int m_numMG
Definition: AMRMultiGrid.H:381
Vector< Copier > m_HOCornerCopier
Definition: AMRFASMultiGrid.H:89
void setAvoidNorms(bool b=true)
Definition: AMRFASMultiGrid.H:50
void clear()
Definition: Vector.H:180
virtual int refToCoarser()=0
#define CH_TIME(name)
Definition: CH_Timer.H:82
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
Vector< Copier > m_HOCopier
Definition: AMRFASMultiGrid.H:86
int m_iter
Definition: AMRFASMultiGrid.H:60
Definition: FASMultiGrid.H:22
Real m_hang
Definition: AMRMultiGrid.H:378
Real m_rnorm
Definition: AMRFASMultiGrid.H:66
int m_iterMin
Definition: AMRMultiGrid.H:380
size_t size() const
Definition: Vector.H:192
Vector< T * > m_correction
Definition: AMRMultiGrid.H:484
Vector< T * > m_resC
Definition: AMRMultiGrid.H:486
AMRFASMultiGrid()
Definition: AMRFASMultiGrid.H:101
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: AMRFASMultiGrid.H:24
void outputAMR(Vector< T *> &a_data, string &a_name, int a_lmax, int a_lbase)
Definition: AMRMultiGrid.H:516
int m_exitStatus
Definition: AMRMultiGrid.H:380
ProblemDomain & refine(int a_refinement_ratio)
Refine this problem domain.
bool m_avoid_norms
Definition: AMRFASMultiGrid.H:84
Definition: LinearSolver.H:156
int m_post
Definition: AMRMultiGrid.H:381
int m_imin
Definition: AMRMultiGrid.H:380
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
#define CH_STOP(tpointer)
Definition: CH_Timer.H:150
virtual void init(const Vector< T *> &a_phi, const Vector< T *> &a_rhs, int l_max, int l_base)
Definition: AMRMultiGrid.H:1212
Vector< Copier > m_resCopier
Definition: AMRMultiGrid.H:487
Definition: AMRFASMultiGrid.H:31
virtual ~AMRFASMultiGrid()
Definition: AMRFASMultiGrid.H:113
Definition: AMRMultiGrid.H:233
virtual void init(const Vector< T *> &a_phi, const Vector< T *> &a_rhs, int l_max, int l_base)
Definition: AMRFASMultiGrid.H:167
NoOpSolver< T > m_nosolve
Definition: AMRMultiGrid.H:490
void setNumVcycles(int n)
Definition: AMRFASMultiGrid.H:54
Real m_normThresh
Definition: AMRMultiGrid.H:378
virtual void FMG(Vector< T *> &a_phi, const Vector< T *> &a_rhs, int l, int l_max, int l_base)
Definition: AMRFASMultiGrid.H:361
void clear_private()
Definition: AMRFASMultiGrid.H:200
Vector< AMRLevelOp< T > * > m_op
Definition: AMRMultiGrid.H:482
void setBottomSolver(int l_max, int l_base)
set up bottom solver for internal MG solver
Definition: AMRMultiGrid.H:1266