From: kanschat Date: Thu, 17 Jun 2010 04:01:18 +0000 (+0000) Subject: documentation X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f3a507d72ef5880863610eae70679a98c9060460;p=dealii-svn.git documentation git-svn-id: https://svn.dealii.org/trunk@21220 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/multigrid/multigrid.h b/deal.II/deal.II/include/multigrid/multigrid.h index 9ffeddc4b0..f48b57a0e1 100644 --- a/deal.II/deal.II/include/multigrid/multigrid.h +++ b/deal.II/deal.II/include/multigrid/multigrid.h @@ -254,7 +254,7 @@ class Multigrid : public Subscriptor * matrices originate from * discontinuous Galerkin methods * (see FE_DGQ etc.), where they - * correspond tu the edge fluxes + * correspond to the edge fluxes * at the refinement edge between * two levels. * diff --git a/deal.II/examples/step-39/doc/results.dox b/deal.II/examples/step-39/doc/results.dox index 91c515a668..20b3dc462d 100644 --- a/deal.II/examples/step-39/doc/results.dox +++ b/deal.II/examples/step-39/doc/results.dox @@ -14,7 +14,7 @@ DEAL::Assemble right hand side DEAL::Solve DEAL:cg::Starting value 27.1275 DEAL:cg::Convergence step 1 value 1.97998e-14 -DEAL::Error 0.161536 +DEAL::Error 0.161172 DEAL::Estimate 1.35839 DEAL::Writing solution to ... DEAL:: @@ -26,8 +26,8 @@ DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve DEAL:cg::Starting value 35.5356 -DEAL:cg::Convergence step 10 value 1.24790e-13 -DEAL::Error 0.164997 +DEAL:cg::Convergence step 14 value 3.21479e-13 +DEAL::Error 0.164760 DEAL::Estimate 1.08528 DEAL::Writing solution to ... DEAL:: @@ -39,8 +39,8 @@ DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve DEAL:cg::Starting value 37.0552 -DEAL:cg::Convergence step 10 value 4.74024e-14 -DEAL::Error 0.113596 +DEAL:cg::Convergence step 14 value 6.05416e-13 +DEAL::Error 0.113503 DEAL::Estimate 0.990460 DEAL::Writing solution to ... @@ -54,8 +54,8 @@ DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve DEAL:cg::Starting value 38.5798 -DEAL:cg::Convergence step 12 value 7.56337e-14 -DEAL::Error 0.0101279 +DEAL:cg::Convergence step 17 value 2.64999e-13 +DEAL::Error 0.0101278 DEAL::Estimate 0.0957571 DEAL::Writing solution to ... DEAL:: @@ -67,21 +67,22 @@ DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve DEAL:cg::Starting value 44.1721 -DEAL:cg::Convergence step 12 value 9.49720e-14 -DEAL::Error 0.00716966 +DEAL:cg::Convergence step 17 value 3.18657e-13 +DEAL::Error 0.00716962 DEAL::Estimate 0.0681646 DEAL::Writing solution to ... +DEAL:: @endcode This log for instance shows that the number of conjugate gradient -iteration steps is constant at approximately 12. +iteration steps is constant at approximately 17.

Postprocessing of the logfile

