Bug in multigrid removed: the old way of restricting the residual was wrong

git-svn-id: https://svn.dealii.org/trunk@3324 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/deal.II/include/multigrid/mg_base.h b/deal.II/deal.II/include/multigrid/mg_base.h
index 71a41e963ec65c4b9df3d5f4897245fdf6c658a1..3eeabaa848afd0b75cbbe5b42cc97123751a7d6f 100644 (file)
--- a/deal.II/deal.II/include/multigrid/mg_base.h
+++ b/deal.II/deal.II/include/multigrid/mg_base.h
@@ -356,6 +356,11 @@ class MGBase : public Subscriptor
                                      */
     MGLevelObject<Vector<double> > solution;
     
+                                    /**
+                                     * Auxiliary vector.
+                                     */
+    MGLevelObject<Vector<double> > t;    
+
                                     /**
                                      * Prolongation and restriction object.
                                      */
@@ -431,11 +436,6 @@ class MGBase : public Subscriptor
                               const Vector<double>& v,
                               const char* name) const = 0;
     
-  private:
-                                    /**
-                                     * Auxiliary vector.
-                                     */
-    Vector<double> t;    
 };
 
 

diff --git a/deal.II/deal.II/include/multigrid/multigrid.templates.h b/deal.II/deal.II/include/multigrid/multigrid.templates.h
index 59eb3f701e11684bd78f16c5ae554928d552e1a4..8ff60ac22fd218db09018bdea44e32f5cb957ab0 100644 (file)
--- a/deal.II/deal.II/include/multigrid/multigrid.templates.h
+++ b/deal.II/deal.II/include/multigrid/multigrid.templates.h
@@ -64,6 +64,7 @@ Multigrid<dim>::copy_to_mg (const Vector<number>& osrc)
       
       solution[level].reinit(system_size);
       defect[level].reinit(system_size);
+      t[level].reinit(system_size);
     };
 
                                   // traverse the grid top-down
@@ -121,14 +122,15 @@ Multigrid<dim>::copy_to_mg (const Vector<number>& osrc)
                };
            };
        };
-
                                       // for that part of the level
                                       // which is further refined:
                                       // get the defect by
                                       // restriction of the defect on
                                       // one level higher
       if (static_cast<unsigned int>(level) < maxlevel)
-       transfer->restrict_and_add (level+1, defect[level], defect[level+1]);
+       {
+         transfer->restrict_and_add (level+1, defect[level], defect[level+1]);
+       }
     };
 };
 

