From eca9b77619f1357613a2000c43df3b11b3906dbd Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 16 Sep 2016 12:09:50 +0200 Subject: [PATCH] Work on comments. Fix and test W and F cycles. --- .../deal.II/multigrid/multigrid.templates.h | 167 +++--- tests/multigrid/cycles.output | 494 ++++++++++------- tests/multigrid/step-16-05.cc | 509 ++++++++++++++++++ tests/multigrid/step-16-05.output | 57 ++ tests/multigrid/step-16-06.cc | 509 ++++++++++++++++++ tests/multigrid/step-16-06.output | 57 ++ 6 files changed, 1489 insertions(+), 304 deletions(-) create mode 100644 tests/multigrid/step-16-05.cc create mode 100644 tests/multigrid/step-16-05.output create mode 100644 tests/multigrid/step-16-06.cc create mode 100644 tests/multigrid/step-16-06.output diff --git a/include/deal.II/multigrid/multigrid.templates.h b/include/deal.II/multigrid/multigrid.templates.h index 8400c63521..d5694e1cb2 100644 --- a/include/deal.II/multigrid/multigrid.templates.h +++ b/include/deal.II/multigrid/multigrid.templates.h @@ -121,71 +121,61 @@ Multigrid::level_v_step (const unsigned int level) (*coarse)(level, solution[level], defect[level]); return; } + + // smoothing of the residual if (debug>1) deallog << "Smoothing on level " << level << std::endl; - // smoothing of the residual by - // modifying s -// defect[level].print(std::cout, 2,false); -// std::cout<smooth(level, solution[level], defect[level]); -// solution[level].print(std::cout, 2,false); + // solution[level].print(std::cout, 2,false); if (debug>2) deallog << "Solution norm " << solution[level].l2_norm() << std::endl; + // compute residual on level, which includes the (CG) edge matrix if (debug>1) deallog << "Residual on level " << level << std::endl; - // t = A*solution[level] matrix->vmult(level, t[level], solution[level]); - - if (debug>2) - deallog << "Residual norm " << t[level].l2_norm() - << std::endl; -// std::cout<vmult_add(level, t[level], solution[level]); if (debug>2) deallog << "Norm t[" << level << "] " << t[level].l2_norm() << std::endl; } - t[level].sadd(-1.0, 1.0, defect[level]); - // transfer the residual to the coarser level + if (debug>2) + deallog << "Residual norm " << t[level].l2_norm() + << std::endl; + + // Get the defect on the next coarser level as part of the (DG) edge matrix + // and then the main part by the restriction of the transfer if (edge_down != 0) { edge_down->vmult(level, t[level-1], solution[level]); defect[level-1] -= t[level-1]; } - transfer->restrict_and_add(level, defect[level-1], t[level]); // do recursion - solution[level-1] = 0.; level_v_step(level-1); - // reset size of the auxiliary - // vector, since it has been - // resized in the recursive call to - // level_v_step directly above - t[level] = 0.; - // do coarse grid correction transfer->prolongate(level, t[level], solution[level-1]); if (debug>2) deallog << "Prolongate norm " << t[level].l2_norm() << std::endl; - solution[level] += t[level]; + // get in contribution from edge matrices to the defect if (edge_in != 0) { edge_in->Tvmult(level, t[level], solution[level]); defect[level] -= t[level]; } - if (edge_up != 0) { edge_up->Tvmult(level, t[level], solution[level-1]); @@ -196,15 +186,10 @@ Multigrid::level_v_step (const unsigned int level) deallog << "V-cycle Defect norm " << defect[level].l2_norm() << std::endl; + // post-smoothing if (debug>1) deallog << "Smoothing on level " << level << std::endl; - // post-smoothing - -// std::cout<smooth(level, solution[level], defect[level]); -// solution[level].print(std::cout, 2,false); -// std::cout<2) deallog << "Solution norm " << solution[level].l2_norm() @@ -236,112 +221,106 @@ Multigrid::level_step(const unsigned int level, } if (debug>0) - deallog << cychar << "-cycle entering level " << level << std::endl; + deallog << cychar << "-cycle entering level " << level << std::endl; - // Not actually the defect yet, but - // we do not want to spend yet - // another vector. - if (level>minlevel) - { - defect2[level-1] = 0.; - transfer->restrict_and_add (level, defect2[level-1], defect2[level]); - } - - // We get an update of the defect - // from the previous level in t and - // from two levels above in - // defect2. This is subtracted from - // the original defect. - t[level].equ(-1.,defect2[level]); - t[level] += defect[level]; + // Combine the defect from the initial copy_to_mg with the one that has come + // from the finer level by the transfer + defect2[level] += defect[level]; + defect[level] = 0; if (debug>2) - deallog << cychar << "-cycle defect norm " << t[level].l2_norm() + deallog << cychar << "-cycle defect norm " << defect2[level].l2_norm() << std::endl; if (level == minlevel) { if (debug>0) - deallog << cychar << "-cycle coarse level " << level << std::endl; + deallog << cychar << "-cycle coarse level " << level << std::endl; - (*coarse)(level, solution[level], t[level]); + (*coarse)(level, solution[level], defect2[level]); return; } + + // smoothing of the residual if (debug>1) deallog << cychar << "-cycle smoothing level " << level << std::endl; - // smoothing of the residual by - // modifying s - pre_smooth->smooth(level, solution[level], t[level]); + pre_smooth->smooth(level, solution[level], defect2[level]); if (debug>2) - deallog << cychar << "-cycle solution norm " + deallog << cychar << "-cycle solution norm " << solution[level].l2_norm() << std::endl; + // compute residual on level, which includes the (CG) edge matrix if (debug>1) - deallog << cychar << "-cycle residual level " << level << std::endl; - // t = A*solution[level] + deallog << cychar << "-cycle residual level " << level << std::endl; matrix->vmult(level, t[level], solution[level]); - // make t rhs of lower level The - // non-refined parts of the - // coarse-level defect already - // contain the global defect, the - // refined parts its restriction. if (edge_out != 0) edge_out->vmult_add(level, t[level], solution[level]); + t[level].sadd(-1.0, 1.0, defect2[level]); + if (debug>2) + deallog << cychar << "-cycle residual norm " << t[level].l2_norm() + << std::endl; + + // Get the defect on the next coarser level as part of the (DG) edge matrix + // and then the main part by the restriction of the transfer if (edge_down != 0) - edge_down->vmult_add(level, defect2[level-1], solution[level]); + edge_down->vmult(level, defect2[level-1], solution[level]); + else + defect2[level-1] = 0; transfer->restrict_and_add (level, defect2[level-1], t[level]); - // do recursion - solution[level-1] = 0.; - // Every cycle need one recursion, - // the V-cycle, which is included - // here for the sake of the - // F-cycle, needs only one, + + // Every cycle needs one recursion, the V-cycle, which is included here for + // the sake of the F-cycle, needs only one. As opposed to the V-cycle above + // where we do not need to set the solution vector to zero (as that has been + // done in the init functions), we also need to set the content in the + // solution from previous cycles to zero. + solution[level-1] = 0; level_step(level-1, cycle); - // If we solve exactly, then we do - // not need a second coarse grid - // step. + + // For W and F-cycle, repeat the process on the next coarser level except + // for the coarse solver which we invoke just once if (level>minlevel+1) { // while the W-cycle repeats itself if (cycle == w_cycle) level_step(level-1, cycle); - // and the F-cycle does a V-cycle - // after an F-cycle. + + // the F-cycle does a V-cycle after an F-cycle. else if (cycle == f_cycle) level_step(level-1, v_cycle); } - // reset size of the auxiliary - // vector, since it has been - // resized in the recursive call to - // level_v_step directly above - t[level] = 0.; // do coarse grid correction transfer->prolongate(level, t[level], solution[level-1]); - solution[level] += t[level]; - + // get in contribution from edge matrices to the defect if (edge_in != 0) - edge_in->Tvmult(level, t[level], solution[level]); + { + edge_in->Tvmult(level, t[level], solution[level]); + defect2[level] -= t[level]; + } if (edge_up != 0) - edge_up->Tvmult(level, t[level], solution[level-1]); - - t[level].sadd(-1.,-1.,defect2[level]); - t[level] += defect[level]; + { + edge_up->Tvmult(level, t[level], solution[level-1]); + defect2[level] -= t[level]; + } if (debug>2) - deallog << cychar << "-cycle Defect norm " << t[level].l2_norm() + deallog << cychar << "-cycle Defect norm " << defect2[level].l2_norm() << std::endl; + // post-smoothing if (debug>1) deallog << cychar << "-cycle smoothing level " << level << std::endl; - // post-smoothing - post_smooth->smooth(level, solution[level], t[level]); + post_smooth->smooth(level, solution[level], defect2[level]); + + if (debug>2) + deallog << cychar << "-cycle solution norm " << solution[level].l2_norm() + << std::endl; if (debug>1) deallog << cychar << "-cycle leaving level " << level << std::endl; @@ -353,9 +332,8 @@ template void Multigrid::cycle() { - // The defect vector has been - // initialized by copy_to_mg. Now - // adjust the other vectors. + // The defect vector has been initialized by copy_to_mg. Now adjust the + // other vectors. solution.resize(minlevel, maxlevel); t.resize(minlevel, maxlevel); if (cycle_type != v_cycle) @@ -381,9 +359,8 @@ template void Multigrid::vcycle() { - // The defect vector has been - // initialized by copy_to_mg. Now - // adjust the other vectors. + // The defect vector has been initialized by copy_to_mg. Now adjust the + // other vectors. solution.resize(minlevel, maxlevel); t.resize(minlevel, maxlevel); diff --git a/tests/multigrid/cycles.output b/tests/multigrid/cycles.output index 4a3ca5cc1d..fa9dab6bb5 100644 --- a/tests/multigrid/cycles.output +++ b/tests/multigrid/cycles.output @@ -47,242 +47,292 @@ DEAL::Smoothing on level 4 DEAL::Solution norm 0 DEAL::V-cycle leaving level 4 DEAL:: -DEAL::W-cycle entering level 4 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle entering level 4 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 4 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 4 -DEAL::W-cycle entering level 3 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 4 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 3 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 3 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 3 -DEAL::W-cycle entering level 2 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 3 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 2 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 2 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 2 -DEAL::W-cycle entering level 1 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 2 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 1 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 1 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 1 -DEAL::W-cycle entering level 0 -DEAL::W-cycle defect norm 0 -DEAL::W-cycle coarse level 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 1 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 0 +DEAL::W-cycle defect norm 0 +DEAL::W-cycle coarse level 0 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 1 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 1 -DEAL::W-cycle entering level 1 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle entering level 1 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 1 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 1 -DEAL::W-cycle entering level 0 -DEAL::W-cycle defect norm 0 -DEAL::W-cycle coarse level 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 1 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 0 +DEAL::W-cycle defect norm 0 +DEAL::W-cycle coarse level 0 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 1 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 1 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 2 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 2 -DEAL::W-cycle entering level 2 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle entering level 2 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 2 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 2 -DEAL::W-cycle entering level 1 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 2 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 1 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 1 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 1 -DEAL::W-cycle entering level 0 -DEAL::W-cycle defect norm 0 -DEAL::W-cycle coarse level 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 1 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 0 +DEAL::W-cycle defect norm 0 +DEAL::W-cycle coarse level 0 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 1 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 1 -DEAL::W-cycle entering level 1 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle entering level 1 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 1 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 1 -DEAL::W-cycle entering level 0 -DEAL::W-cycle defect norm 0 -DEAL::W-cycle coarse level 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 1 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 0 +DEAL::W-cycle defect norm 0 +DEAL::W-cycle coarse level 0 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 1 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 1 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 2 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 2 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 3 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 3 -DEAL::W-cycle entering level 3 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle entering level 3 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 3 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 3 -DEAL::W-cycle entering level 2 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 3 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 2 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 2 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 2 -DEAL::W-cycle entering level 1 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 2 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 1 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 1 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 1 -DEAL::W-cycle entering level 0 -DEAL::W-cycle defect norm 0 -DEAL::W-cycle coarse level 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 1 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 0 +DEAL::W-cycle defect norm 0 +DEAL::W-cycle coarse level 0 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 1 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 1 -DEAL::W-cycle entering level 1 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle entering level 1 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 1 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 1 -DEAL::W-cycle entering level 0 -DEAL::W-cycle defect norm 0 -DEAL::W-cycle coarse level 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 1 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 0 +DEAL::W-cycle defect norm 0 +DEAL::W-cycle coarse level 0 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 1 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 1 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 2 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 2 -DEAL::W-cycle entering level 2 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle entering level 2 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 2 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 2 -DEAL::W-cycle entering level 1 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 2 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 1 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 1 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 1 -DEAL::W-cycle entering level 0 -DEAL::W-cycle defect norm 0 -DEAL::W-cycle coarse level 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 1 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 0 +DEAL::W-cycle defect norm 0 +DEAL::W-cycle coarse level 0 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 1 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 1 -DEAL::W-cycle entering level 1 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle entering level 1 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 1 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 1 -DEAL::W-cycle entering level 0 -DEAL::W-cycle defect norm 0 -DEAL::W-cycle coarse level 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 1 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 0 +DEAL::W-cycle defect norm 0 +DEAL::W-cycle coarse level 0 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 1 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 1 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 2 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 2 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 3 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 3 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 4 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 4 DEAL:: -DEAL::F-cycle entering level 4 -DEAL::F-cycle defect norm 0 +DEAL::F-cycle entering level 4 +DEAL::F-cycle defect norm 0 DEAL::F-cycle smoothing level 4 -DEAL::F-cycle solution norm 0 -DEAL::F-cycle residual level 4 -DEAL::F-cycle entering level 3 -DEAL::F-cycle defect norm 0 +DEAL::F-cycle solution norm 0 +DEAL::F-cycle residual level 4 +DEAL::F-cycle residual norm 0 +DEAL::F-cycle entering level 3 +DEAL::F-cycle defect norm 0 DEAL::F-cycle smoothing level 3 -DEAL::F-cycle solution norm 0 -DEAL::F-cycle residual level 3 -DEAL::F-cycle entering level 2 -DEAL::F-cycle defect norm 0 +DEAL::F-cycle solution norm 0 +DEAL::F-cycle residual level 3 +DEAL::F-cycle residual norm 0 +DEAL::F-cycle entering level 2 +DEAL::F-cycle defect norm 0 DEAL::F-cycle smoothing level 2 -DEAL::F-cycle solution norm 0 -DEAL::F-cycle residual level 2 -DEAL::F-cycle entering level 1 -DEAL::F-cycle defect norm 0 +DEAL::F-cycle solution norm 0 +DEAL::F-cycle residual level 2 +DEAL::F-cycle residual norm 0 +DEAL::F-cycle entering level 1 +DEAL::F-cycle defect norm 0 DEAL::F-cycle smoothing level 1 -DEAL::F-cycle solution norm 0 -DEAL::F-cycle residual level 1 -DEAL::F-cycle entering level 0 -DEAL::F-cycle defect norm 0 -DEAL::F-cycle coarse level 0 +DEAL::F-cycle solution norm 0 +DEAL::F-cycle residual level 1 +DEAL::F-cycle residual norm 0 +DEAL::F-cycle entering level 0 +DEAL::F-cycle defect norm 0 +DEAL::F-cycle coarse level 0 DEAL::F-cycle Defect norm 0 DEAL::F-cycle smoothing level 1 +DEAL::F-cycle solution norm 0 DEAL::F-cycle leaving level 1 -DEAL::V-cycle entering level 1 -DEAL::V-cycle defect norm 0 +DEAL::V-cycle entering level 1 +DEAL::V-cycle defect norm 0 DEAL::V-cycle smoothing level 1 -DEAL::V-cycle solution norm 0 -DEAL::V-cycle residual level 1 -DEAL::V-cycle entering level 0 -DEAL::V-cycle defect norm 0 -DEAL::V-cycle coarse level 0 +DEAL::V-cycle solution norm 0 +DEAL::V-cycle residual level 1 +DEAL::V-cycle residual norm 0 +DEAL::V-cycle entering level 0 +DEAL::V-cycle defect norm 0 +DEAL::V-cycle coarse level 0 DEAL::V-cycle Defect norm 0 DEAL::V-cycle smoothing level 1 +DEAL::V-cycle solution norm 0 DEAL::V-cycle leaving level 1 DEAL::F-cycle Defect norm 0 DEAL::F-cycle smoothing level 2 +DEAL::F-cycle solution norm 0 DEAL::F-cycle leaving level 2 -DEAL::V-cycle entering level 2 -DEAL::V-cycle defect norm 0 +DEAL::V-cycle entering level 2 +DEAL::V-cycle defect norm 0 DEAL::V-cycle smoothing level 2 -DEAL::V-cycle solution norm 0 -DEAL::V-cycle residual level 2 -DEAL::V-cycle entering level 1 -DEAL::V-cycle defect norm 0 +DEAL::V-cycle solution norm 0 +DEAL::V-cycle residual level 2 +DEAL::V-cycle residual norm 0 +DEAL::V-cycle entering level 1 +DEAL::V-cycle defect norm 0 DEAL::V-cycle smoothing level 1 -DEAL::V-cycle solution norm 0 -DEAL::V-cycle residual level 1 -DEAL::V-cycle entering level 0 -DEAL::V-cycle defect norm 0 -DEAL::V-cycle coarse level 0 +DEAL::V-cycle solution norm 0 +DEAL::V-cycle residual level 1 +DEAL::V-cycle residual norm 0 +DEAL::V-cycle entering level 0 +DEAL::V-cycle defect norm 0 +DEAL::V-cycle coarse level 0 DEAL::V-cycle Defect norm 0 DEAL::V-cycle smoothing level 1 +DEAL::V-cycle solution norm 0 DEAL::V-cycle leaving level 1 DEAL::V-cycle Defect norm 0 DEAL::V-cycle smoothing level 2 +DEAL::V-cycle solution norm 0 DEAL::V-cycle leaving level 2 DEAL::F-cycle Defect norm 0 DEAL::F-cycle smoothing level 3 +DEAL::F-cycle solution norm 0 DEAL::F-cycle leaving level 3 -DEAL::V-cycle entering level 3 -DEAL::V-cycle defect norm 0 +DEAL::V-cycle entering level 3 +DEAL::V-cycle defect norm 0 DEAL::V-cycle smoothing level 3 -DEAL::V-cycle solution norm 0 -DEAL::V-cycle residual level 3 -DEAL::V-cycle entering level 2 -DEAL::V-cycle defect norm 0 +DEAL::V-cycle solution norm 0 +DEAL::V-cycle residual level 3 +DEAL::V-cycle residual norm 0 +DEAL::V-cycle entering level 2 +DEAL::V-cycle defect norm 0 DEAL::V-cycle smoothing level 2 -DEAL::V-cycle solution norm 0 -DEAL::V-cycle residual level 2 -DEAL::V-cycle entering level 1 -DEAL::V-cycle defect norm 0 +DEAL::V-cycle solution norm 0 +DEAL::V-cycle residual level 2 +DEAL::V-cycle residual norm 0 +DEAL::V-cycle entering level 1 +DEAL::V-cycle defect norm 0 DEAL::V-cycle smoothing level 1 -DEAL::V-cycle solution norm 0 -DEAL::V-cycle residual level 1 -DEAL::V-cycle entering level 0 -DEAL::V-cycle defect norm 0 -DEAL::V-cycle coarse level 0 +DEAL::V-cycle solution norm 0 +DEAL::V-cycle residual level 1 +DEAL::V-cycle residual norm 0 +DEAL::V-cycle entering level 0 +DEAL::V-cycle defect norm 0 +DEAL::V-cycle coarse level 0 DEAL::V-cycle Defect norm 0 DEAL::V-cycle smoothing level 1 +DEAL::V-cycle solution norm 0 DEAL::V-cycle leaving level 1 DEAL::V-cycle Defect norm 0 DEAL::V-cycle smoothing level 2 +DEAL::V-cycle solution norm 0 DEAL::V-cycle leaving level 2 DEAL::V-cycle Defect norm 0 DEAL::V-cycle smoothing level 3 +DEAL::V-cycle solution norm 0 DEAL::V-cycle leaving level 3 DEAL::F-cycle Defect norm 0 DEAL::F-cycle smoothing level 4 +DEAL::F-cycle solution norm 0 DEAL::F-cycle leaving level 4 DEAL::V-cycle entering level 5 DEAL::V-cycle Defect norm 0 @@ -321,129 +371,155 @@ DEAL::Smoothing on level 5 DEAL::Solution norm 0 DEAL::V-cycle leaving level 5 DEAL:: -DEAL::W-cycle entering level 5 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle entering level 5 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 5 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 5 -DEAL::W-cycle entering level 4 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 5 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 4 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 4 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 4 -DEAL::W-cycle entering level 3 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 4 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 3 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 3 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 3 -DEAL::W-cycle entering level 2 -DEAL::W-cycle defect norm 0 -DEAL::W-cycle coarse level 2 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 3 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 2 +DEAL::W-cycle defect norm 0 +DEAL::W-cycle coarse level 2 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 3 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 3 -DEAL::W-cycle entering level 3 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle entering level 3 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 3 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 3 -DEAL::W-cycle entering level 2 -DEAL::W-cycle defect norm 0 -DEAL::W-cycle coarse level 2 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 3 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 2 +DEAL::W-cycle defect norm 0 +DEAL::W-cycle coarse level 2 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 3 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 3 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 4 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 4 -DEAL::W-cycle entering level 4 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle entering level 4 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 4 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 4 -DEAL::W-cycle entering level 3 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 4 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 3 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 3 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 3 -DEAL::W-cycle entering level 2 -DEAL::W-cycle defect norm 0 -DEAL::W-cycle coarse level 2 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 3 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 2 +DEAL::W-cycle defect norm 0 +DEAL::W-cycle coarse level 2 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 3 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 3 -DEAL::W-cycle entering level 3 -DEAL::W-cycle defect norm 0 +DEAL::W-cycle entering level 3 +DEAL::W-cycle defect norm 0 DEAL::W-cycle smoothing level 3 -DEAL::W-cycle solution norm 0 -DEAL::W-cycle residual level 3 -DEAL::W-cycle entering level 2 -DEAL::W-cycle defect norm 0 -DEAL::W-cycle coarse level 2 +DEAL::W-cycle solution norm 0 +DEAL::W-cycle residual level 3 +DEAL::W-cycle residual norm 0 +DEAL::W-cycle entering level 2 +DEAL::W-cycle defect norm 0 +DEAL::W-cycle coarse level 2 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 3 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 3 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 4 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 4 DEAL::W-cycle Defect norm 0 DEAL::W-cycle smoothing level 5 +DEAL::W-cycle solution norm 0 DEAL::W-cycle leaving level 5 DEAL:: -DEAL::F-cycle entering level 5 -DEAL::F-cycle defect norm 0 +DEAL::F-cycle entering level 5 +DEAL::F-cycle defect norm 0 DEAL::F-cycle smoothing level 5 -DEAL::F-cycle solution norm 0 -DEAL::F-cycle residual level 5 -DEAL::F-cycle entering level 4 -DEAL::F-cycle defect norm 0 +DEAL::F-cycle solution norm 0 +DEAL::F-cycle residual level 5 +DEAL::F-cycle residual norm 0 +DEAL::F-cycle entering level 4 +DEAL::F-cycle defect norm 0 DEAL::F-cycle smoothing level 4 -DEAL::F-cycle solution norm 0 -DEAL::F-cycle residual level 4 -DEAL::F-cycle entering level 3 -DEAL::F-cycle defect norm 0 +DEAL::F-cycle solution norm 0 +DEAL::F-cycle residual level 4 +DEAL::F-cycle residual norm 0 +DEAL::F-cycle entering level 3 +DEAL::F-cycle defect norm 0 DEAL::F-cycle smoothing level 3 -DEAL::F-cycle solution norm 0 -DEAL::F-cycle residual level 3 -DEAL::F-cycle entering level 2 -DEAL::F-cycle defect norm 0 -DEAL::F-cycle coarse level 2 +DEAL::F-cycle solution norm 0 +DEAL::F-cycle residual level 3 +DEAL::F-cycle residual norm 0 +DEAL::F-cycle entering level 2 +DEAL::F-cycle defect norm 0 +DEAL::F-cycle coarse level 2 DEAL::F-cycle Defect norm 0 DEAL::F-cycle smoothing level 3 +DEAL::F-cycle solution norm 0 DEAL::F-cycle leaving level 3 -DEAL::V-cycle entering level 3 -DEAL::V-cycle defect norm 0 +DEAL::V-cycle entering level 3 +DEAL::V-cycle defect norm 0 DEAL::V-cycle smoothing level 3 -DEAL::V-cycle solution norm 0 -DEAL::V-cycle residual level 3 -DEAL::V-cycle entering level 2 -DEAL::V-cycle defect norm 0 -DEAL::V-cycle coarse level 2 +DEAL::V-cycle solution norm 0 +DEAL::V-cycle residual level 3 +DEAL::V-cycle residual norm 0 +DEAL::V-cycle entering level 2 +DEAL::V-cycle defect norm 0 +DEAL::V-cycle coarse level 2 DEAL::V-cycle Defect norm 0 DEAL::V-cycle smoothing level 3 +DEAL::V-cycle solution norm 0 DEAL::V-cycle leaving level 3 DEAL::F-cycle Defect norm 0 DEAL::F-cycle smoothing level 4 +DEAL::F-cycle solution norm 0 DEAL::F-cycle leaving level 4 -DEAL::V-cycle entering level 4 -DEAL::V-cycle defect norm 0 +DEAL::V-cycle entering level 4 +DEAL::V-cycle defect norm 0 DEAL::V-cycle smoothing level 4 -DEAL::V-cycle solution norm 0 -DEAL::V-cycle residual level 4 -DEAL::V-cycle entering level 3 -DEAL::V-cycle defect norm 0 +DEAL::V-cycle solution norm 0 +DEAL::V-cycle residual level 4 +DEAL::V-cycle residual norm 0 +DEAL::V-cycle entering level 3 +DEAL::V-cycle defect norm 0 DEAL::V-cycle smoothing level 3 -DEAL::V-cycle solution norm 0 -DEAL::V-cycle residual level 3 -DEAL::V-cycle entering level 2 -DEAL::V-cycle defect norm 0 -DEAL::V-cycle coarse level 2 +DEAL::V-cycle solution norm 0 +DEAL::V-cycle residual level 3 +DEAL::V-cycle residual norm 0 +DEAL::V-cycle entering level 2 +DEAL::V-cycle defect norm 0 +DEAL::V-cycle coarse level 2 DEAL::V-cycle Defect norm 0 DEAL::V-cycle smoothing level 3 +DEAL::V-cycle solution norm 0 DEAL::V-cycle leaving level 3 DEAL::V-cycle Defect norm 0 DEAL::V-cycle smoothing level 4 +DEAL::V-cycle solution norm 0 DEAL::V-cycle leaving level 4 DEAL::F-cycle Defect norm 0 DEAL::F-cycle smoothing level 5 +DEAL::F-cycle solution norm 0 DEAL::F-cycle leaving level 5 diff --git a/tests/multigrid/step-16-05.cc b/tests/multigrid/step-16-05.cc new file mode 100644 index 0000000000..d0b5bdda23 --- /dev/null +++ b/tests/multigrid/step-16-05.cc @@ -0,0 +1,509 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2003 - 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Similar to step-16-04 but using a W cycle instead of a V-cycle and no +// adaptive meshes + +#include "../tests.h" +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +using namespace dealii; + +template +class LaplaceProblem +{ +public: + LaplaceProblem (const unsigned int deg); + void run (); + +private: + void setup_system (); + void assemble_system (); + void assemble_multigrid (); + void solve (); + + Triangulation triangulation; + FE_Q fe; + DoFHandler mg_dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + ConstraintMatrix hanging_node_constraints; + ConstraintMatrix constraints; + + LinearAlgebra::distributed::Vector solution; + LinearAlgebra::distributed::Vector system_rhs; + + const unsigned int degree; + const unsigned int min_level; + + MGLevelObject mg_sparsity_patterns; + MGLevelObject > mg_matrices; + MGLevelObject > mg_interface_matrices; + MGConstrainedDoFs mg_constrained_dofs; +}; + + +template +class Coefficient : public Function +{ +public: + Coefficient () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component = 0) const; +}; + + + +template +double Coefficient::value (const Point &p, + const unsigned int) const +{ + if (p.square() < 0.5*0.5) + return 20; + else + return 1; +} + + + +template +void Coefficient::value_list (const std::vector > &points, + std::vector &values, + const unsigned int component) const +{ + const unsigned int n_points = points.size(); + + Assert (values.size() == n_points, + ExcDimensionMismatch (values.size(), n_points)); + + Assert (component == 0, + ExcIndexRange (component, 0, 1)); + + for (unsigned int i=0; i::value (points[i]); +} + + + +template +LaplaceProblem::LaplaceProblem (const unsigned int degree) + : + triangulation (Triangulation:: + limit_level_difference_at_vertices), + fe (degree), + mg_dof_handler (triangulation), + degree(degree), + min_level(2) +{} + + + +template +void LaplaceProblem::setup_system () +{ + mg_dof_handler.distribute_dofs(fe); + mg_dof_handler.distribute_mg_dofs (fe); + + sparsity_pattern.reinit (mg_dof_handler.n_dofs(), + mg_dof_handler.n_dofs(), + mg_dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (mg_dof_handler, sparsity_pattern); + + solution.reinit (mg_dof_handler.n_dofs()); + system_rhs.reinit (mg_dof_handler.n_dofs()); + + constraints.clear (); + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints); + DoFTools::make_hanging_node_constraints (mg_dof_handler, hanging_node_constraints); + typename FunctionMap::type dirichlet_boundary; + ZeroFunction homogeneous_dirichlet_bc (1); + dirichlet_boundary[0] = &homogeneous_dirichlet_bc; + MappingQGeneric mapping(1); + VectorTools::interpolate_boundary_values (mapping, + mg_dof_handler, + dirichlet_boundary, + constraints); + constraints.close (); + hanging_node_constraints.close (); + constraints.condense (sparsity_pattern); + sparsity_pattern.compress(); + system_matrix.reinit (sparsity_pattern); + + mg_constrained_dofs.clear(); + mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary); + const unsigned int n_levels = triangulation.n_levels(); + + mg_interface_matrices.resize(min_level, n_levels-1); + mg_interface_matrices.clear_elements (); + mg_matrices.resize(min_level, n_levels-1); + mg_matrices.clear_elements (); + mg_sparsity_patterns.resize(min_level, n_levels-1); + + for (unsigned int level=min_level; level +void LaplaceProblem::assemble_system () +{ + Vector tmp(system_rhs.size()); + const QGauss quadrature_formula(degree+1); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + typename DoFHandler::active_cell_iterator + cell = mg_dof_handler.begin_active(), + endc = mg_dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit (cell); + + coefficient.value_list (fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, cell_rhs, + local_dof_indices, + system_matrix, tmp); + } + for (unsigned int i=0; i +void LaplaceProblem::assemble_multigrid () +{ + QGauss quadrature_formula(1+degree); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + std::vector boundary_constraints (triangulation.n_levels()); + ConstraintMatrix empty_constraints; + for (unsigned int level=min_level; level::cell_iterator cell = mg_dof_handler.begin(min_level), + endc = mg_dof_handler.end(); + + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + fe_values.reinit (cell); + + coefficient.value_list (fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_point=0; q_pointget_mg_dof_indices (local_dof_indices); + + boundary_constraints[cell->level()] + .distribute_local_to_global (cell_matrix, + local_dof_indices, + mg_matrices[cell->level()]); + + const unsigned int lvl = cell->level(); + + for (unsigned int i=0; ilevel()]); + } +} + + + +template +void LaplaceProblem::solve () +{ + MGTransferPrebuilt > + mg_transfer(hanging_node_constraints, mg_constrained_dofs); + mg_transfer.build_matrices(mg_dof_handler); + + SolverControl coarse_solver_control (1000, 1e-10, false, false); + SolverCG > coarse_solver(coarse_solver_control); + PreconditionIdentity id; + MGCoarseGridLACIteration >,LinearAlgebra::distributed::Vector > + coarse_grid_solver(coarse_solver, mg_matrices[min_level], id); + deallog << " Size of coarse grid matrix: " << mg_matrices[min_level].m() << std::endl; + + typedef PreconditionChebyshev,LinearAlgebra::distributed::Vector > Smoother; + GrowingVectorMemory > vector_memory; + MGSmootherPrecondition, Smoother, LinearAlgebra::distributed::Vector > + mg_smoother; + typename Smoother::AdditionalData smoother_data; + smoother_data.smoothing_range = 20.; + smoother_data.degree = 2; + smoother_data.eig_cg_n_iterations = 20; + mg_smoother.initialize(mg_matrices, smoother_data); + + mg::Matrix > mg_matrix(mg_matrices); + mg::Matrix > mg_interface_up(mg_interface_matrices); + mg::Matrix > mg_interface_down(mg_interface_matrices); + + Multigrid > mg(mg_matrix, + coarse_grid_solver, + mg_transfer, + mg_smoother, + mg_smoother, + min_level, + triangulation.n_global_levels()-1, + Multigrid >::w_cycle); + mg.set_edge_matrices(mg_interface_down, mg_interface_up); + + PreconditionMG, MGTransferPrebuilt > > + preconditioner(mg_dof_handler, mg, mg_transfer); + + SolverControl solver_control (1000, 1e-12); + SolverCG > cg (solver_control); + + solution = 0; + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + constraints.distribute (solution); + + deallog << " " << solver_control.last_step() + << " CG iterations needed to obtain convergence." + << std::endl; +} + + + +template +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<4; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (min_level+1); + } + else + triangulation.refine_global(1); + + + deallog << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system (); + + deallog << " Number of degrees of freedom: " + << mg_dof_handler.n_dofs() + << " (by level: "; + for (unsigned int level=min_level; level laplace_problem(1); + laplace_problem.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/multigrid/step-16-05.output b/tests/multigrid/step-16-05.output new file mode 100644 index 0000000000..b7f18bfa52 --- /dev/null +++ b/tests/multigrid/step-16-05.output @@ -0,0 +1,57 @@ + +DEAL::Cycle 0: +DEAL:: Number of active cells: 64 +DEAL:: Number of degrees of freedom: 81 (by level: 25, 81) +DEAL:: Size of coarse grid matrix: 25 +DEAL:cg::Starting value 0.9798 +DEAL:cg::Convergence step 7 value 0 +DEAL:cg::Starting value 0.9938 +DEAL:cg::Convergence step 18 value 0 +DEAL:cg::Starting value 0.1094 +DEAL:cg::Convergence step 10 value 0 +DEAL:: 10 CG iterations needed to obtain convergence. +DEAL::Cycle 1: +DEAL:: Number of active cells: 256 +DEAL:: Number of degrees of freedom: 289 (by level: 25, 81, 289) +DEAL:: Size of coarse grid matrix: 25 +DEAL:cg::Starting value 0.9798 +DEAL:cg::Convergence step 7 value 0 +DEAL:cg::Starting value 0.9938 +DEAL:cg::Convergence step 18 value 0 +DEAL:cg::Starting value 0.9983 +DEAL:cg::Failure step 20 value 0.001167 +DEAL:cg::Starting value 0.05859 +DEAL:cg::Convergence step 10 value 0 +DEAL:: 10 CG iterations needed to obtain convergence. +DEAL::Cycle 2: +DEAL:: Number of active cells: 1024 +DEAL:: Number of degrees of freedom: 1089 (by level: 25, 81, 289, 1089) +DEAL:: Size of coarse grid matrix: 25 +DEAL:cg::Starting value 0.9798 +DEAL:cg::Convergence step 7 value 0 +DEAL:cg::Starting value 0.9938 +DEAL:cg::Convergence step 18 value 0 +DEAL:cg::Starting value 0.9983 +DEAL:cg::Failure step 20 value 0.001167 +DEAL:cg::Starting value 0.9995 +DEAL:cg::Failure step 20 value 0.1128 +DEAL:cg::Starting value 0.03027 +DEAL:cg::Convergence step 11 value 0 +DEAL:: 11 CG iterations needed to obtain convergence. +DEAL::Cycle 3: +DEAL:: Number of active cells: 4096 +DEAL:: Number of degrees of freedom: 4225 (by level: 25, 81, 289, 1089, 4225) +DEAL:: Size of coarse grid matrix: 25 +DEAL:cg::Starting value 0.9798 +DEAL:cg::Convergence step 7 value 0 +DEAL:cg::Starting value 0.9938 +DEAL:cg::Convergence step 18 value 0 +DEAL:cg::Starting value 0.9983 +DEAL:cg::Failure step 20 value 0.001167 +DEAL:cg::Starting value 0.9995 +DEAL:cg::Failure step 20 value 0.1128 +DEAL:cg::Starting value 0.9999 +DEAL:cg::Failure step 20 value 1.273 +DEAL:cg::Starting value 0.01538 +DEAL:cg::Convergence step 11 value 0 +DEAL:: 11 CG iterations needed to obtain convergence. diff --git a/tests/multigrid/step-16-06.cc b/tests/multigrid/step-16-06.cc new file mode 100644 index 0000000000..878058162c --- /dev/null +++ b/tests/multigrid/step-16-06.cc @@ -0,0 +1,509 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2003 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Similar to step-16-04 but using a multigrid F cycle and without adaptive +// grid + +#include "../tests.h" +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +using namespace dealii; + +template +class LaplaceProblem +{ +public: + LaplaceProblem (const unsigned int deg); + void run (); + +private: + void setup_system (); + void assemble_system (); + void assemble_multigrid (); + void solve (); + + Triangulation triangulation; + FE_Q fe; + DoFHandler mg_dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + ConstraintMatrix hanging_node_constraints; + ConstraintMatrix constraints; + + LinearAlgebra::distributed::Vector solution; + LinearAlgebra::distributed::Vector system_rhs; + + const unsigned int degree; + const unsigned int min_level; + + MGLevelObject mg_sparsity_patterns; + MGLevelObject > mg_matrices; + MGLevelObject > mg_interface_matrices; + MGConstrainedDoFs mg_constrained_dofs; +}; + + +template +class Coefficient : public Function +{ +public: + Coefficient () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component = 0) const; +}; + + + +template +double Coefficient::value (const Point &p, + const unsigned int) const +{ + if (p.square() < 0.5*0.5) + return 20; + else + return 1; +} + + + +template +void Coefficient::value_list (const std::vector > &points, + std::vector &values, + const unsigned int component) const +{ + const unsigned int n_points = points.size(); + + Assert (values.size() == n_points, + ExcDimensionMismatch (values.size(), n_points)); + + Assert (component == 0, + ExcIndexRange (component, 0, 1)); + + for (unsigned int i=0; i::value (points[i]); +} + + + +template +LaplaceProblem::LaplaceProblem (const unsigned int degree) + : + triangulation (Triangulation:: + limit_level_difference_at_vertices), + fe (degree), + mg_dof_handler (triangulation), + degree(degree), + min_level(2) +{} + + + +template +void LaplaceProblem::setup_system () +{ + mg_dof_handler.distribute_dofs(fe); + mg_dof_handler.distribute_mg_dofs (fe); + + sparsity_pattern.reinit (mg_dof_handler.n_dofs(), + mg_dof_handler.n_dofs(), + mg_dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (mg_dof_handler, sparsity_pattern); + + solution.reinit (mg_dof_handler.n_dofs()); + system_rhs.reinit (mg_dof_handler.n_dofs()); + + constraints.clear (); + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints); + DoFTools::make_hanging_node_constraints (mg_dof_handler, hanging_node_constraints); + typename FunctionMap::type dirichlet_boundary; + ZeroFunction homogeneous_dirichlet_bc (1); + dirichlet_boundary[0] = &homogeneous_dirichlet_bc; + MappingQGeneric mapping(1); + VectorTools::interpolate_boundary_values (mapping, + mg_dof_handler, + dirichlet_boundary, + constraints); + constraints.close (); + hanging_node_constraints.close (); + constraints.condense (sparsity_pattern); + sparsity_pattern.compress(); + system_matrix.reinit (sparsity_pattern); + + mg_constrained_dofs.clear(); + mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary); + const unsigned int n_levels = triangulation.n_levels(); + + mg_interface_matrices.resize(min_level, n_levels-1); + mg_interface_matrices.clear_elements (); + mg_matrices.resize(min_level, n_levels-1); + mg_matrices.clear_elements (); + mg_sparsity_patterns.resize(min_level, n_levels-1); + + for (unsigned int level=min_level; level +void LaplaceProblem::assemble_system () +{ + Vector tmp(system_rhs.size()); + const QGauss quadrature_formula(degree+1); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + typename DoFHandler::active_cell_iterator + cell = mg_dof_handler.begin_active(), + endc = mg_dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit (cell); + + coefficient.value_list (fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, cell_rhs, + local_dof_indices, + system_matrix, tmp); + } + for (unsigned int i=0; i +void LaplaceProblem::assemble_multigrid () +{ + QGauss quadrature_formula(1+degree); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + std::vector boundary_constraints (triangulation.n_levels()); + ConstraintMatrix empty_constraints; + for (unsigned int level=min_level; level::cell_iterator cell = mg_dof_handler.begin(min_level), + endc = mg_dof_handler.end(); + + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + fe_values.reinit (cell); + + coefficient.value_list (fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_point=0; q_pointget_mg_dof_indices (local_dof_indices); + + boundary_constraints[cell->level()] + .distribute_local_to_global (cell_matrix, + local_dof_indices, + mg_matrices[cell->level()]); + + const unsigned int lvl = cell->level(); + + for (unsigned int i=0; ilevel()]); + } +} + + + +template +void LaplaceProblem::solve () +{ + MGTransferPrebuilt > + mg_transfer(hanging_node_constraints, mg_constrained_dofs); + mg_transfer.build_matrices(mg_dof_handler); + + SolverControl coarse_solver_control (1000, 1e-10, false, false); + SolverCG > coarse_solver(coarse_solver_control); + PreconditionIdentity id; + MGCoarseGridLACIteration >,LinearAlgebra::distributed::Vector > + coarse_grid_solver(coarse_solver, mg_matrices[min_level], id); + deallog << " Size of coarse grid matrix: " << mg_matrices[min_level].m() << std::endl; + + typedef PreconditionChebyshev,LinearAlgebra::distributed::Vector > Smoother; + GrowingVectorMemory > vector_memory; + MGSmootherPrecondition, Smoother, LinearAlgebra::distributed::Vector > + mg_smoother; + typename Smoother::AdditionalData smoother_data; + smoother_data.smoothing_range = 20.; + smoother_data.degree = 2; + smoother_data.eig_cg_n_iterations = 20; + mg_smoother.initialize(mg_matrices, smoother_data); + + mg::Matrix > mg_matrix(mg_matrices); + mg::Matrix > mg_interface_up(mg_interface_matrices); + mg::Matrix > mg_interface_down(mg_interface_matrices); + + Multigrid > mg(mg_matrix, + coarse_grid_solver, + mg_transfer, + mg_smoother, + mg_smoother, + min_level, + triangulation.n_global_levels()-1, + Multigrid >::f_cycle); + mg.set_edge_matrices(mg_interface_down, mg_interface_up); + + PreconditionMG, MGTransferPrebuilt > > + preconditioner(mg_dof_handler, mg, mg_transfer); + + SolverControl solver_control (1000, 1e-12); + SolverCG > cg (solver_control); + + solution = 0; + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + constraints.distribute (solution); + + deallog << " " << solver_control.last_step() + << " CG iterations needed to obtain convergence." + << std::endl; +} + + + +template +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<4; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (min_level+1); + } + else + triangulation.refine_global (1); + + + deallog << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system (); + + deallog << " Number of degrees of freedom: " + << mg_dof_handler.n_dofs() + << " (by level: "; + for (unsigned int level=min_level; level laplace_problem(1); + laplace_problem.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/multigrid/step-16-06.output b/tests/multigrid/step-16-06.output new file mode 100644 index 0000000000..b7f18bfa52 --- /dev/null +++ b/tests/multigrid/step-16-06.output @@ -0,0 +1,57 @@ + +DEAL::Cycle 0: +DEAL:: Number of active cells: 64 +DEAL:: Number of degrees of freedom: 81 (by level: 25, 81) +DEAL:: Size of coarse grid matrix: 25 +DEAL:cg::Starting value 0.9798 +DEAL:cg::Convergence step 7 value 0 +DEAL:cg::Starting value 0.9938 +DEAL:cg::Convergence step 18 value 0 +DEAL:cg::Starting value 0.1094 +DEAL:cg::Convergence step 10 value 0 +DEAL:: 10 CG iterations needed to obtain convergence. +DEAL::Cycle 1: +DEAL:: Number of active cells: 256 +DEAL:: Number of degrees of freedom: 289 (by level: 25, 81, 289) +DEAL:: Size of coarse grid matrix: 25 +DEAL:cg::Starting value 0.9798 +DEAL:cg::Convergence step 7 value 0 +DEAL:cg::Starting value 0.9938 +DEAL:cg::Convergence step 18 value 0 +DEAL:cg::Starting value 0.9983 +DEAL:cg::Failure step 20 value 0.001167 +DEAL:cg::Starting value 0.05859 +DEAL:cg::Convergence step 10 value 0 +DEAL:: 10 CG iterations needed to obtain convergence. +DEAL::Cycle 2: +DEAL:: Number of active cells: 1024 +DEAL:: Number of degrees of freedom: 1089 (by level: 25, 81, 289, 1089) +DEAL:: Size of coarse grid matrix: 25 +DEAL:cg::Starting value 0.9798 +DEAL:cg::Convergence step 7 value 0 +DEAL:cg::Starting value 0.9938 +DEAL:cg::Convergence step 18 value 0 +DEAL:cg::Starting value 0.9983 +DEAL:cg::Failure step 20 value 0.001167 +DEAL:cg::Starting value 0.9995 +DEAL:cg::Failure step 20 value 0.1128 +DEAL:cg::Starting value 0.03027 +DEAL:cg::Convergence step 11 value 0 +DEAL:: 11 CG iterations needed to obtain convergence. +DEAL::Cycle 3: +DEAL:: Number of active cells: 4096 +DEAL:: Number of degrees of freedom: 4225 (by level: 25, 81, 289, 1089, 4225) +DEAL:: Size of coarse grid matrix: 25 +DEAL:cg::Starting value 0.9798 +DEAL:cg::Convergence step 7 value 0 +DEAL:cg::Starting value 0.9938 +DEAL:cg::Convergence step 18 value 0 +DEAL:cg::Starting value 0.9983 +DEAL:cg::Failure step 20 value 0.001167 +DEAL:cg::Starting value 0.9995 +DEAL:cg::Failure step 20 value 0.1128 +DEAL:cg::Starting value 0.9999 +DEAL:cg::Failure step 20 value 1.273 +DEAL:cg::Starting value 0.01538 +DEAL:cg::Convergence step 11 value 0 +DEAL:: 11 CG iterations needed to obtain convergence. -- 2.39.5