From 7907537c251b67e2296ee9ee543dacd834445cd3 Mon Sep 17 00:00:00 2001 From: Ryan Grove Date: Sat, 23 Jan 2016 23:22:37 -0500 Subject: [PATCH] Step-55: Geometric Multigrid for Stokes --- examples/step-55/CMakeLists.txt | 50 ++ examples/step-55/doc/builds-on | 1 + examples/step-55/doc/intro.dox | 71 ++ examples/step-55/doc/kind | 1 + examples/step-55/doc/results.dox | 162 ++++ examples/step-55/doc/tooltip | 1 + examples/step-55/step-55.cc | 1181 ++++++++++++++++++++++++++++++ 7 files changed, 1467 insertions(+) create mode 100644 examples/step-55/CMakeLists.txt create mode 100644 examples/step-55/doc/builds-on create mode 100644 examples/step-55/doc/intro.dox create mode 100644 examples/step-55/doc/kind create mode 100644 examples/step-55/doc/results.dox create mode 100644 examples/step-55/doc/tooltip create mode 100644 examples/step-55/step-55.cc diff --git a/examples/step-55/CMakeLists.txt b/examples/step-55/CMakeLists.txt new file mode 100644 index 0000000000..5a32dc83cd --- /dev/null +++ b/examples/step-55/CMakeLists.txt @@ -0,0 +1,50 @@ +## +# CMake script for the step-55 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "step-55") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) + +FIND_PACKAGE(deal.II 8.3 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +# +# Are all dependencies fulfilled? +# +IF(NOT DEAL_II_WITH_UMFPACK) + MESSAGE(FATAL_ERROR " +Error! The deal.II library found at ${DEAL_II_PATH} was not configured with + DEAL_II_WITH_UMFPACK = ON +One or all of these are OFF in your installation but are required for this tutorial step." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/examples/step-55/doc/builds-on b/examples/step-55/doc/builds-on new file mode 100644 index 0000000000..4d43c07055 --- /dev/null +++ b/examples/step-55/doc/builds-on @@ -0,0 +1 @@ +step-16,step-22 diff --git a/examples/step-55/doc/intro.dox b/examples/step-55/doc/intro.dox new file mode 100644 index 0000000000..f4923c9c77 --- /dev/null +++ b/examples/step-55/doc/intro.dox @@ -0,0 +1,71 @@ + +

Introduction

+ +

Stokes Problem

+ +The purpose of this tutorial is to create an efficient linear solver for the Stokes equation and compare it to alternative approaches. Using FGMRES with geometric multigrid as a precondtioner for the velocity block, we see that the linear solvers used in Step-22 cannot keep up since multigrid is the only way to get $O(n)$ solve time. Using the Timer class, we collect some statistics to compare setup times, solve times, and number of iterations. We also compute errors to make sure that what we have implemented is correct. + +Let $u \in H_0^1 = \{ u \in H^1(\Omega), u|_{\partial \Omega} = 0 \}$ and $p \in L_*^2 = \{ p_f \in L^2(\Omega), \int_\Omega p_f = 0 \}$. The Stokes equations read as follows in non-dimensionalized form: + +@f{eqnarray*} + - 2 \text{div} \frac {1}{2} \left[ (\nabla \textbf{u}) + (\nabla \textbf{u})^T\right] + \nabla p & =& f \\ + - \nabla \cdot u &=& 0 +@f} + +Note that we are using the deformation tensor instead of $\Delta u$ (a detailed desription of the difference between the two can be found in Step-22, but in summary, the deformation tensor is more physical as well as more expensive). + +

Linear solver and preconditioning issues