diff --git a/deal.II/deal.II/source/multigrid/mg_base.cc b/deal.II/deal.II/source/multigrid/mg_base.cc
index f1ab40ffe2927ae76930b7234792dc23bd8cc817..516e607d2ac176601f5b041ffd8d6270b8bbdaa4 100644 (file)
--- a/deal.II/deal.II/source/multigrid/mg_base.cc
+++ b/deal.II/deal.II/source/multigrid/mg_base.cc
@@ -30,6 +30,7 @@ MGBase::MGBase(const MGTransferBase &transfer,
                minlevel(minlevel),
                defect(minlevel,maxlevel),
                solution(minlevel,maxlevel),
+               t(minlevel,maxlevel),
                transfer(&transfer)
 {
   Assert(minlevel <= maxlevel,
@@ -73,24 +74,31 @@ MGBase::level_mgstep(const unsigned int        level,
 
                           // smoothing of the residual by modifying s
   pre_smooth.smooth(level, solution[level], defect[level]);
-                                  // t = d-As
 
 #ifdef MG_DEBUG
   sprintf(name, "MG%d-pre",level);
   print_vector(level, solution[level], name);
 #endif
   
-  t.reinit(solution[level].size());
-  level_vmult(level, t, solution[level], defect[level]);
+  t[level] = 0.;
+
+                                  // t = -A*solution[level]
+  level_vmult(level, t[level], solution[level], defect[level]);
   
                                   // make t rhs of lower level
                                   // The non-refined parts of the
                                   // coarse-level defect already contain
-                                  // the global defect.
-  transfer->restrict_and_add (level, defect[level-1], t);
+                                  // the global defect, the refined parts
+                                  // its restriction.
+  for (unsigned int l = level;l>minlevel;--l)
+    {
+      t[l-1] = 0.;
+      transfer->restrict_and_add (l, t[l-1], t[l]);
+      defect[l-1] += t[l-1];
+    }
 
                                   // add additional DG contribution
-  edge_vmult(level, defect[level-1], defect[level]);
+//  edge_vmult(level, defect[level-1], defect[level]);
   
                                   // do recursion
   level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver);
@@ -99,21 +107,20 @@ MGBase::level_mgstep(const unsigned int        level,
                                   // vector, since it has been
                                   // resized in the recursive call to
                                   // level_mgstep directly above
-  t.reinit(solution[level].size());
+  t[level] = 0.;
 
                                   // do coarse grid correction
 
-  transfer->prolongate(level, t, solution[level-1]);
+  transfer->prolongate(level, t[level], solution[level-1]);
 
 #ifdef MG_DEBUG
   sprintf(name, "MG%d-cgc",level);
   print_vector(level, t, name);
 #endif
 
-  solution[level] += t;
+  solution[level] += t[level];
   
                                   // post-smoothing
-
   post_smooth.smooth(level, solution[level], defect[level]);
 
 #ifdef MG_DEBUG

diff --git a/deal.II/deal.II/source/multigrid/multigrid.cc b/deal.II/deal.II/source/multigrid/multigrid.cc
index a6878404c0726356867bc1df3438abd5a9e3a205..652b36ad0d5bc2f1a35908f02a83ff6aeca4282c 100644 (file)
--- a/deal.II/deal.II/source/multigrid/multigrid.cc
+++ b/deal.II/deal.II/source/multigrid/multigrid.cc
@@ -134,7 +134,6 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler<dim> &mg_dof)
       prolongation_matrices.push_back (SparseMatrix<double>());
       prolongation_matrices[level].reinit (prolongation_sparsities[level]);
 
-
 // now actually build the matrices
       for (MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
           cell != mg_dof.end(level); ++cell)

diff --git a/tests/deal.II/mg.checked b/tests/deal.II/mg.checked
index 43091d5a02213ac86875c6a0937f688e3565be02..272ed178a21306afb26eae9b3b2849637c108515 100644 (file)
--- a/tests/deal.II/mg.checked
+++ b/tests/deal.II/mg.checked
@@ -6,32 +6,32 @@ DEAL:Q1:cg::Convergence step 4
 DEAL:Q1::DoFs 25
 DEAL:Q1::Levels: 3
 DEAL:Q1:cg::Starting 
-DEAL:Q1:cg::Convergence step 7 
+DEAL:Q1:cg::Convergence step 5 
 DEAL:Q1::DoFs 81
 DEAL:Q1::Levels: 4
 DEAL:Q1:cg::Starting 
-DEAL:Q1:cg::Convergence step 9 
+DEAL:Q1:cg::Convergence step 6 
 DEAL:Q2::DoFs 25
 DEAL:Q2::Levels: 2
 DEAL:Q2:cg::Starting 
-DEAL:Q2:cg::Convergence step 6 
+DEAL:Q2:cg::Convergence step 5 
 DEAL:Q2::DoFs 81
 DEAL:Q2::Levels: 3
 DEAL:Q2:cg::Starting 
-DEAL:Q2:cg::Convergence step 10 
+DEAL:Q2:cg::Convergence step 7 
 DEAL:Q2::DoFs 289
 DEAL:Q2::Levels: 4
 DEAL:Q2:cg::Starting 
-DEAL:Q2:cg::Convergence step 9 
+DEAL:Q2:cg::Convergence step 7 
 DEAL:Q3::DoFs 49
 DEAL:Q3::Levels: 2
 DEAL:Q3:cg::Starting 
-DEAL:Q3:cg::Convergence step 7 
+DEAL:Q3:cg::Convergence step 8 
 DEAL:Q3::DoFs 169
 DEAL:Q3::Levels: 3
 DEAL:Q3:cg::Starting 
-DEAL:Q3:cg::Convergence step 10 
+DEAL:Q3:cg::Convergence step 7 
 DEAL:Q3::DoFs 625
 DEAL:Q3::Levels: 4
 DEAL:Q3:cg::Starting 
-DEAL:Q3:cg::Convergence step 11 
+DEAL:Q3:cg::Convergence step 7