From: Timo Heister Date: Tue, 24 Jul 2018 15:42:57 +0000 (+0200) Subject: port step-16 to mesh_loop X-Git-Tag: v9.1.0-rc1~739^2~6 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0d5c0ff317168fd56bbdca6706a31945b0a40203;p=dealii.git port step-16 to mesh_loop --- diff --git a/examples/step-16/step-16.cc b/examples/step-16/step-16.cc index 5446333291..f28d688d3c 100644 --- a/examples/step-16/step-16.cc +++ b/examples/step-16/step-16.cc @@ -16,6 +16,7 @@ * Authors: Guido Kanschat, University of Heidelberg, 2003 * Baerbel Janssen, University of Heidelberg, 2010 * Wolfgang Bangerth, Texas A&M University, 2010 + * Timo Heister, Clemson University, 2018 */ @@ -68,20 +69,10 @@ #include #include -// Finally we include the MeshWorker framework. This framework through its -// function loop() and integration_loop(), automates loops over cells and -// assembling of data into vectors, matrices, etc. It obeys constraints -// automatically. Since we have to build several matrices and have to be aware -// of several sets of constraints, this will save us a lot of headache. -#include -#include -#include -#include -#include - -// In order to save effort, we use the pre-implemented Laplacian found in -#include -#include +// We will be using MeshWorker::mesh_loop to loop over the cells, so include it +// here: +#include + // This is C++: #include @@ -91,85 +82,59 @@ using namespace dealii; namespace Step16 { - // @sect3{The integrator on each cell} - - // The MeshWorker::integration_loop() expects a class that provides functions - // for integration on cells and boundary and interior faces. This is done by - // the following class. In the constructor, we tell the loop that cell - // integrals should be computed (the 'true'), but integrals should not be - // computed on boundary and interior faces (the two 'false'). Accordingly, we - // only need a cell function, but none for the faces. + // @sect3{The Scratch and Copy objects} + // + // We use MeshWorker::mesh_loop() to assemble our matrices. For this, we need + // a ScratchData object to store temporary data on each cell (this is just the + // FEValues object) and a CopyData object that will contain the output of each + // cell assembly. template - class LaplaceIntegrator : public MeshWorker::LocalIntegrator + struct ScratchData { - public: - LaplaceIntegrator(); - virtual void cell(MeshWorker::DoFInfo & dinfo, - MeshWorker::IntegrationInfo &info) const override; + ScratchData(const Mapping & mapping, + const FiniteElement &fe, + const unsigned int quadrature_degree, + const UpdateFlags update_flags) + : fe_values(mapping, fe, QGauss(quadrature_degree), update_flags) + {} + + ScratchData(const ScratchData &scratch_data) + : fe_values(scratch_data.fe_values.get_mapping(), + scratch_data.fe_values.get_fe(), + scratch_data.fe_values.get_quadrature(), + scratch_data.fe_values.get_update_flags()) + {} + + FEValues fe_values; }; - - template - LaplaceIntegrator::LaplaceIntegrator() - : MeshWorker::LocalIntegrator(true, false, false) - {} - - - // Next the actual integrator on each cell. We solve a Poisson problem with a - // coefficient one in the right half plane and one tenth in the left - // half plane. - - // The MeshWorker::LocalResults base class of MeshWorker::DoFInfo contains - // objects that can be filled in this local integrator. How many objects are - // created is determined inside the MeshWorker framework by the assembler - // class. Here, we test for instance that one matrix is required - // (MeshWorker::LocalResults::n_matrices()). The matrices are accessed through - // MeshWorker::LocalResults::matrix(), which takes the number of the matrix as - // its first argument. The second argument is only used for integrals over - // faces when there are two matrices for each test function used. Then, a - // second matrix with indicator 'true' would exist with the same index. - - // MeshWorker::IntegrationInfo provides one or several FEValues objects, which - // below are used by LocalIntegrators::Laplace::cell_matrix() or - // LocalIntegrators::L2::L2(). Since we are assembling only a single PDE, - // there is also only one of these objects with index zero. - - // In addition, we note that this integrator serves to compute the matrices - // for the multilevel preconditioner as well as the matrix and the right hand - // side for the global system. Since the assembler for a system requires an - // additional vector, MeshWorker::LocalResults::n_vectors() is returning a - // nonzero value. Accordingly, we fill a right hand side vector at the end of - // this function. Since LocalResults can deal with several BlockVector - // objects, but we are again in the simplest case here, we enter the - // information into block zero of vector zero. - template - void - LaplaceIntegrator::cell(MeshWorker::DoFInfo & dinfo, - MeshWorker::IntegrationInfo &info) const + struct CopyData { - AssertDimension(dinfo.n_matrices(), 1); - const double coefficient = (dinfo.cell->center()(0) > 0.) ? .1 : 1.; + unsigned int level; + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; - LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0, false).matrix, - info.fe_values(0), - coefficient); - - if (dinfo.n_vectors() > 0) - { - std::vector rhs(info.fe_values(0).n_quadrature_points, 1.); - LocalIntegrators::L2::L2(dinfo.vector(0).block(0), - info.fe_values(0), - rhs); - } - } + template + void reinit(const Iterator &cell, unsigned int dofs_per_cell) + { + cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_rhs.reinit(dofs_per_cell); + local_dof_indices.resize(dofs_per_cell); + cell->get_active_or_mg_dof_indices(local_dof_indices); + level = cell->level(); + } + }; // @sect3{The LaplaceProblem class template} - // This main class is basically the same class as in step-6. As far as - // member functions is concerned, the only addition is the - // assemble_multigrid function that assembles the matrices that - // correspond to the discrete operators on intermediate levels: + // This main class is similar to the same class in step-6. As far as + // member functions is concerned, the only additions are: + // - The assemble_multigrid function that assembles the matrices + // that correspond to the discrete operators on intermediate levels. + // - The cell_worker function that assembles our PDE on a single + // cell. template class LaplaceProblem { @@ -178,6 +143,11 @@ namespace Step16 void run(); private: + template + void cell_worker(const Iterator & cell, + ScratchData &scratch_data, + CopyData & copy_data); + void setup_system(); void assemble_system(); void assemble_multigrid(); @@ -200,9 +170,9 @@ namespace Step16 const unsigned int degree; // The following members are the essential data structures for the multigrid - // method. The first two represent the sparsity patterns and the matrices on - // individual levels of the multilevel hierarchy, very much like the objects - // for the global mesh above. + // method. The first four represent the sparsity patterns and the matrices + // on individual levels of the multilevel hierarchy, very much like the + // objects for the global mesh above. // // Then we have two new matrices only needed for multigrid methods with // local smoothing on adaptive meshes. They convey data between the interior @@ -213,10 +183,11 @@ namespace Step16 // level and information about indices lying on a refinement edge between // two different refinement levels. It thus serves a similar purpose as // AffineConstraints, but on each level. - MGLevelObject mg_sparsity_patterns; + MGLevelObject mg_sparsity_patterns; + MGLevelObject mg_interface_sparsity_patterns; + MGLevelObject> mg_matrices; - MGLevelObject> mg_interface_in; - MGLevelObject> mg_interface_out; + MGLevelObject> mg_interface_matrices; MGConstrainedDoFs mg_constrained_dofs; }; @@ -257,15 +228,13 @@ namespace Step16 dof_handler.distribute_dofs(fe); dof_handler.distribute_mg_dofs(); - deallog << " Number of degrees of freedom: " << dof_handler.n_dofs() - << " (by level: "; + std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs() + << " (by level: "; for (unsigned int level = 0; level < triangulation.n_levels(); ++level) - deallog << dof_handler.n_dofs(level) - << (level == triangulation.n_levels() - 1 ? ")" : ", "); - deallog << std::endl; + std::cout << dof_handler.n_dofs(level) + << (level == triangulation.n_levels() - 1 ? ")" : ", "); + std::cout << std::endl; - DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern(dof_handler, dsp); solution.reinit(dof_handler.n_dofs()); system_rhs.reinit(dof_handler.n_dofs()); @@ -278,18 +247,20 @@ namespace Step16 const std::map *> dirichlet_boundary_functions = { {types::boundary_id(0), &homogeneous_dirichlet_bc}}; - VectorTools::interpolate_boundary_values( - static_cast &>(dof_handler), - dirichlet_boundary_functions, - constraints); + VectorTools::interpolate_boundary_values(dof_handler, + dirichlet_boundary_functions, + constraints); constraints.close(); - constraints.condense(dsp); - sparsity_pattern.copy_from(dsp); + + { + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints); + sparsity_pattern.copy_from(dsp); + } system_matrix.reinit(sparsity_pattern); // The multigrid constraints have to be initialized. They need to know - // about the boundary values as well, so we pass the - // dirichlet_boundary here as well. + // where Dirichlet boundary conditions are prescribed. mg_constrained_dofs.clear(); mg_constrained_dofs.initialize(dof_handler); mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, @@ -306,143 +277,218 @@ namespace Step16 // their SparsityPattern before the can be destroyed upon resizing. const unsigned int n_levels = triangulation.n_levels(); - mg_interface_in.resize(0, n_levels - 1); - mg_interface_in.clear_elements(); - mg_interface_out.resize(0, n_levels - 1); - mg_interface_out.clear_elements(); + mg_interface_matrices.resize(0, n_levels - 1); mg_matrices.resize(0, n_levels - 1); - mg_matrices.clear_elements(); mg_sparsity_patterns.resize(0, n_levels - 1); + mg_interface_sparsity_patterns.resize(0, n_levels - 1); // Now, we have to provide a matrix on each level. To this end, we first use // the MGTools::make_sparsity_pattern function to generate a preliminary // compressed sparsity pattern on each level (see the @ref Sparsity module // for more information on this topic) and then copy it over to the one we - // really want. The next step is to initialize both kinds of level matrices - // with these sparsity patterns. + // really want. The next step is to initialize the interface matrices with + // the fitting sparsity pattern. // // It may be worth pointing out that the interface matrices only have // entries for degrees of freedom that sit at or next to the interface // between coarser and finer levels of the mesh. They are therefore even // sparser than the matrices on the individual levels of our multigrid - // hierarchy. If we were more concerned about memory usage (and possibly the - // speed with which we can multiply with these matrices), we should use - // separate and different sparsity patterns for these two kinds of matrices. + // hierarchy. Therefore, we use a function specifically build for this + // purpose to generate it. for (unsigned int level = 0; level < n_levels; ++level) { - DynamicSparsityPattern dsp(dof_handler.n_dofs(level), - dof_handler.n_dofs(level)); - MGTools::make_sparsity_pattern(dof_handler, dsp, level); + { + DynamicSparsityPattern dsp(dof_handler.n_dofs(level), + dof_handler.n_dofs(level)); + MGTools::make_sparsity_pattern(dof_handler, dsp, level); + + mg_sparsity_patterns[level].copy_from(dsp); + mg_matrices[level].reinit(mg_sparsity_patterns[level]); + } + { + DynamicSparsityPattern dsp(dof_handler.n_dofs(level), + dof_handler.n_dofs(level)); + MGTools::make_interface_sparsity_pattern(dof_handler, + mg_constrained_dofs, + dsp, + level); + mg_interface_sparsity_patterns[level].copy_from(dsp); + mg_interface_matrices[level].reinit( + mg_interface_sparsity_patterns[level]); + } + } + } + + + // @sect4{LaplaceProblem::cell_worker} - mg_sparsity_patterns[level].copy_from(dsp); + // The cell_worker function is used to assemble the matrix and right-hand side + // on the given cell. This function is used for the active cells to generate + // the system_matrix and on each level to build the level matrices. + // + // Note that we also assemble a right-hand side when called from + // assemble_multigrid() even though it is not used. + template + template + void LaplaceProblem::cell_worker(const Iterator & cell, + ScratchData &scratch_data, + CopyData & copy_data) + { + FEValues &fe_values = scratch_data.fe_values; + fe_values.reinit(cell); + + const unsigned int dofs_per_cell = fe_values.get_fe().dofs_per_cell; + const unsigned int n_q_points = fe_values.get_quadrature().size(); + + copy_data.reinit(cell, dofs_per_cell); + + const std::vector &JxW = fe_values.get_JxW_values(); + + for (unsigned int q = 0; q < n_q_points; ++q) + { + const double coefficient = + (fe_values.get_quadrature_points()[q][0] < 0.0) ? 1.0 : 0.1; + //(cell->center().square() < 0.5 * 0.5) ? 10.0:1.0; - mg_matrices[level].reinit(mg_sparsity_patterns[level]); - mg_interface_in[level].reinit(mg_sparsity_patterns[level]); - mg_interface_out[level].reinit(mg_sparsity_patterns[level]); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + copy_data.cell_matrix(i, j) += + coefficient * + (fe_values.shape_grad(i, q) * fe_values.shape_grad(j, q)) * + JxW[q]; + } + copy_data.cell_rhs(i) += 1.0 * fe_values.shape_value(i, q) * JxW[q]; + } } } + // @sect4{LaplaceProblem::assemble_system} - // The following function assembles the linear system on the finest level of - // the mesh. Since we want to reuse the code here for the level assembling - // below, we use the local integrator class LaplaceIntegrator and leave the - // loops to the MeshWorker framework. Thus, this function first sets up the - // objects necessary for this framework, namely - // - a MeshWorker::IntegrationInfoBox object, which will provide all the - // required data in quadrature points on the cell. This object can be seen - // as an extension of FEValues, providing a lot more useful information, - // - a MeshWorker::DoFInfo object, which on the one hand side extends the - // functionality of cell iterators, but also provides space for return - // values in its base class LocalResults, - // - an assembler, in this case for the whole system. The term 'simple' here - // refers to the fact that the global system does not have a block - // structure, - // - the local integrator, which implements the actual forms. - // - // After the loop has combined all of these into a matrix and a right hand - // side, there is one thing left to do: the assemblers leave matrix rows and - // columns of constrained degrees of freedom untouched. Therefore, we put a - // one on the diagonal to make the whole system well posed. The value one, or - // any fixed value has the advantage, that its effect on the spectrum of the - // matrix is easily understood. Since the corresponding eigenvectors form an - // invariant subspace, the value chosen does not affect the convergence of - // Krylov space solvers. + // The following function assembles the linear system on the active cells of + // the mesh. For this, we pass two lambda functions to the mesh_loop() + // function. The cell_worker function redirects to the class member function + // of the same name, while the copier is specific to this function and copies + // local matrix and vector to the corresponding global ones using the + // constraints. template void LaplaceProblem::assemble_system() { - MappingQ1 mapping; - MeshWorker::IntegrationInfoBox info_box; - UpdateFlags update_flags = - update_values | update_gradients | update_hessians; - info_box.add_update_flags_all(update_flags); - info_box.initialize(fe, mapping); - - MeshWorker::DoFInfo dof_info(dof_handler); - - MeshWorker::Assembler::SystemSimple, Vector> - assembler; - assembler.initialize(constraints); - assembler.initialize(system_matrix, system_rhs); - - LaplaceIntegrator matrix_integrator; - MeshWorker::integration_loop(dof_handler.begin_active(), - dof_handler.end(), - dof_info, - info_box, - matrix_integrator, - assembler); - - for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i) - if (constraints.is_constrained(i)) - system_matrix.set(i, i, 1.); + MappingQ1 mapping; + + typedef decltype(dof_handler.begin_active()) Iterator; + + auto cell_worker = [&](const Iterator & cell, + ScratchData &scratch_data, + CopyData & copy_data) { + this->cell_worker(cell, scratch_data, copy_data); + }; + + auto copier = [&](const CopyData &cd) { + this->constraints.distribute_local_to_global(cd.cell_matrix, + cd.cell_rhs, + cd.local_dof_indices, + system_matrix, + system_rhs); + }; + + const unsigned int n_gauss_points = degree + 1; + + ScratchData scratch_data(mapping, + fe, + n_gauss_points, + update_values | update_gradients | + update_JxW_values | + update_quadrature_points); + + MeshWorker::mesh_loop(dof_handler.begin_active(), + dof_handler.end(), + cell_worker, + copier, + scratch_data, + CopyData(), + MeshWorker::assemble_own_cells); } // @sect4{LaplaceProblem::assemble_multigrid} - // The next function is the one that builds the linear operators (matrices) + // The next function is the one that builds the matrices // that define the multigrid method on each level of the mesh. The integration // core is the same as above, but the loop below will go over all existing // cells instead of just the active ones, and the results must be entered into // the correct level matrices. Fortunately, MeshWorker hides most of that from // us, and thus the difference between this function and the previous lies // only in the setup of the assembler and the different iterators in the loop. - // Also, fixing up the matrices in the end is a little more complicated. + // + // We generate an AffineConstraints<> object template void LaplaceProblem::assemble_multigrid() { - MappingQ1 mapping; - MeshWorker::IntegrationInfoBox info_box; - UpdateFlags update_flags = - update_values | update_gradients | update_hessians; - info_box.add_update_flags_all(update_flags); - info_box.initialize(fe, mapping); - - MeshWorker::DoFInfo dof_info(dof_handler); - - MeshWorker::Assembler::MGMatrixSimple> assembler; - assembler.initialize(mg_constrained_dofs); - assembler.initialize(mg_matrices); - assembler.initialize_interfaces(mg_interface_in, mg_interface_out); - - LaplaceIntegrator matrix_integrator; - MeshWorker::integration_loop(dof_handler.begin_mg(), - dof_handler.end_mg(), - dof_info, - info_box, - matrix_integrator, - assembler); - - const unsigned int nlevels = triangulation.n_levels(); - for (unsigned int level = 0; level < nlevels; ++level) + MappingQ1 mapping; + const unsigned int n_levels = triangulation.n_levels(); + + std::vector> boundary_constraints(n_levels); + for (unsigned int level = 0; level < n_levels; ++level) { - for (unsigned int i = 0; i < dof_handler.n_dofs(level); ++i) - if (mg_constrained_dofs.is_boundary_index(level, i) || - mg_constrained_dofs.at_refinement_edge(level, i)) - mg_matrices[level].set(i, i, 1.); + IndexSet dofset; + DoFTools::extract_locally_relevant_level_dofs(dof_handler, + level, + dofset); + boundary_constraints[level].reinit(dofset); + boundary_constraints[level].add_lines( + mg_constrained_dofs.get_refinement_edge_indices(level)); + boundary_constraints[level].add_lines( + mg_constrained_dofs.get_boundary_indices(level)); + boundary_constraints[level].close(); } + + typedef decltype(dof_handler.begin_mg()) Iterator; + + auto cell_worker = [&](const Iterator & cell, + ScratchData &scratch_data, + CopyData & copy_data) { + this->cell_worker(cell, scratch_data, copy_data); + }; + + auto copier = [&](const CopyData &cd) { + boundary_constraints[cd.level].distribute_local_to_global( + cd.cell_matrix, cd.local_dof_indices, mg_matrices[cd.level]); + + const unsigned int dofs_per_cell = cd.local_dof_indices.size(); + + // TODO EXPLAIN: + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + if (mg_constrained_dofs.is_interface_matrix_entry( + cd.level, cd.local_dof_indices[i], cd.local_dof_indices[j])) + { + mg_interface_matrices[cd.level].add(cd.local_dof_indices[i], + cd.local_dof_indices[j], + cd.cell_matrix(i, j)); + } + }; + + const unsigned int n_gauss_points = degree + 1; + + ScratchData scratch_data(mapping, + fe, + n_gauss_points, + update_values | update_gradients | + update_JxW_values | + update_quadrature_points); + + MeshWorker::mesh_loop(dof_handler.begin_mg(), + dof_handler.end_mg(), + cell_worker, + copier, + scratch_data, + CopyData(), + MeshWorker::assemble_own_cells); } @@ -521,8 +567,8 @@ namespace Step16 // initialize both up and down versions of the operator with the matrices // we already built: mg::Matrix> mg_matrix(mg_matrices); - mg::Matrix> mg_interface_up(mg_interface_in); - mg::Matrix> mg_interface_down(mg_interface_out); + mg::Matrix> mg_interface_up(mg_interface_matrices); + mg::Matrix> mg_interface_down(mg_interface_matrices); // Now, we are ready to set up the V-cycle operator and the multilevel // preconditioner. @@ -541,6 +587,9 @@ namespace Step16 solution = 0; solver.solve(system_matrix, solution, system_rhs, preconditioner); + std::cout << " Number of CG iterations: " << solver_control.last_step() + << "\n" + << std::endl; constraints.distribute(solution); } @@ -551,10 +600,7 @@ namespace Step16 // The following two functions postprocess a solution once it is // computed. In particular, the first one refines the mesh at the beginning // of each cycle while the second one outputs results at the end of each - // such cycle. The functions are almost unchanged from those in step-6, with - // the exception of one minor difference: we generate output in VTK - // format, to use the more modern visualization programs available today - // compared to those that were available when step-6 was written. + // such cycle. The functions are almost unchanged from those in step-6. template void LaplaceProblem::refine_grid() { @@ -562,7 +608,7 @@ namespace Step16 KellyErrorEstimator::estimate( dof_handler, - QGauss(3), + QGauss(degree + 2), std::map *>(), solution, estimated_error_per_cell); @@ -600,18 +646,18 @@ namespace Step16 { for (unsigned int cycle = 0; cycle < 8; ++cycle) { - deallog << "Cycle " << cycle << std::endl; + std::cout << "Cycle " << cycle << std::endl; if (cycle == 0) { GridGenerator::hyper_ball(triangulation); - triangulation.refine_global(1); + triangulation.refine_global(2); } else refine_grid(); - deallog << " Number of active cells: " - << triangulation.n_active_cells() << std::endl; + std::cout << " Number of active cells: " + << triangulation.n_active_cells() << std::endl; setup_system(); @@ -634,8 +680,6 @@ int main() { using namespace Step16; - deallog.depth_console(2); - LaplaceProblem<2> laplace_problem(1); laplace_problem.run(); }