From: Timo Heister Date: Thu, 26 Jul 2018 12:46:20 +0000 (+0200) Subject: add test for examples/step-56 X-Git-Tag: v9.1.0-rc1~877^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7cf9676958961c8e6497facab1c90ff0ee8293f9;p=dealii.git add test for examples/step-56 --- diff --git a/tests/examples/CMakeLists.txt b/tests/examples/CMakeLists.txt new file mode 100644 index 0000000000..d33eafa20b --- /dev/null +++ b/tests/examples/CMakeLists.txt @@ -0,0 +1,4 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) +INCLUDE(../setup_testsubproject.cmake) +PROJECT(testsuite CXX) +DEAL_II_PICKUP_TESTS() diff --git a/tests/examples/step-56.cc b/tests/examples/step-56.cc new file mode 100644 index 0000000000..1d7bbca5fa --- /dev/null +++ b/tests/examples/step-56.cc @@ -0,0 +1,1165 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2016 - 2018 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.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + + * Author: Ryan Grove, Clemson University + * Timo Heister, Clemson University + */ + +// @sect3{Include files} + +#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 + +#include +#include + +// We need to include the following file to do timings: +#include + +// This includes the files necessary for us to use geometric Multigrid +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace Step56 +{ + using namespace dealii; + + // In order to make it easy to switch between the different solvers that are + // being used, we declare an enum that can be passed as an argument to the + // constructor of the main class. + struct SolverType + { + enum type + { + FGMRES_ILU, + FGMRES_GMG, + UMFPACK + }; + }; + + // @sect3{Functions for Solution and Righthand side} + // + // The class Solution is used to define the boundary conditions and to + // compute errors of the numerical solution. Note that we need to define the + // values and gradients in order to compute L2 and H1 errors. Here we + // decided to separate the implementations for 2d and 3d using template + // specialization. + // + // Note that the first dim components are the velocity components + // and the last is the pressure. + template + class Solution : public Function + { + public: + Solution() + : Function(dim + 1) + {} + virtual double value(const Point & p, + const unsigned int component = 0) const override; + virtual Tensor<1, dim> + gradient(const Point & p, + const unsigned int component = 0) const override; + }; + + template <> + double Solution<2>::value(const Point<2> & p, + const unsigned int component) const + { + Assert(component <= 2 + 1, ExcIndexRange(component, 0, 2 + 1)); + + using numbers::PI; + const double x = p(0); + const double y = p(1); + + if (component == 0) + return sin(PI * x); + if (component == 1) + return -PI * y * cos(PI * x); + if (component == 2) + return sin(PI * x) * cos(PI * y); + + return 0; + } + + template <> + double Solution<3>::value(const Point<3> & p, + const unsigned int component) const + { + Assert(component <= 3 + 1, ExcIndexRange(component, 0, 3 + 1)); + + using numbers::PI; + const double x = p(0); + const double y = p(1); + const double z = p(2); + + if (component == 0) + return 2.0 * sin(PI * x); + if (component == 1) + return -PI * y * cos(PI * x); + if (component == 2) + return -PI * z * cos(PI * x); + if (component == 3) + return sin(PI * x) * cos(PI * y) * sin(PI * z); + + return 0; + } + + // Note that for the gradient we need to return a Tensor<1,dim> + template <> + Tensor<1, 2> Solution<2>::gradient(const Point<2> & p, + const unsigned int component) const + { + Assert(component <= 2, ExcIndexRange(component, 0, 2 + 1)); + + using numbers::PI; + const double x = p(0); + const double y = p(1); + + Tensor<1, 2> return_value; + if (component == 0) + { + return_value[0] = PI * cos(PI * x); + return_value[1] = 0.0; + } + else if (component == 1) + { + return_value[0] = y * PI * PI * sin(PI * x); + return_value[1] = -PI * cos(PI * x); + } + else if (component == 2) + { + return_value[0] = PI * cos(PI * x) * cos(PI * y); + return_value[1] = -PI * sin(PI * x) * sin(PI * y); + } + + return return_value; + } + + template <> + Tensor<1, 3> Solution<3>::gradient(const Point<3> & p, + const unsigned int component) const + { + Assert(component <= 3, ExcIndexRange(component, 0, 3 + 1)); + + using numbers::PI; + const double x = p(0); + const double y = p(1); + const double z = p(2); + + Tensor<1, 3> return_value; + if (component == 0) + { + return_value[0] = 2 * PI * cos(PI * x); + return_value[1] = 0.0; + return_value[2] = 0.0; + } + else if (component == 1) + { + return_value[0] = y * PI * PI * sin(PI * x); + return_value[1] = -PI * cos(PI * x); + return_value[2] = 0.0; + } + else if (component == 2) + { + return_value[0] = z * PI * PI * sin(PI * x); + return_value[1] = 0.0; + return_value[2] = -PI * cos(PI * x); + } + else if (component == 3) + { + return_value[0] = PI * cos(PI * x) * cos(PI * y) * sin(PI * z); + return_value[1] = -PI * sin(PI * x) * sin(PI * y) * sin(PI * z); + return_value[2] = PI * sin(PI * x) * cos(PI * y) * cos(PI * z); + } + + return return_value; + } + + // Implementation of $f$. See the introduction for more information. + template + class RightHandSide : public Function + { + public: + RightHandSide() + : Function(dim + 1) + {} + + virtual double value(const Point & p, + const unsigned int component = 0) const override; + }; + + template <> + double RightHandSide<2>::value(const Point<2> & p, + const unsigned int component) const + { + Assert(component <= 2, ExcIndexRange(component, 0, 2 + 1)); + + using numbers::PI; + double x = p(0); + double y = p(1); + if (component == 0) + return PI * PI * sin(PI * x) + PI * cos(PI * x) * cos(PI * y); + if (component == 1) + return -PI * PI * PI * y * cos(PI * x) - PI * sin(PI * y) * sin(PI * x); + if (component == 2) + return 0; + + return 0; + } + + template <> + double RightHandSide<3>::value(const Point<3> & p, + const unsigned int component) const + { + Assert(component <= 3, ExcIndexRange(component, 0, 3 + 1)); + + using numbers::PI; + double x = p(0); + double y = p(1); + double z = p(2); + if (component == 0) + return 2 * PI * PI * sin(PI * x) + + PI * cos(PI * x) * cos(PI * y) * sin(PI * z); + if (component == 1) + return -PI * PI * PI * y * cos(PI * x) + + PI * (-1) * sin(PI * y) * sin(PI * x) * sin(PI * z); + if (component == 2) + return -PI * PI * PI * z * cos(PI * x) + + PI * cos(PI * z) * sin(PI * x) * cos(PI * y); + if (component == 3) + return 0; + + return 0; + } + + + + // @sect3{ASPECT BlockSchurPreconditioner} + + // In the following, we will implement a preconditioner that expands + // on the ideas discussed in the Results section of step-22. + // Specifically, we + // 1. use an upper block-triangular preconditioner because we want to + // use right preconditioning. + // 2. optionally allow using an inner solver for the velocity block instead + // of a single preconditioner application. + // 3. do not use InverseMatrix but explicitly call SolverCG. + // This approach is also used in the ASPECT code + // (see http://aspect.dealii.org) that solves the Stokes equations in + // the context of simulating convection in the earth mantle, and which + // has been used to solve problems on many thousands of processors. + // + // The bool flag @p do_solve_A in the constructor allows us to either + // apply the preconditioner for the velocity block once or use an inner + // iterative solver for a more accurate approximation instead. + // + // Notice how we keep track of the sum of the inner iterations + // (preconditioner applications). + template + class BlockSchurPreconditioner : public Subscriptor + { + public: + BlockSchurPreconditioner( + const BlockSparseMatrix &system_matrix, + const SparseMatrix & schur_complement_matrix, + const PreconditionerAType & preconditioner_A, + const PreconditionerSType & preconditioner_S, + const bool do_solve_A); + + void vmult(BlockVector &dst, const BlockVector &src) const; + + mutable unsigned int n_iterations_A; + mutable unsigned int n_iterations_S; + + private: + const BlockSparseMatrix &system_matrix; + const SparseMatrix & schur_complement_matrix; + const PreconditionerAType & preconditioner_A; + const PreconditionerSType & preconditioner_S; + + const bool do_solve_A; + }; + + template + BlockSchurPreconditioner:: + BlockSchurPreconditioner( + const BlockSparseMatrix &system_matrix, + const SparseMatrix & schur_complement_matrix, + const PreconditionerAType & preconditioner_A, + const PreconditionerSType & preconditioner_S, + const bool do_solve_A) + : n_iterations_A(0) + , n_iterations_S(0) + , system_matrix(system_matrix) + , schur_complement_matrix(schur_complement_matrix) + , preconditioner_A(preconditioner_A) + , preconditioner_S(preconditioner_S) + , do_solve_A(do_solve_A) + {} + + + + template + void + BlockSchurPreconditioner::vmult( + BlockVector & dst, + const BlockVector &src) const + { + Vector utmp(src.block(0)); + + // First solve with the approximation for S + { + SolverControl solver_control(1000, 1e-6 * src.block(1).l2_norm()); + SolverCG<> cg(solver_control); + + dst.block(1) = 0.0; + cg.solve(schur_complement_matrix, + dst.block(1), + src.block(1), + preconditioner_S); + + n_iterations_S += solver_control.last_step(); + dst.block(1) *= -1.0; + } + + // Second, apply the top right block (B^T) + { + system_matrix.block(0, 1).vmult(utmp, dst.block(1)); + utmp *= -1.0; + utmp += src.block(0); + } + + // Finally, either solve with the top left block + // or just apply one preconditioner sweep + if (do_solve_A == true) + { + SolverControl solver_control(10000, utmp.l2_norm() * 1e-4); + SolverCG<> cg(solver_control); + + dst.block(0) = 0.0; + cg.solve(system_matrix.block(0, 0), + dst.block(0), + utmp, + preconditioner_A); + + n_iterations_A += solver_control.last_step(); + } + else + { + preconditioner_A.vmult(dst.block(0), utmp); + n_iterations_A += 1; + } + } + + // @sect3{The StokesProblem class} + // + // This is the main class of the problem. + template + class StokesProblem + { + public: + StokesProblem(const unsigned int pressure_degree, + SolverType::type solver_type); + void run(); + + private: + void setup_dofs(); + void assemble_system(); + void assemble_multigrid(); + void solve(); + void compute_errors(); + void output_results(const unsigned int refinement_cycle) const; + + const unsigned int pressure_degree; + SolverType::type solver_type; + + Triangulation triangulation; + FESystem velocity_fe; + FESystem fe; + DoFHandler dof_handler; + DoFHandler velocity_dof_handler; + + ConstraintMatrix constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + SparseMatrix pressure_mass_matrix; + + BlockVector solution; + BlockVector system_rhs; + + MGLevelObject mg_sparsity_patterns; + MGLevelObject> mg_matrices; + MGLevelObject> mg_interface_matrices; + MGConstrainedDoFs mg_constrained_dofs; + + TimerOutput computing_timer; + }; + + + + template + StokesProblem::StokesProblem(const unsigned int pressure_degree, + SolverType::type solver_type) + : pressure_degree(pressure_degree) + , solver_type(solver_type) + , triangulation(Triangulation::maximum_smoothing) + , + // Finite element for the velocity only: + velocity_fe(FE_Q(pressure_degree + 1), dim) + , + // Finite element for the whole system: + fe(velocity_fe, 1, FE_Q(pressure_degree), 1) + , dof_handler(triangulation) + , velocity_dof_handler(triangulation) + , computing_timer(std::cout, TimerOutput::never, TimerOutput::wall_times) + {} + + + + // @sect4{StokesProblem::setup_dofs} + + // This function sets up the DoFHandler, matrices, vectors, and Multigrid + // structures (if needed). + template + void StokesProblem::setup_dofs() + { + TimerOutput::Scope scope(computing_timer, "Setup"); + + system_matrix.clear(); + pressure_mass_matrix.clear(); + + // The main DoFHandler only needs active DoFs, so we are not calling + // distribute_mg_dofs() here + dof_handler.distribute_dofs(fe); + + // This block structure separates the dim velocity components from + // the pressure component (used for reordering). Note that we have + // 2 instead of dim+1 blocks like in step-22, because our FESystem + // is nested and the dim velocity components appear as one block. + std::vector block_component(2); + block_component[0] = 0; + block_component[1] = 1; + + // Velocities start at component 0: + const FEValuesExtractors::Vector velocities(0); + + // ILU behaves better if we apply a reordering to reduce fillin. There + // is no advantage in doing this for the other solvers. + if (solver_type == SolverType::FGMRES_ILU) + { + TimerOutput::Scope ilu_specific(computing_timer, "(ILU specific)"); + DoFRenumbering::Cuthill_McKee(dof_handler); + } + + // This ensures that all velocities DoFs are enumerated before the + // pressure unknowns. This allows us to use blocks for vectors and + // matrices and allows us to get the same DoF numbering for + // dof_handler and velocity_dof_handler. + DoFRenumbering::block_wise(dof_handler); + + if (solver_type == SolverType::FGMRES_GMG) + { + TimerOutput::Scope multigrid_specific(computing_timer, + "(Multigrid specific)"); + TimerOutput::Scope setup_multigrid(computing_timer, + "Setup - Multigrid"); + + // This distributes the active dofs and multigrid dofs for the + // velocity space in a separate DoFHandler as described in the + // introduction. + velocity_dof_handler.distribute_dofs(velocity_fe); + velocity_dof_handler.distribute_mg_dofs(); + + // The following block of code initializes the MGConstrainedDofs + // (using the boundary conditions for the velocity), and the + // sparsity patterns and matrices for each level. The resize() + // function of MGLevelObject will destroy all existing contained + // objects. + std::set zero_boundary_ids; + zero_boundary_ids.insert(0); + + mg_constrained_dofs.clear(); + mg_constrained_dofs.initialize(velocity_dof_handler); + mg_constrained_dofs.make_zero_boundary_constraints(velocity_dof_handler, + zero_boundary_ids); + const unsigned int n_levels = triangulation.n_levels(); + + mg_interface_matrices.resize(0, n_levels - 1); + mg_matrices.resize(0, n_levels - 1); + mg_sparsity_patterns.resize(0, n_levels - 1); + + for (unsigned int level = 0; level < n_levels; ++level) + { + DynamicSparsityPattern csp(velocity_dof_handler.n_dofs(level), + velocity_dof_handler.n_dofs(level)); + MGTools::make_sparsity_pattern(velocity_dof_handler, csp, level); + mg_sparsity_patterns[level].copy_from(csp); + + mg_matrices[level].reinit(mg_sparsity_patterns[level]); + mg_interface_matrices[level].reinit(mg_sparsity_patterns[level]); + } + } + + std::vector dofs_per_block(2); + DoFTools::count_dofs_per_block(dof_handler, + dofs_per_block, + block_component); + const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1]; + + { + constraints.clear(); + // The following makes use of a component mask for interpolation of the + // boundary values for the velocity only, which is further explained in + // the vector valued dealii step-20 tutorial. + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Solution(), + constraints, + fe.component_mask(velocities)); + + // As discussed in the introduction, we need to fix one degree of freedom + // of the pressure variable to ensure solvability of the problem. We do + // this here by marking the first pressure dof, which has index n_u as a + // constrained dof. + if (solver_type == SolverType::UMFPACK) + constraints.add_line(n_u); + + constraints.close(); + } + + std::cout << "\tNumber of active cells: " << triangulation.n_active_cells() + << std::endl + << "\tNumber of degrees of freedom: " << dof_handler.n_dofs() + << " (" << n_u << '+' << n_p << ')' << std::endl; + + { + BlockDynamicSparsityPattern csp(dofs_per_block, dofs_per_block); + DoFTools::make_sparsity_pattern(dof_handler, csp, constraints, false); + sparsity_pattern.copy_from(csp); + } + system_matrix.reinit(sparsity_pattern); + + solution.reinit(dofs_per_block); + system_rhs.reinit(dofs_per_block); + } + + + // @sect4{StokesProblem::assemble_system} + + // In this function, the system matrix is assembled. We assemble the pressure + // mass matrix in the (1,1) block (if needed) and move it out of this location + // at the end of this function. + template + void StokesProblem::assemble_system() + { + TimerOutput::Scope assemble(computing_timer, "Assemble"); + system_matrix = 0; + system_rhs = 0; + + // If true, we will assemble the pressure mass matrix in the (1,1) block: + const bool assemble_pressure_mass_matrix = + (solver_type == SolverType::UMFPACK) ? false : true; + + QGauss quadrature_formula(pressure_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); + Vector local_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + const RightHandSide right_hand_side; + std::vector> rhs_values(n_q_points, Vector(dim + 1)); + + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + + std::vector> symgrad_phi_u(dofs_per_cell); + std::vector div_phi_u(dofs_per_cell); + std::vector phi_p(dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell != endc; ++cell) + { + fe_values.reinit(cell); + local_matrix = 0; + local_rhs = 0; + + right_hand_side.vector_value_list(fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + symgrad_phi_u[k] = + fe_values[velocities].symmetric_gradient(k, q); + div_phi_u[k] = fe_values[velocities].divergence(k, q); + phi_p[k] = fe_values[pressure].value(k, q); + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j <= i; ++j) + { + local_matrix(i, j) += + (2 * (symgrad_phi_u[i] * symgrad_phi_u[j]) - + div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j] + + (assemble_pressure_mass_matrix ? phi_p[i] * phi_p[j] : + 0)) * + fe_values.JxW(q); + } + + const unsigned int component_i = + fe.system_to_component_index(i).first; + local_rhs(i) += fe_values.shape_value(i, q) * + rhs_values[q](component_i) * fe_values.JxW(q); + } + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = i + 1; j < dofs_per_cell; ++j) + local_matrix(i, j) = local_matrix(j, i); + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(local_matrix, + local_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } + + if (solver_type != SolverType::UMFPACK) + { + 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; + } + } + + // @sect4{StokesProblem::assemble_multigrid} + + // Here, like in step-16, we have a function that assembles the level + // and interface matrices necessary for the multigrid preconditioner. + template + void StokesProblem::assemble_multigrid() + { + TimerOutput::Scope multigrid_specific(computing_timer, + "(Multigrid specific)"); + TimerOutput::Scope assemble_multigrid(computing_timer, + "Assemble Multigrid"); + + mg_matrices = 0.; + + QGauss quadrature_formula(pressure_degree + 2); + + FEValues fe_values(velocity_fe, + quadrature_formula, + update_values | update_quadrature_points | + update_JxW_values | update_gradients); + + const unsigned int dofs_per_cell = velocity_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 FEValuesExtractors::Vector velocities(0); + + std::vector> symgrad_phi_u(dofs_per_cell); + + std::vector boundary_constraints( + triangulation.n_levels()); + std::vector boundary_interface_constraints( + triangulation.n_levels()); + for (unsigned int level = 0; level < triangulation.n_levels(); ++level) + { + 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(); + + IndexSet idx = mg_constrained_dofs.get_refinement_edge_indices(level) & + mg_constrained_dofs.get_boundary_indices(level); + + boundary_interface_constraints[level].add_lines(idx); + boundary_interface_constraints[level].close(); + } + + // This iterator goes over all cells (not just active) + typename DoFHandler::cell_iterator cell = velocity_dof_handler.begin(), + endc = velocity_dof_handler.end(); + + for (; cell != endc; ++cell) + { + fe_values.reinit(cell); + cell_matrix = 0; + + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int k = 0; k < dofs_per_cell; ++k) + symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient(k, q); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j <= i; ++j) + { + cell_matrix(i, j) += + (symgrad_phi_u[i] * symgrad_phi_u[j]) * fe_values.JxW(q); + } + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = i + 1; j < dofs_per_cell; ++j) + cell_matrix(i, j) = cell_matrix(j, i); + + cell->get_mg_dof_indices(local_dof_indices); + + boundary_constraints[cell->level()].distribute_local_to_global( + cell_matrix, local_dof_indices, mg_matrices[cell->level()]); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + if (!mg_constrained_dofs.at_refinement_edge(cell->level(), + local_dof_indices[i]) || + mg_constrained_dofs.at_refinement_edge(cell->level(), + local_dof_indices[j])) + cell_matrix(i, j) = 0; + + boundary_interface_constraints[cell->level()] + .distribute_local_to_global(cell_matrix, + local_dof_indices, + mg_interface_matrices[cell->level()]); + } + } + + // @sect4{StokesProblem::solve} + + // This function sets up things differently based on if you want to use ILU + // or GMG as a preconditioner. Both methods share the same solver (FGMRES) + // but require a different preconditioner to be initialized. Here we time not + // only the entire solve function, but we separately time the setup of the + // preconditioner as well as the solve itself. + template + void StokesProblem::solve() + { + TimerOutput::Scope solve(computing_timer, "Solve"); + constraints.set_zero(solution); + + if (solver_type == SolverType::UMFPACK) + { + computing_timer.enter_subsection("(UMFPACK specific)"); + computing_timer.enter_subsection("Solve - Initialize"); + + SparseDirectUMFPACK A_direct; + A_direct.initialize(system_matrix); + + computing_timer.leave_subsection(); + computing_timer.leave_subsection(); + + { + TimerOutput::Scope solve_backslash(computing_timer, + "Solve - Backslash"); + A_direct.vmult(solution, system_rhs); + } + + constraints.distribute(solution); + return; + } + + // Here we must make sure to solve for the residual with "good enough" + // accuracy + SolverControl solver_control(system_matrix.m(), + 1e-10 * system_rhs.l2_norm()); + unsigned int n_iterations_A; + unsigned int n_iterations_S; + + // This is used to pass whether or not we want to solve for A inside + // the preconditioner. One could change this to false to see if + // there is still convergence and if so does the program then run + // faster or slower + const bool use_expensive = true; + + SolverFGMRES> solver(solver_control); + + if (solver_type == SolverType::FGMRES_ILU) + { + computing_timer.enter_subsection("(ILU specific)"); + computing_timer.enter_subsection("Solve - Set-up Preconditioner"); + + std::cout << " Computing preconditioner..." << std::endl + << std::flush; + + SparseILU A_preconditioner; + A_preconditioner.initialize(system_matrix.block(0, 0)); + + SparseILU S_preconditioner; + S_preconditioner.initialize(pressure_mass_matrix); + + const BlockSchurPreconditioner, SparseILU> + preconditioner(system_matrix, + pressure_mass_matrix, + A_preconditioner, + S_preconditioner, + use_expensive); + + computing_timer.leave_subsection(); + computing_timer.leave_subsection(); + + { + TimerOutput::Scope solve_fmgres(computing_timer, "Solve - FGMRES"); + + solver.solve(system_matrix, solution, system_rhs, preconditioner); + n_iterations_A = preconditioner.n_iterations_A; + n_iterations_S = preconditioner.n_iterations_S; + } + } + else + { + computing_timer.enter_subsection("(Multigrid specific)"); + computing_timer.enter_subsection("Solve - Set-up Preconditioner"); + + // Transfer operators between levels + MGTransferPrebuilt> mg_transfer(mg_constrained_dofs); + mg_transfer.build_matrices(velocity_dof_handler); + + // Setup coarse grid solver + FullMatrix coarse_matrix; + coarse_matrix.copy_from(mg_matrices[0]); + MGCoarseGridHouseholder<> coarse_grid_solver; + coarse_grid_solver.initialize(coarse_matrix); + + using Smoother = PreconditionSOR>; + mg::SmootherRelaxation> mg_smoother; + mg_smoother.initialize(mg_matrices); + mg_smoother.set_steps(2); + + // Multigrid, when used as a preconditioner for CG, needs to be a + // symmetric operator, so the smoother must be symmetric + mg_smoother.set_symmetric(true); + + mg::Matrix> mg_matrix(mg_matrices); + 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. + Multigrid> mg( + mg_matrix, coarse_grid_solver, mg_transfer, mg_smoother, mg_smoother); + mg.set_edge_matrices(mg_interface_down, mg_interface_up); + + PreconditionMG, MGTransferPrebuilt>> + A_Multigrid(velocity_dof_handler, mg, mg_transfer); + + SparseILU S_preconditioner; + S_preconditioner.initialize(pressure_mass_matrix, + SparseILU::AdditionalData()); + + const BlockSchurPreconditioner< + PreconditionMG, + MGTransferPrebuilt>>, + SparseILU> + preconditioner(system_matrix, + pressure_mass_matrix, + A_Multigrid, + S_preconditioner, + use_expensive); + + computing_timer.leave_subsection(); + computing_timer.leave_subsection(); + + { + TimerOutput::Scope solve_fmgres(computing_timer, "Solve - FGMRES"); + solver.solve(system_matrix, solution, system_rhs, preconditioner); + n_iterations_A = preconditioner.n_iterations_A; + n_iterations_S = preconditioner.n_iterations_S; + } + } + + constraints.distribute(solution); + + std::cout + << std::endl + << "\tNumber of FGMRES iterations: " << solver_control.last_step() + << std::endl + << "\tTotal number of iterations used for approximation of A inverse: " + << n_iterations_A << std::endl + << "\tTotal number of iterations used for approximation of S inverse: " + << n_iterations_S << std::endl + << std::endl; + } + + + // @sect4{StokesProblem::process_solution} + + // This function computes the L2 and H1 errors of the solution. For this, + // we need to make sure the pressure has mean zero. + template + void StokesProblem::compute_errors() + { + // Compute the mean pressure $\frac{1}{\Omega} \int_{\Omega} p(x) dx $ + // and then subtract it from each pressure coefficient. This will result + // in a pressure with mean value zero. Here we make use of the fact that + // the pressure is component $dim$ and that the finite element space + // is nodal. + const double mean_pressure = VectorTools::compute_mean_value( + dof_handler, QGauss(pressure_degree + 2), solution, dim); + solution.block(1).add(-mean_pressure); + std::cout << " Note: The mean value was adjusted by " << -mean_pressure + << std::endl; + + const ComponentSelectFunction pressure_mask(dim, dim + 1); + const ComponentSelectFunction velocity_mask(std::make_pair(0, dim), + dim + 1); + + Vector difference_per_cell(triangulation.n_active_cells()); + VectorTools::integrate_difference(dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss(pressure_degree + 2), + VectorTools::L2_norm, + &velocity_mask); + + const double Velocity_L2_error = + VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::L2_norm); + + VectorTools::integrate_difference(dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss(pressure_degree + 2), + VectorTools::L2_norm, + &pressure_mask); + + const double Pressure_L2_error = + VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::L2_norm); + + VectorTools::integrate_difference(dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss(pressure_degree + 2), + VectorTools::H1_norm, + &velocity_mask); + + const double Velocity_H1_error = + VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::H1_norm); + + std::cout << std::endl + << " Velocity L2 Error: " << Velocity_L2_error << std::endl + << " Pressure L2 Error: " << Pressure_L2_error << std::endl + << " Velocity H1 Error: " << Velocity_H1_error << std::endl; + } + + + // @sect4{StokesProblem::output_results} + + // This function generates graphical output like it is done in step-22. + template + void + StokesProblem::output_results(const unsigned int refinement_cycle) const + { + std::vector solution_names(dim, "velocity"); + solution_names.emplace_back("pressure"); + + std::vector + data_component_interpretation( + dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, + solution_names, + DataOut::type_dof_data, + data_component_interpretation); + data_out.build_patches(); + + std::ofstream output( + "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtk"); + data_out.write_vtk(output); + } + + + + // @sect4{StokesProblem::run} + + // The last step in the Stokes class is, as usual, the function that + // generates the initial grid and calls the other functions in the + // respective order. + template + void StokesProblem::run() + { + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(2); + + if (solver_type == SolverType::FGMRES_ILU) + std::cout << "Now running with ILU" << std::endl; + else if (solver_type == SolverType::FGMRES_GMG) + std::cout << "Now running with Multigrid" << std::endl; + else + std::cout << "Now running with UMFPACK" << std::endl; + + + for (unsigned int refinement_cycle = 0; refinement_cycle < 2; + ++refinement_cycle) + { + std::cout << "Refinement cycle " << refinement_cycle << std::endl; + + if (refinement_cycle > 0) + triangulation.refine_global(1); + + std::cout << " Set-up..." << std::endl; + setup_dofs(); + + std::cout << " Assembling..." << std::endl; + assemble_system(); + + if (solver_type == SolverType::FGMRES_GMG) + { + std::cout << " Assembling Multigrid..." << std::endl; + + assemble_multigrid(); + } + + std::cout << " Solving..." << std::flush; + solve(); + + compute_errors(); + + output_results(refinement_cycle); + + computing_timer.reset(); + } + } +} // namespace Step56 + +// @sect3{The main function} +int main() +{ + try + { + using namespace dealii; + using namespace Step56; + + const int degree = 1; + const int dim = 3; + // options for SolverType: UMFPACK FGMRES_ILU FGMRES_GMG + { + + StokesProblem flow_problem(degree, SolverType::FGMRES_GMG); + flow_problem.run(); + } + { + StokesProblem flow_problem(degree, SolverType::FGMRES_ILU); + flow_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/examples/step-56.output b/tests/examples/step-56.output new file mode 100644 index 0000000000..b4f4b729b1 --- /dev/null +++ b/tests/examples/step-56.output @@ -0,0 +1,66 @@ +Now running with Multigrid +Refinement cycle 0 + Set-up... + Number of active cells: 64 + Number of degrees of freedom: 2312 (2187+125) + Assembling... + Assembling Multigrid... + Solving... + Number of FGMRES iterations: 22 + Total number of iterations used for approximation of A inverse: 69 + Total number of iterations used for approximation of S inverse: 23 + + Note: The mean value was adjusted by -3.36485e-15 + + Velocity L2 Error: 0.00539774 + Pressure L2 Error: 0.019341 + Velocity H1 Error: 0.165265 +Refinement cycle 1 + Set-up... + Number of active cells: 512 + Number of degrees of freedom: 15468 (14739+729) + Assembling... + Assembling Multigrid... + Solving... + Number of FGMRES iterations: 21 + Total number of iterations used for approximation of A inverse: 67 + Total number of iterations used for approximation of S inverse: 22 + + Note: The mean value was adjusted by -1.51964e-13 + + Velocity L2 Error: 0.000670888 + Pressure L2 Error: 0.0036533 + Velocity H1 Error: 0.0414704 +Now running with ILU +Refinement cycle 0 + Set-up... + Number of active cells: 64 + Number of degrees of freedom: 2312 (2187+125) + Assembling... + Solving... Computing preconditioner... + + Number of FGMRES iterations: 22 + Total number of iterations used for approximation of A inverse: 113 + Total number of iterations used for approximation of S inverse: 23 + + Note: The mean value was adjusted by 1.55695e-14 + + Velocity L2 Error: 0.00539774 + Pressure L2 Error: 0.019341 + Velocity H1 Error: 0.165265 +Refinement cycle 1 + Set-up... + Number of active cells: 512 + Number of degrees of freedom: 15468 (14739+729) + Assembling... + Solving... Computing preconditioner... + + Number of FGMRES iterations: 21 + Total number of iterations used for approximation of A inverse: 180 + Total number of iterations used for approximation of S inverse: 22 + + Note: The mean value was adjusted by -1.52194e-13 + + Velocity L2 Error: 0.000670888 + Pressure L2 Error: 0.0036533 + Velocity H1 Error: 0.0414704