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