@image html "step-39-convergence.png" Using the perl script postprocess.pl, we extract relevant data into output.dat, which can be used to plot graphs with -gnuplot. The graph on the right for instance was produced with +gnuplot. The graph above for instance was produced with @code set style data linespoints diff --git a/deal.II/examples/step-39/step-39.cc b/deal.II/examples/step-39/step-39.cc index 8563ea4101..2028ed8672 100644 --- a/deal.II/examples/step-39/step-39.cc +++ b/deal.II/examples/step-39/step-39.cc @@ -88,9 +88,16 @@ Functions::SlitSingularityFunction<2> exact_solution; // the information needed for the // local integration is provided by // MeshWorker::IntegrationInfo. Note - // that this public interface cannot + // that the signature of the functions cannot // be changed, because it is expected // by MeshWorker::integration_loop(). + + // The first class defining local + // integrators is responsible for + // computing cell and face + // matrices. It is used to assemble + // the global matrix as well as the + // level matrices. template class MatrixIntegrator : public Subscriptor { @@ -150,7 +157,8 @@ void MatrixIntegrator::bdry( * fe.JxW(k); } - + // Interior faces use the interior + // penalty method template void MatrixIntegrator::face( MeshWorker::DoFInfo& dinfo1, @@ -200,7 +208,12 @@ const double penalty = penalty1 + penalty2; } } - + // The second local integrator builds + // the right hand side. In our + // example, the right hand side + // function is zero, such that only + // the boundary condition is set here + // in weak form. template class RHSIntegrator : public Subscriptor { @@ -247,6 +260,11 @@ void RHSIntegrator::face(MeshWorker::DoFInfo&, {} + // The third local integrator is + // responsible for the contributions + // to the error estimate. This is the + // standard energy estimator due to + // Karakashian and Pascal (2003). template class Estimator : public Subscriptor { @@ -260,6 +278,10 @@ class Estimator : public Subscriptor }; + // The cell contribution is the + // Laplacian of the discrete + // solution, since the right hand + // side is zero. template void Estimator::cell(MeshWorker::DoFInfo& dinfo, typename MeshWorker::IntegrationInfo& info) { @@ -274,7 +296,12 @@ void Estimator::cell(MeshWorker::DoFInfo& dinfo, typename MeshWorker:: dinfo.value(0) = std::sqrt(dinfo.value(0)); } - + // At the boundary, we use simply a + // weighted form of the boundary + // residual, namely the norm of the + // difference between the finite + // element solution and the correct + // boundary condition. template void Estimator::bdry(MeshWorker::DoFInfo& dinfo, typename MeshWorker::IntegrationInfo& info) { @@ -295,6 +322,10 @@ void Estimator::bdry(MeshWorker::DoFInfo& dinfo, typename MeshWorker:: } + // Finally, on interior faces, the + // estimator consists of the jumps of + // the solution and its normal + // derivative, weighted appropriately. template void Estimator::face(MeshWorker::DoFInfo& dinfo1, MeshWorker::DoFInfo& dinfo2, @@ -326,11 +357,17 @@ void Estimator::face(MeshWorker::DoFInfo& dinfo1, // @sect3{The main class} + + // This class does the main job, like + // in previous examples. For a + // description of the functions + // declared here, please refer to + // the implementation below. template class Step39 { public: - typedef typename MeshWorker::IntegrationInfo CellInfo; + typedef MeshWorker::IntegrationInfo CellInfo; Step39(const FiniteElement& fe); @@ -346,27 +383,66 @@ class Step39 void solve (); void refine_grid (); void output_results (const unsigned int cycle) const; - + + // The member objects related to + // the discretization are here. Triangulation triangulation; const MappingQ1 mapping; const FiniteElement& fe; MGDoFHandler mg_dof_handler; DoFHandler& dof_handler; + // Then, we have the matrices and + // vectors related to the global + // discrete system. SparsityPattern sparsity; SparseMatrix matrix; Vector solution; Vector right_hand_side; BlockVector estimates; - + + // Finally, we have a group of + // sparsity patterns and sparse + // matrices related to the + // multilevel preconditioner. + // First, we have a level matrix + // and its sparsity pattern. MGLevelObject mg_sparsity; - MGLevelObject mg_sparsity_dg_interface; MGLevelObject > mg_matrix; - MGLevelObject > mg_matrix_dg_up; + + // When we perform multigrid with + // local smoothing on locally + // refined meshes, additional + // matrices are required; see + // Kanschat (2004). Here is the + // sparsity pattern for these + // edge matrices. We only need + // one, because the pattern of + // the up matrix is the + // transpose of that of the down + // matrix. Actually, we do not + // care too much about these + // details, since the MeshWorker + // is filling these matrices. + MGLevelObject mg_sparsity_dg_interface; + // The flux matrix at the + // refinement edge, coupling fine + // level degrees of freedom to + // coarse level. MGLevelObject > mg_matrix_dg_down; + // The transpose of the flux + // matrix at the refinement edge, + // coupling coarse level degrees + // of freedom to fine level. + MGLevelObject > mg_matrix_dg_up; }; + // The constructor simply sets up the + // coarse grid and the + // DoFHandler. The FiniteElement is + // provided as a parameter to allow + // flexibility. template Step39::Step39(const FiniteElement& fe) : @@ -379,49 +455,93 @@ Step39::Step39(const FiniteElement& fe) } + // In this function, we set up the + // dimension of the linear system and + // the sparsity patterns for the + // global matrix as well as the level + // matrices. template void Step39::setup_system() { + // First, we use the finite element + // to distribute degrees of + // freedom over the mesh and number + // them. dof_handler.distribute_dofs(fe); unsigned int n_dofs = dof_handler.n_dofs(); - + // Then, we already know the size + // of the vectors representing + // finite element functions. + solution.reinit(n_dofs); + right_hand_side.reinit(n_dofs); + + // Next, we set up the sparsity + // pattern for the global + // matrix. Since we do not know the + // row sizes in advance, we first + // fill a temporary + // CompressedSparsityPattern object + // and copy it to the regular + // SparsityPattern once it is + // complete. CompressedSparsityPattern c_sparsity(n_dofs); DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity); sparsity.copy_from(c_sparsity); matrix.reinit(sparsity); - - solution.reinit(n_dofs); - right_hand_side.reinit(n_dofs); const unsigned int n_levels = triangulation.n_levels(); - + // The global system is set up, now + // we attend to the level + // matrices. We resize all matrix + // objects to hold one matrix per level. mg_matrix.resize(0, n_levels-1); mg_matrix.clear(); mg_matrix_dg_up.resize(0, n_levels-1); mg_matrix_dg_up.clear(); mg_matrix_dg_down.resize(0, n_levels-1); mg_matrix_dg_down.clear(); - + // It is important to update the + // sparsity patterns after + // clear() was called for + // the level matrices, since the + // matrices lock the sparsity + // pattern through the Smartpointer + // ans Subscriptor mechanism. mg_sparsity.resize(0, n_levels-1); mg_sparsity_dg_interface.resize(0, n_levels-1); - + + // Now all objects are prepared to + // hold one sparsity pattern or + // matrix per level. What's left is + // setting up the sparsity patterns + // on each level. for (unsigned int level=mg_sparsity.get_minlevel(); level<=mg_sparsity.get_maxlevel();++level) { - CompressedSparsityPattern c_sparsity(mg_dof_handler.n_dofs(level)); - CompressedSparsityPattern ci_sparsity; - if (level>0) - ci_sparsity.reinit(mg_dof_handler.n_dofs(level-1), mg_dof_handler.n_dofs(level)); - + // These are roughly the same + // lines as above for the + // global matrix, now for each + // level. + CompressedSparsityPattern c_sparsity(mg_dof_handler.n_dofs(level)); MGTools::make_flux_sparsity_pattern(mg_dof_handler, c_sparsity, level); - if (level>0) - MGTools::make_flux_sparsity_pattern_edge(mg_dof_handler, ci_sparsity, level); - mg_sparsity[level].copy_from(c_sparsity); mg_matrix[level].reinit(mg_sparsity[level]); + + // Additionally, we need to + // initialize the transfer + // matrices at the refinement + // edge between levels. They + // are stored at the index + // referring to the finer of + // the two indices, thus there + // is no such object on level + // 0. if (level>0) { + CompressedSparsityPattern ci_sparsity; + ci_sparsity.reinit(mg_dof_handler.n_dofs(level-1), mg_dof_handler.n_dofs(level)); + MGTools::make_flux_sparsity_pattern_edge(mg_dof_handler, ci_sparsity, level); mg_sparsity_dg_interface[level].copy_from(ci_sparsity); mg_matrix_dg_up[level].reinit(mg_sparsity_dg_interface[level]); mg_matrix_dg_down[level].reinit(mg_sparsity_dg_interface[level]); @@ -430,23 +550,61 @@ Step39::setup_system() } + // In this function, we assemble the + // global system matrix, where by + // global we indicate that this is + // the matrix of the discrete system + // we solve and it is covering the + // whole mesh. template void Step39::assemble_matrix() { + // First, we need t set up the + // object providing the values we + // integrate. This object contains + // all FEValues and FEFaceValues + // objects needed and also + // maintains them automatically + // such that they always point to + // the current cell. To this end, + // we need to tell it first, where + // and what to compute, MeshWorker::IntegrationInfoBox info_box; + // namely, which quadrature + // formulas to use and const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points); + // which values to update in these + // points. We call + // initialize_update_flags + // first in order to set default + // values. Then, we add what we + // need additionally. info_box.initialize_update_flags(); UpdateFlags update_flags = update_values | update_gradients; info_box.add_update_flags(update_flags, true, true, true, true); info_box.initialize(fe, mapping); - + + // This is the object into which we + // integrate local data. MeshWorker::DoFInfo dof_info(dof_handler); - + + // Finally, we need an object that + // assembles the local matrix into + // the global matrix. MeshWorker::Assembler::MatrixSimple > assembler; assembler.initialize(matrix); - + + // Now, we throw everything into a + // MeshWorker::loop(), which here + // traverses all active cells of + // the mesh, computes cell and face + // matrices and assembles them into + // the global matrix. We use the + // variable dof_handler + // here in order to use the global + // numbering of degrees of freedom. MeshWorker::integration_loop( dof_handler.begin_active(), dof_handler.end(), dof_info, info_box, @@ -457,6 +615,11 @@ Step39::assemble_matrix() } + // Now, we do the same for the level + // matrices. Not too surprisingly, + // this function looks like a twin of + // the previous one. Indeed, there + // are only two minor differences. template void Step39::assemble_mg_matrix() @@ -471,10 +634,23 @@ Step39::assemble_mg_matrix() MeshWorker::DoFInfo dof_info(mg_dof_handler); + // Obviously, the assembler needs + // to be replaced by one filling + // level matrices. Note that it + // automatically fills the edge + // matrices as well. MeshWorker::Assembler::MGMatrixSimple > assembler; assembler.initialize(mg_matrix); assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down); + // Here is the other difference to + // the previous function: we run + // over all cells, not only the + // active ones. And we use + // mg_dof_handler, since + // we need the degrees of freedom + // on each level, not the global + // numbering. MeshWorker::integration_loop ( mg_dof_handler.begin(), mg_dof_handler.end(), dof_info, info_box, @@ -485,6 +661,11 @@ Step39::assemble_mg_matrix() } + // Here we have another clone of the + // assemble function. The difference + // to assembling the system matrix + // consists in that we assemble a + // vector here. template void Step39::assemble_right_hand_side() @@ -498,7 +679,17 @@ Step39::assemble_right_hand_side() info_box.initialize(fe, mapping); MeshWorker::DoFInfo dof_info(dof_handler); - + + // Since this assembler alows us to + // fill several vectors, the + // interface is a little more + // complicated as above. The + // pointers to the vectors have to + // be stored in a NamedData + // object. While this seems to + // cause two extra lines of code + // here, it actually comes handy in + // more complex applications. MeshWorker::Assembler::ResidualSimple > assembler; NamedData* > data; Vector* rhs = &right_hand_side; @@ -517,27 +708,52 @@ Step39::assemble_right_hand_side() } + // Now that we have coded all + // functions building the discrete + // linear system, it is about time + // that we actually solve it. template void Step39::solve() { + // The solver of choice is + // conjugate gradient. SolverControl control(1000, 1.e-12); SolverCG > cg(control); - GrowingVectorMemory > mem; + // Now we are setting up the + // components of the multilevel + // preconditioner. First, we need + // transfer between grid + // levels. The object we are using + // here generates sparse matrices + // for these transfers. MGTransferPrebuilt > mg_transfer; mg_transfer.build_matrices(mg_dof_handler); + + // Then, we need an exact solver + // for the matrix on the coarsest + // level. FullMatrix coarse_matrix; coarse_matrix.copy_from (mg_matrix[0]); MGCoarseGridHouseholder > mg_coarse; mg_coarse.initialize(coarse_matrix); - typedef PreconditionSSOR > RELAXATION; + + // While transfer and coarse grid + // solver are pretty much generic, + // more flexibility is offered for + // the smoother. First, we choose + // Gauss-Seidel as our smoothing + // method. + GrowingVectorMemory > mem; + typedef PreconditionSOR > RELAXATION; MGSmootherRelaxation, RELAXATION, Vector > mg_smoother(mem); RELAXATION::AdditionalData smoother_data(1.); mg_smoother.initialize(mg_matrix, smoother_data); - // Do two smoothing steps per level + // Do two smoothing steps on each + // level. mg_smoother.set_steps(2); // Since the SOR method is not // symmetric, but we use conjugate @@ -547,15 +763,18 @@ Step39::solve() // symmetric operator even for // nonsymmetric smoothers. mg_smoother.set_symmetric(true); + // The smoother class optionally + // implements the variable V-cycle, + // which we do not want here. mg_smoother.set_variable(false); - // We must wrap our matrices in an - // object having the required - // multiplication functions. + // Finally, we must wrap our + // matrices in an object having the + // required multiplication + // functions. MGMatrix, Vector > mgmatrix(&mg_matrix); MGMatrix, Vector > mgdown(&mg_matrix_dg_down); MGMatrix, Vector > mgup(&mg_matrix_dg_up); - // Now, we are ready to set up the // V-cycle operator and the @@ -563,18 +782,26 @@ Step39::solve() Multigrid > mg(mg_dof_handler, mgmatrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother); + // Let us not forget the edge + // matrices needed because of the + // adaptive refinement. mg.set_edge_flux_matrices(mgdown, mgup); - mg.set_debug(0); - mg_smoother.set_debug(0); - + + // After all preparations, wrap the + // Multigrid object into another + // object, which can be used as a + // regular preconditioner, PreconditionMG, MGTransferPrebuilt > > preconditioner(mg_dof_handler, mg, mg_transfer); - + // and use it to solve the system. cg.solve(matrix, solution, right_hand_side, preconditioner); } - + // Here we compare our finite element + // solution with the (known) exact + // solution and compute the mean + // quadratic error of the gradient. template void Step39::error() @@ -584,42 +811,71 @@ Step39::error() QGauss quadrature(n_gauss_points); VectorTools::integrate_difference(mapping, dof_handler, solution, exact_solution, - cell_errors, quadrature, VectorTools::H1_norm); + cell_errors, quadrature, VectorTools::H1_seminorm); deallog << "Error " << cell_errors.l2_norm() << std::endl; } + // Another clone of the assemble + // function. The big difference to + // the previous ones is here that we + // also have an input vector. template double Step39::estimate() { + // The results of the estimator are + // stored in a vector with one + // entry per cell. Since cells in + // deal.II are not numbered, we + // have to create our own numbering + // in order to use this vector. estimates.block(0).reinit(triangulation.n_active_cells()); unsigned int i=0; for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); cell != triangulation.end();++cell,++i) cell->set_user_index(i); - + + // This starts like before, MeshWorker::IntegrationInfoBox info_box; const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); - + + // but now we need to notify the + // info box of the finite element + // functio we want to evaluate in + // the quadrature points. First, we + // create a NamedData object with + // this vector, which is the + // solution we just computed. NamedData* > solution_data; solution_data.add(&solution, "solution"); - - MeshWorker::VectorSelector cs; - MeshWorker::VectorSelector fs; - cs.add("solution", true, true, true); - fs.add("solution", true, true, false); - - info_box.cell_selector = cs; - info_box.boundary_selector = fs; - info_box.face_selector = fs; + // Then, we tell the Meshworker::VectorSelector + // for cells, that we need the + // second derivatives of this + // solution (to compute the Laplacian). + info_box.cell_selector.add("solution", false, false, true); + // On interior and boundary faces, + // we need the function values and + // the first derivatives. + info_box.boundary_selector.add("solution", true, true, false); + info_box.face_selector.add("solution", true, true, false); + + // And we continue as before, with + // the exception that the default + // update flags are already + // adjusted to the values and + // derivatives we requested above. info_box.initialize_update_flags(); info_box.add_update_flags(update_quadrature_points, false, true, false, false); info_box.initialize(fe, mapping, solution_data); MeshWorker::DoFInfo dof_info(dof_handler); - + + // The assembler stores one number + // per cell, but else this is the + // same as in the computation of + // the right hand side. MeshWorker::Assembler::CellsAndFaces assembler; NamedData* > out_data; BlockVector* est = &estimates; @@ -637,6 +893,7 @@ Step39::estimate() } + // Some graphical output template void Step39::output_results (const unsigned int cycle) const { @@ -661,7 +918,9 @@ void Step39::output_results (const unsigned int cycle) const data_out.write_gnuplot(gnuplot_output); } - + // And finally the adaptive loop, + // more or less like in previous + // examples. template void Step39::run(unsigned int n_steps)