From 2607a7fd7edb506dd5105841593069d1f3aaa0ce Mon Sep 17 00:00:00 2001 From: janssen Date: Tue, 5 Jan 2010 09:55:39 +0000 Subject: [PATCH] change the smoother for the multigrid git-svn-id: https://svn.dealii.org/trunk@20284 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/step-42.cc | 363 +++++++++++++++++++++++++--- 1 file changed, 323 insertions(+), 40 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index a7feadc094..de3a193ba8 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -56,6 +56,14 @@ #include +#include +#include +#include +#include +#include +#include +#include +#include #include #include @@ -92,9 +100,13 @@ class StokesProblem private: void setup_dofs (); void assemble_system (); + void assemble_multigrid (); void solve_schur (); void solve_block (); + void solve_mg (); + void find_dofs_on_lower_level (std::vector > &lower_dofs, + std::vector > &boundary_dofs); void output_results (const unsigned int refinement_cycle) const; void refine_mesh (); @@ -102,7 +114,7 @@ class StokesProblem Triangulation triangulation; FESystem fe; - DoFHandler dof_handler; + MGDoFHandler dof_handler; ConstraintMatrix constraints; @@ -112,6 +124,12 @@ class StokesProblem BlockVector solution; BlockVector system_rhs; + MGLevelObject mg_constraints; + MGLevelObject mg_sparsity; + MGLevelObject > mg_matrices; + + MGLevelObject > mg_interface_matrices; + std::vector > mg_dofs_per_component; std_cxx1x::shared_ptr::type> A_preconditioner; }; @@ -350,7 +368,7 @@ void StokesProblem::setup_dofs () system_matrix.clear (); dof_handler.distribute_dofs (fe); - DoFRenumbering::Cuthill_McKee (dof_handler); +// DoFRenumbering::Cuthill_McKee (dof_handler); std::vector block_component (dim+1,0); block_component[dim] = 1; @@ -399,7 +417,9 @@ void StokesProblem::setup_dofs () csp.collect_sizes(); - DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false); + DoFTools::make_sparsity_pattern ( + static_cast&>(dof_handler), + csp, constraints, false); sparsity_pattern.copy_from (csp); } @@ -415,6 +435,28 @@ void StokesProblem::setup_dofs () system_rhs.block(0).reinit (n_u); system_rhs.block(1).reinit (n_p); system_rhs.collect_sizes (); + + //now setup stuff for mg + const unsigned int nlevels = triangulation.n_levels(); + + mg_matrices.resize(0, nlevels-1); + mg_matrices.clear (); + mg_interface_matrices.resize(0, nlevels-1); + mg_interface_matrices.clear (); + mg_sparsity.resize(0, nlevels-1); + + MGTools::count_dofs_per_block (dof_handler, mg_dofs_per_component); + + for (unsigned int level=0;level0) + mg_interface_matrices[level].reinit (mg_sparsity[level]); + } } @@ -456,7 +498,7 @@ void StokesProblem::assemble_system () std::vector div_phi_u (dofs_per_cell); std::vector phi_p (dofs_per_cell); - typename DoFHandler::active_cell_iterator + typename MGDoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) @@ -517,52 +559,221 @@ void StokesProblem::assemble_system () } +template +void StokesProblem::assemble_multigrid () +{ + QGauss quadrature_formula(degree+2); + FEValues fe_values (fe, quadrature_formula, + update_values | + update_quadrature_points | + update_JxW_values | + update_gradients); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + + + std::vector > phi_grads_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); + + std::vector > interface_dofs; + std::vector > boundary_interface_dofs; + for (unsigned int level = 0; level tmp (dof_handler.n_dofs(level)); + interface_dofs.push_back (tmp); + boundary_interface_dofs.push_back (tmp); + } + MGTools::extract_inner_interface_dofs (dof_handler, + interface_dofs, boundary_interface_dofs); + + typename FunctionMap::type dirichlet_boundary; + BoundaryValues dirichlet_bc; + dirichlet_boundary[0] = &dirichlet_bc; + + std::vector > boundary_indices(triangulation.n_levels()); + std::vector component_mask (dim+1, true); + component_mask[dim] = false; + MGTools::make_boundary_list (dof_handler, dirichlet_boundary, + boundary_indices, component_mask); + + std::vector boundary_constraints (triangulation.n_levels()); + std::vector boundary_interface_constraints (triangulation.n_levels()); + for (unsigned int level=0; level::cell_iterator + cell = dof_handler.begin(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + // Remember the level of the + // current cell. + const unsigned int level = cell->level(); + // Compute the values specified + // by update flags above. + fe_values.reinit (cell); + local_matrix = 0; + + for (unsigned int q=0; qget_mg_dof_indices (local_dof_indices); + boundary_constraints[level] + .distribute_local_to_global (local_matrix, + local_dof_indices, + mg_matrices[level]); + + for (unsigned int i=0; i void StokesProblem::solve_block () { - SparseMatrix pressure_mass_matrix; - pressure_mass_matrix.reinit(sparsity_pattern.block(1,1)); - pressure_mass_matrix.copy_from(system_matrix.block(1,1)); - system_matrix.block(1,1) = 0; - - SparseILU pmass_preconditioner; - pmass_preconditioner.initialize (pressure_mass_matrix, - SparseILU::AdditionalData()); - - InverseMatrix,SparseILU > - m_inverse (pressure_mass_matrix, pmass_preconditioner); + SparseMatrix pressure_mass_matrix; + pressure_mass_matrix.reinit(sparsity_pattern.block(1,1)); + pressure_mass_matrix.copy_from(system_matrix.block(1,1)); + system_matrix.block(1,1) = 0; - BlockSchurPreconditioner::type, - SparseILU > - preconditioner (system_matrix, m_inverse, *A_preconditioner); - - SolverControl solver_control (system_matrix.m(), - 1e-6*system_rhs.l2_norm()); - GrowingVectorMemory > vector_memory; - SolverGMRES >::AdditionalData gmres_data; - gmres_data.max_n_tmp_vectors = 100; - - SolverGMRES > gmres(solver_control, vector_memory, - gmres_data); - - gmres.solve(system_matrix, solution, system_rhs, - preconditioner); - - constraints.distribute (solution); - - std::cout << " " - << solver_control.last_step() - << " block GMRES iterations"; + SparseILU pmass_preconditioner; + pmass_preconditioner.initialize (pressure_mass_matrix, + SparseILU::AdditionalData()); + + InverseMatrix,SparseILU > + m_inverse (pressure_mass_matrix, pmass_preconditioner); + + BlockSchurPreconditioner::type, + SparseILU > + preconditioner (system_matrix, m_inverse, *A_preconditioner); + + SolverControl solver_control (system_matrix.m(), + 1e-6*system_rhs.l2_norm()); + GrowingVectorMemory > vector_memory; + SolverGMRES >::AdditionalData gmres_data; + gmres_data.max_n_tmp_vectors = 100; + + SolverGMRES > gmres(solver_control, vector_memory, + gmres_data); + + gmres.solve(system_matrix, solution, system_rhs, + preconditioner); + + constraints.distribute (solution); + + std::cout << " " + << solver_control.last_step() + << " block GMRES iterations"; } -template -void StokesProblem::solve_schur () + template +void StokesProblem::solve_mg () { + assemble_multigrid (); + typedef PreconditionMG, MGTransferPrebuilt > > + MGPREC; + + GrowingVectorMemory > mg_vector_memory; + + MGTransferPrebuilt > mg_transfer(constraints); + mg_transfer.build_matrices(dof_handler); + + FullMatrix mg_coarse_matrix; + mg_coarse_matrix.copy_from (mg_matrices[0]); + MGCoarseGridHouseholder > mg_coarse; + mg_coarse.initialize(mg_coarse_matrix); + + MGMatrix, BlockVector > + mg_matrix(&mg_matrices); + + typedef PreconditionSSOR > REL; + MGSmootherRelaxation, REL, BlockVector > + mg_smoother(mg_vector_memory); + REL::AdditionalData smoother_data; + mg_smoother.initialize(mg_matrices, smoother_data); + mg_smoother.set_steps(2); + Multigrid > mg(dof_handler, + mg_matrix, + mg_coarse, + mg_transfer, + mg_smoother, + mg_smoother); + mg.set_debug(3); + + MGPREC preconditioner(dof_handler, mg, mg_transfer); + + SolverControl solver_control (system_matrix.m(), + 1e-6*system_rhs.l2_norm()); + GrowingVectorMemory > vector_memory; + SolverGMRES >::AdditionalData gmres_data; + gmres_data.max_n_tmp_vectors = 100; + + SolverGMRES > gmres(solver_control, vector_memory, + gmres_data); + +// PreconditionIdentity precondition_identity; + gmres.solve(system_matrix, solution, system_rhs, + preconditioner); + //gmres.solve(system_matrix, solution, system_rhs, + // precondition_identity); + + constraints.distribute (solution); + + std::cout << std::endl << " " + << solver_control.last_step() + << " block GMRES iterations"; +} + + +template +void StokesProblem::solve_schur () +{ const InverseMatrix, typename InnerPreconditioner::type> A_inverse (system_matrix.block(0,0), *A_preconditioner); @@ -622,6 +833,78 @@ void StokesProblem::solve_schur () } +template +void StokesProblem::find_dofs_on_lower_level (std::vector > &dof_lower_level, + std::vector > &boundary_dof_lower_level) +{ + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int dofs_per_face = fe.dofs_per_face; + + std::vector local_dof_indices (dofs_per_cell); + std::vector face_dof_indices (dofs_per_face); + + typename MGDoFHandler::cell_iterator cell = dof_handler.begin(), + endc = dof_handler.end(); + + for (; cell!=endc; ++cell) + { + std::vector cell_dofs(dofs_per_cell); + std::vector boundary_cell_dofs(dofs_per_cell); + const unsigned int level = cell->level(); + cell->get_mg_dof_indices (local_dof_indices); + for (unsigned int face_nr=0; + face_nr::faces_per_cell; ++face_nr) + { + typename DoFHandler::face_iterator face = cell->face(face_nr); + if(!cell->at_boundary(face_nr)) + { + //interior face + typename MGDoFHandler::cell_iterator neighbor + = cell->neighbor(face_nr); + // Do refinement face + // from the coarse side + if (neighbor->level() < cell->level()) + { + for(unsigned int j=0; j::faces_per_cell; ++inner_face_nr) + { + const unsigned int neighbor_face_nr = (face_nr+inner_face_nr)%GeometryInfo::faces_per_cell; + if(!cell->at_boundary(neighbor_face_nr)) + { + //other face is interior + typename MGDoFHandler::cell_iterator neighbor + = cell->neighbor(neighbor_face_nr); + // Do refinement face + // from the coarse side + if (neighbor->level() < cell->level()) + { + for(unsigned int j=0; j void @@ -662,7 +945,7 @@ StokesProblem::refine_mesh () std::vector component_mask (dim+1, false); component_mask[dim] = true; - KellyErrorEstimator::estimate (dof_handler, + KellyErrorEstimator::estimate (static_cast&>(dof_handler), QGauss(degree+1), typename FunctionMap::type(), solution, @@ -724,7 +1007,7 @@ void StokesProblem::run () assemble_system (); std::cout << " Solving..." << std::flush; - solve_block (); + solve_mg (); output_results (refinement_cycle); -- 2.39.5