+The weak form of the discrete equations naturally leads to the following linear system for the nodal values of the velocity and pressure fields: +@f{eqnarray*} +\left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) \left(\begin{array}{c} U \\ P \end{array}\right) = \left(\begin{array}{c} F \\ 0 \end{array}\right), +@f} + +Our goal is to compare several solver approaches. In contrast to the way in which step-22 solves the Stokes equation, we instead attack the block system at once using a direct solver or FMGRES with an efficient preconditioner. The idea is as follows: if we find a block preconditioner $P$ such that the matrix +@f{eqnarray*} +\left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) P^{-1} +@f} + +is simple, then an iterative solver with that preconditioner will converge in a few iterations. Notice that we are doing right preconditioning for this. Using the Schur complement $S=BA^{-1}B^T$, we find that + +@f{eqnarray*} + P^{-1} = \left(\begin{array}{cc} \hat{A} & B^T \\ 0 & \hat{S} \end{array}\right)^{-1} +@f} + +is a good choice. It is important to note that +@f{eqnarray*} + P = \left(\begin{array}{cc} A^{-1} & 0 \\ 0 & I \end{array}\right) \left(\begin{array}{cc} I & B^T \\ 0 & -I \end{array}\right) \left(\begin{array}{cc} I & 0 \\ 0 & S^{-1} \end{array}\right). +@f} + +Since $P$ is aimed to be a preconditioner only, we shall use approximations to the inverse of the Schur complement $S$ and the matrix $A$. Therefore, in the above equations, $-M_p=\hat{S} \approx S$, where $M_p$ is the pressure mass matrix and is solved by using CG + ILU, and $\hat{A}$ is an approximation of $A$ obtained by one of multiple methods: CG + ILU, just using ILU, CG + GMG (Geometric Multigrid as described in Step-16), or just performing a few V-cycles of GMG. The inclusion of CG is more expensive, in general. + +As a comparison, instead of FGMRES, we also use the direct solver UMFPACK to compare our results to. If you want to use UMFPACK as a solver, it is important to note that since you have a singular system (since the integral of mean pressure being equal to zero not implemented), we set the first pressure node equal to zero since the direct solver can not handle the singular system like the other methods could. + +

Reference Solution

+ +The domain, right hand side, and boundary conditions we implemented were chosen for their simplicity and the fact that they made it possible for us to compute errors using a reference solution. We apply Dirichlet boundary condtions for the whole velocity on the whole boundary of the domain Ω=[0,1]×[0,1]×[0,1]. To enforce the boundary conditions we can just use our reference solution that we will now define. + +Let $u=(u_1,u_2,u_3)=(2\sin (\pi x), - \pi y \cos (\pi x),- \pi z \cos (\pi x))$ and $p = \sin (\pi x)\cos (\pi y)\sin (\pi z)$. + +If you look up in the deal.ii manual what is needed to create a class inherited from Function@, you will find not only a value function, but vector_value, value_list, etc. Different things you use in your code will require one of these particular functions. This can be confusing at first, but luckily the only thing you actually need to implement is value. The other ones have default implementations inside deal.ii and will be called on their own as long as you implement value correctly. + +Notice that our reference solution fulfills $\nabla \cdot u = 0$. In addition, the pressure is chosen to have a mean value of zero. For the Method of Manufactured Solutions of Step-7, we need to find $\bf f$ such that: + +@f{align*} +{\bf f} = - 2 \text{div} \frac {1}{2} \left[ (\nabla \textbf{u}) + (\nabla \textbf{u})^T\right] + \nabla p. +@f} + +Using the reference solution above, we obtain: + +@f{eqnarray*} +{\bf f} &=& (2 \pi^2 \sin (\pi x),- \pi^3 y \cos(\pi x),- \pi^3 z \cos(\pi x))\\ +& & + (\pi \cos(\pi x) \cos(\pi y) \sin(\pi z) ,- \pi \sin(\pi y) \sin(\pi x) \sin(\pi z), \pi \cos(\pi z) \sin(\pi x) \cos(\pi y)) +@f} + +

Computing Errors

+Because we do not enforce the mean pressure to be zero for our numerical solution in the linear system, we need to postprocess the solution after solving. To do this we use the compute_mean_value function to compute the mean value of the pressure to subtract it from the pressure. + +

DoF Handlers

+Geometric multigrid needs to know about the finite element system for the velocity. Since this is now part of the entire system, it is no longer easy to access. The reason for this is that there is currently no way in deal.ii to ask, "May I have just part of a DoF handler?" So in order to answer this request for our needs, we have to create a new DoF handler for just the velocites and assure that it has the same ordering as the DoF Handler for the entire system so that you can copy over one to the other. + +

Differences from Step-22

+The main difference between Step-55 and Step-22 is that we use block solvers instead of the Schur Complement approach used in step-22. Details of this approach can be found under the Block Schur complement preconditioner subsection of the Possible Extensions section of Step-22. For the preconditioner of the velocity block, we borrow a class from ASPECT called BlockSchurPreconditioner that has the option to solve for the inverse of $A$ or just apply one preconditioner sweep for it instead, which provides us with an expensive and cheap approach, respectively. \ No newline at end of file diff --git a/examples/step-55/doc/kind b/examples/step-55/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-55/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-55/doc/results.dox b/examples/step-55/doc/results.dox new file mode 100644 index 0000000000..af80c1a67a --- /dev/null +++ b/examples/step-55/doc/results.dox @@ -0,0 +1,162 @@ +

Results

+ +

Errors

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
 L2 VelocityRateL2 PressureRateH1 VelocityRate
