Multigrid extension for NeuroFEM 
  last change: 04.03.2008  

New Solver - Full Multigrid

The following code shows the full multigrid implementation. Due to the fact that we are only using V-Cycles, the recursion is unrolled in our code to some loops. Additional the functionality of ParExPDE allows an easy writing of code for doing e.g. Gauss-Seidel smoothing steps. The outermost loop represents the initialization phase: we are starting on the coarsest level, solve the problem with multigrid and then proceed on the next finer level. Inside this loop a number of V-Cycle iterations is done, either til convergence or til the maximum number of cycles is reached. Inside each of those iterations a single V-Cycle is performed. This means that we start on the fine grid and proceed down to the coarest grid with some presmoothing steps in between. On the coarsest grid the problem is solved with a high number of Gauss-Seidel operations. Afterwards we proceed upwars again to the finest level including some postsmoothings. Finally in each iteration the residual is calculated. When taking a closer look at the following piece of code, you can find all those parts inside.
Code from file "neurofem/fmg/fmgsolver.cpp"
void FMGSolver::SolveProblem() {
  // *** Running the multigrid (V cycle) ***
  for (int level = 0; level < number_levels-1; ++level) {
    (*f).doRestrict();
    (*u).doRestrict();
    (*r).levelDown();
  }
  *u = 0.0;
  for (int lvl = number_levels-1; lvl >= 0; lvl--) {
    double resnorm = eps+1;
    int smooth;
    int i = 1;
    while ((i <= vcycle[lvl]) && (eps < resnorm)) {
      i++;
      // Downwards
      for (int level = lvl; level < number_levels-1; ++level) {
        // Presmoothing
        for (int smooth = 0; smooth < presmooth; ++smooth) {
          *u = *u + (*f - op(*u)) / ParExPDE::Diag(op) | ParExPDE::interior_points;
          counter[level]++;
        }
        // Restriction
        *r = *f - op(*u) | ParExPDE::interior_points;
        *r = 0.0 | ParExPDE::boundary_points;
        (*r).doRestrict();
        (*f).levelDown();
        *f = *r;
        (*u).levelDown();
        *u = 0.0;
      }
      // Solve on coarsest grid
      smooth = 0;
      *r = *f - op(*u) | ParExPDE::interior_points;
      double resnorm2 = sqrt(ParExPDE::ScalarProduct(*r,*r,&ParExPDE::interior_points) / normi[number_levels-1]);
      while ((smooth < cgridcorrections) && (eps < resnorm2)) {
        smooth++;
        *u = *u + (*f - op(*u)) / ParExPDE::Diag(op) | ParExPDE::interior_points;
        *r = *f - op(*u) | ParExPDE::interior_points;
        resnorm2 = sqrt(ParExPDE::ScalarProduct(*r,*r,&ParExPDE::interior_points) / normi[number_levels-1]);
        counter[number_levels-1]++;
      }
      // Upwards
      for (int level = number_levels - 2; level >= lvl; level--) {
        // Prolongation
        (*f).levelUp();
        *r = *u;
        (*r).doProlongate();
        (*u).levelUp();
        *u = *u + *r | ParExPDE::interior_points;
        // Postsmoothing
        for (int smooth = 0; smooth < postsmooth; ++smooth) {
          *u = *u + (*f - op(*u)) / ParExPDE::Diag(op) | ParExPDE::interior_points;
          counter[level]++;
        }
      }
      // Calculate residual
      *r = *f - op(*u) | ParExPDE::interior_points;
      resnorm = sqrt(ParExPDE::ScalarProduct(*r,*r,&ParExPDE::interior_points) / normi[0]);
    }
    // Go one level up
    if (lvl != 0) {
      (*f).levelUp();
      (*u).doProlongate();
      (*r).levelUp();
    }
  }
}