3D, 3 global refinements0.000670888-0.0036533-0.0414704-
3D, 4 global refinements8.38E-0058.00732356780.000884944.12830248380.01037813.9959530164
3D, 5 global refinements1.05E-0058.00086899230.0002202534.01783403630.002595193.9989750269
+ +

Timing Results

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
   UMFPACKILUGMG
    Timings TimingsIterations TimingsIterations
 Set-upAssembleVM PeakSetupSolveVM PeakSetupSolveOuterInner (A)Inner (S)VM PeakSetupSolveOuterInner (A)Inner (S)
3D, 3 global refinements  53059082.16s2.2s50133080.37s0.305s271072850360920.37s1.04s234724
3D, 4 global refinements  11674656231s225s53656763.15s8.21s291863056899522.72s9.14s234824
3D, 5 global refinements  ---863676029.4s152s32361341102529621s78.1s224623
+ +As can be seen from the table, + +1. UMFPACK uses large amounts of memory, especially in 3d. Also, UMFPACK timings do not scale with problem size. + +2. The number of iterations for $A$ increase for ILU with refinement leading to worse then linear scaling in solve time. In contrast, the number of inner iterations for $A$ stay constant with GMG leading to nearly perfect scaling in solve time. + +3. GMG needs slightly more memory than ILU. + +

Possible extensions

+ +

Extension 1

+ diff --git a/examples/step-55/doc/tooltip b/examples/step-55/doc/tooltip new file mode 100644 index 0000000000..2c54b81f63 --- /dev/null +++ b/examples/step-55/doc/tooltip @@ -0,0 +1 @@ +Precondtioners: Geometric multigrid vs ILU diff --git a/examples/step-55/step-55.cc b/examples/step-55/step-55.cc new file mode 100644 index 0000000000..e2c021e839 --- /dev/null +++ b/examples/step-55/step-55.cc @@ -0,0 +1,1181 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 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. + * + * --------------------------------------------------------------------- + + * 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 +#include + +// We need to include this class for all of the timings between ILU and Multigrid +#include + +// This includes the files necessary for us to use geometric Multigrid +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace Step55 +{ + using namespace dealii; + + // In order to make it easy to switch between the different solvers that are being + // used in Step-55, an enum was created 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. We do this to make it easier for us to debug + // as well as its aesthetic value. + template + class Solution : public Function + { + public: + Solution () : Function(dim+1) {} + virtual double value (const Point &p, + const unsigned int component) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; + }; + + template <> + double + Solution<2>::value (const Point<2> &p, + const unsigned int component) const + { + using numbers::PI; + double x = p(0); + 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); + + Assert (false, ExcMessage ("Component out of range in Solution")); + return 0; + } + + template <> + double + Solution<3>::value (const Point<3> &p, + const unsigned int component) const + { + using numbers::PI; + double x = p(0); + double y = p(1); + 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); + + Assert (false, ExcMessage ("Component out of range in Solution")); + 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 + { + using numbers::PI; + double x = p(0); + 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); + } + else + Assert (false, ExcMessage ("Component out of range in Solution")); + return return_value; + } + + template <> + Tensor<1,3> + Solution<3>::gradient (const Point<3> &p, + const unsigned int component) const + { + using numbers::PI; + double x = p(0); + double y = p(1); + 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); + } + else + Assert (false, ExcMessage ("Component out of range in Solution")); + 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; + }; + + template <> + double + RightHandSide<2>::value (const Point<2> &p, + const unsigned int component) const + { + 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; + + Assert (false, ExcMessage ("Component out of range in RightHandSide")); + return 0; + + } + + template <> + double + RightHandSide<3>::value (const Point<3> &p, + const unsigned int component) const + { + 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; + + Assert (false, ExcMessage ("Component out of range in RightHandSide")); + return 0; + } + + // Sadly, we need a separate function for the boundary conditions + // to be used in the geometric multigrid. This is because it needs + // to be a function with $dim$ components, whereas Solution has + // $dim+1$ components. Rather than copying the implementation of + // Solution, we forward the calls to Solution::value. For that we need + // an instance of the class Solution, which you can find as a private + // member. + template + class BoundaryValuesForVelocity : public Function + { + public: + BoundaryValuesForVelocity () : Function(dim) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + private: + Solution solution; + }; + + + template + double + BoundaryValuesForVelocity::value (const Point &p, + const unsigned int component) const + { + Assert (component < this->n_components, + ExcIndexRange (component, 0, this->n_components)); + + return solution.value(p, component); + } + + + + // @sect3{ASPECT BlockSchurPreconditioner} + + // This class, which is taken from ASPECT and then slightly modified, + // implements the block Schur preconditioner for the Stokes system discussed + // above. It is templated on the types + // for the preconditioner blocks for velocity and schur complement. + // + // 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 (if do_solve_A==true) + // or just apply one preconditioner sweep + if (do_solve_A == true) + { + SolverControl solver_control(10000, utmp.l2_norm()*1e-2, true); + 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 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 degree; + SolverType::type solver_type; + + Triangulation triangulation; + FESystem fe; + FESystem velocity_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 degree, SolverType::type solver_type) + : + degree (degree), + solver_type (solver_type), + triangulation (Triangulation::maximum_smoothing), + fe (FE_Q(degree+1), dim, // Finite element for whole system + FE_Q(degree), 1), + velocity_fe (FE_Q(degree+1), dim), // Finite element for velocity-only + dof_handler (triangulation), + velocity_dof_handler (triangulation), + computing_timer (std::cout, TimerOutput::summary, + TimerOutput::wall_times) + {} + + + +// @sect4{StokesProblem::setup_dofs} + +// This function sets up things differently based on if you want to use ILU or GMG as a preconditioner. + template + void StokesProblem::setup_dofs () + { + TimerOutput::Scope scope(computing_timer, "Setup"); + + system_matrix.clear (); + pressure_mass_matrix.clear (); + + // We don't need the multigrid dofs for whole problem finite element + dof_handler.distribute_dofs(fe); + + // This first creates and array (0,0,1) which means that it first does everything with index 0 and then 1 + std::vector block_component (dim+1,0); + block_component[dim] = 1; + + // This always knows how to use the dim (start at 0 one) + FEValuesExtractors::Vector velocities(0); + + if (solver_type == SolverType::FGMRES_ILU) + { + TimerOutput::Scope ilu_specific(computing_timer, "(ILU specific)"); + DoFRenumbering::Cuthill_McKee (dof_handler); + } + + DoFRenumbering::component_wise (dof_handler, block_component); + + if (solver_type == SolverType::FGMRES_GMG) + { + TimerOutput::Scope multigrid_specific(computing_timer, "(Multigrid specific)"); + TimerOutput::Scope setup_multigrid(computing_timer, "Setup - Multigrid"); + + // Distribute only the dofs for velocity finite element + velocity_dof_handler.distribute_dofs(velocity_fe); + + // Multigrid only needs the dofs for velocity + velocity_dof_handler.distribute_mg_dofs(velocity_fe); + + typename FunctionMap::type boundary_condition_function_map; + BoundaryValuesForVelocity velocity_boundary_condition; + boundary_condition_function_map[0] = &velocity_boundary_condition; + + mg_constrained_dofs.clear(); + mg_constrained_dofs.initialize(velocity_dof_handler, boundary_condition_function_map); + const unsigned int n_levels = triangulation.n_levels(); + + mg_interface_matrices.resize(0, n_levels-1); + mg_interface_matrices.clear (); + mg_matrices.resize(0, n_levels-1); + mg_matrices.clear (); + mg_sparsity_patterns.resize(0, n_levels-1); + + for (unsigned int level=0; level(), + constraints, + fe.component_mask(velocities)); + } + + 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]; + + if (solver_type == SolverType::UMFPACK) + { + TimerOutput::Scope umfpack_specific(computing_timer, "(UMFPACK specific)"); + 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 (2,2); + + csp.block(0,0).reinit (n_u, n_u); + csp.block(1,0).reinit (n_p, n_u); + csp.block(0,1).reinit (n_u, n_p); + csp.block(1,1).reinit (n_p, n_p); + + csp.collect_sizes(); + + DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false); + sparsity_pattern.copy_from (csp); + + } + system_matrix.reinit (sparsity_pattern); + + solution.reinit (2); + solution.block(0).reinit (n_u); + solution.block(1).reinit (n_p); + solution.collect_sizes (); + + system_rhs.reinit (2); + system_rhs.block(0).reinit (n_u); + system_rhs.block(1).reinit (n_p); + system_rhs.collect_sizes (); + } + + +// @sect4{StokesProblem::assemble_system} + +// In this function, the system matrix is assembled the same regardless of using ILU and GMG + template + void StokesProblem::assemble_system () + { + TimerOutput::Scope assemble(computing_timer, "Assemble"); + system_matrix=0; + system_rhs=0; + + double mass_factor = (solver_type == SolverType::UMFPACK) ? 0.0 : 1.0; + + 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); + 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; qget_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 step-40, we have a function that assembles everything 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(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 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); + + std::vector boundary_constraints (triangulation.n_levels()); + std::vector boundary_interface_constraints (triangulation.n_levels()); + for (unsigned int level=0; level::cell_iterator cell = velocity_dof_handler.begin(), + endc = velocity_dof_handler.end(); + + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + cell_matrix = 0; + + right_hand_side.vector_value_list(fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q=0; qget_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; ilevel(), + 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 (GMRES) but require a different preconditioner to be assembled. Here we time not only the entire solve +// function, but we separately time the set-up of the preconditioner as well as the GMRES solve. + 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 + 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 pmass_preconditioner; + pmass_preconditioner.initialize (pressure_mass_matrix); + + const BlockSchurPreconditioner, SparseILU > + preconditioner (system_matrix, + pressure_mass_matrix, + A_preconditioner, + pmass_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(constraints, mg_constrained_dofs); + mg_transfer.build_matrices(velocity_dof_handler); + + // Coarse grid solver + // Timo: TODO: should we use something like LACIteration? + FullMatrix coarse_matrix; + coarse_matrix.copy_from (mg_matrices[0]); + MGCoarseGridHouseholder<> coarse_grid_solver; + coarse_grid_solver.initialize (coarse_matrix); + + typedef PreconditionSOR > Smoother; + mg::SmootherRelaxation > mg_smoother; + mg_smoother.initialize(mg_matrices); + mg_smoother.set_steps(2); + + // Multigrid, when used as a preconditioner for CG, expects the + // smoother to be symmetric, and this takes care of that + 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(velocity_dof_handler, + 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 pmass_preconditioner; + pmass_preconditioner.initialize (pressure_mass_matrix, + SparseILU::AdditionalData()); + + const BlockSchurPreconditioner< + PreconditionMG, MGTransferPrebuilt > >, + SparseILU > + preconditioner (system_matrix, + pressure_mass_matrix, + A_Multigrid, + pmass_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 iterations used for block GMRES iterations: " + << solver_control.last_step() << std::endl + << "\tNumber of iterations used for approximation of A inverse: " + << n_iterations_A << std::endl + << "\tNumber of iterations used for approximation of S inverse: " + << n_iterations_S << std::endl + << std::endl; + } + + +// @sect4{StokesProblem::process_solution} + + 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. + double mean_pressure = VectorTools::compute_mean_value (dof_handler, + QGauss(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); + // Timo: TODO: find a better way to do this inside deal.II, maybe with component_mask + const ComponentSelectFunction velocity_mask(std::make_pair(0, dim), dim+1); + + /* + Timo: TODO: think about this + Extractor::Vector velocities(0); + const ComponentSelectFunction + velocity_selector(velocities, dim+1); + */ + + Vector difference_per_cell (triangulation.n_active_cells()); + + VectorTools::integrate_difference (dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss(degree+2), + VectorTools::L2_norm, + &velocity_mask); + + const double Velocity_L2_error = difference_per_cell.l2_norm(); + + VectorTools::integrate_difference (dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss(degree+2), + VectorTools::L2_norm, + &pressure_mask); + + const double Pressure_L2_error = difference_per_cell.l2_norm(); + + VectorTools::integrate_difference (dof_handler, + solution, + Solution(), + difference_per_cell, + QGauss(degree+2), + VectorTools::H1_norm, + &velocity_mask); + + const double Velocity_H1_error = difference_per_cell.l2_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} + + template + void + StokesProblem::output_results (const unsigned int refinement_cycle) const + { + std::vector solution_names (dim, "velocity"); + solution_names.push_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::ostringstream filename; + filename << "solution-" + << Utilities::int_to_string (refinement_cycle, 2) + << ".vtk"; + + std::ofstream output (filename.str().c_str()); + 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 (6-dim); + + 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<3; + ++refinement_cycle) + { + std::cout << "Refinement cycle " << refinement_cycle << std::endl; + + if (refinement_cycle > 0) + triangulation.refine_global (1); + + std::cout << " Set-up..." << std::endl << std::flush; + setup_dofs(); + + std::cout << " Assembling..." << std::endl << std::flush; + assemble_system (); + + if (solver_type == SolverType::FGMRES_GMG) + { + std::cout << " Assembling Multigrid..." << std::endl << std::flush; + + assemble_multigrid (); + } + + std::cout << " Solving..." << std::flush; + solve (); + + compute_errors (); + + Utilities::System::MemoryStats mem; + Utilities::System::get_memory_stats(mem); + std::cout << " VM Peak: " << mem.VmPeak << std::endl; + + computing_timer.print_summary (); + computing_timer.reset (); + output_results (refinement_cycle); + } + } +} + +// @sect3{The main function} +int main () +{ + try + { + using namespace dealii; + using namespace Step55; + + 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 (); + } + 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; +} + -- 2.39.5