From: Jie Cheng Date: Mon, 13 Nov 2017 05:57:57 +0000 (-0500) Subject: Time-dependent Navier-Stokes X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f63511691ac8977087d52144b1508f83ffee5642;p=code-gallery.git Time-dependent Navier-Stokes --- diff --git a/time_dependent_navier_stokes/CMakeLists.txt b/time_dependent_navier_stokes/CMakeLists.txt new file mode 100644 index 0000000..8eac2e3 --- /dev/null +++ b/time_dependent_navier_stokes/CMakeLists.txt @@ -0,0 +1,53 @@ +## +# CMake script for the time-dependent Navier-Stokes solver: +## + +# Set the name of the project and target: +SET(TARGET "fluid") + +# 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 + time_dependent_navier_stokes.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) + +FIND_PACKAGE(deal.II 9.0.0 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) # keep in one line + MESSAGE(FATAL_ERROR " +Error! This tutorial requires a deal.II library that was configured with the following options: + DEAL_II_WITH_UMFPACK = ON +However, the deal.II library found at ${DEAL_II_PATH} was configured with these options + DEAL_II_WITH_UMFPACK = ${DEAL_II_WITH_UMFPACK} +which conflict with the requirements." + ) +ENDIF() + + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/time_dependent_navier_stokes/Readme.md b/time_dependent_navier_stokes/Readme.md new file mode 100644 index 0000000..989029b --- /dev/null +++ b/time_dependent_navier_stokes/Readme.md @@ -0,0 +1,138 @@ +Time-dependent Navier-Stokes +------------------------------------------ + +### General description of the problem ### + +We solve the time-dependent Navier-Stokes equations with implicit-explicit (IMEX) scheme. +Here is the equations we want to solve: +@f{eqnarray*} +{\bold{u}}_{,t} - \nu {\nabla}^2\bold{u} + (\bold{u}\cdot\nabla)\bold{u} + \nabla p = \bold{f} +@f} + +The idea is as follows: we use backward Euler time for time discretization. The diffusion term +is treated implicitly and the convection term is treated explicitly. Let $(u, p)$ denote the +velocity and pressure, respectively and $(v, q)$ denote the corresponding test functions, we +end up with the following linear system: +@f{eqnarray*} +m(u^{n+1}, v) + \Delta{t}\cdot a((u^{n+1}, p^{n+1}), (v, q))=m(u^n, v)-\Delta{t}c(u^n;u^n, v) +@f} + +where $a((u, p), (v, q))$ is the bilinear form of the diffusion term: +@f{eqnarray*} +a((u, p), (v, q)) = \int_\Omega \nu\nabla{u}\nabla{v}-p\nabla\cdot v-q\nabla\cdot ud\Omega +@f} + +$m(u, v)$ is the mass matrix: +@f{eqnarray*} +m(u, v) = \int_{\Omega} u \cdot v d\Omega +@f} + +and $c(u;u, v)$ is the convection term: +@f{eqnarray*} +c(u;u, v) = \int_{\Omega} (u \cdot \nabla u) \cdot v d\Omega +@f} + +Substracting $m(u^n, v) + \Delta{t}a((u^n, p^n), (v, q))$ from both sides of the equation, +we have the incremental form: +@f{eqnarray*} +m(\Delta{u}, v) + \Delta{t}\cdot a((\Delta{u}, \Delta{p}), (v, q)) = \Delta{t}(-a(u^n, p^n), (q, v)) - \Delta{t}c(u^n;u^n, v) +@f} + + +The system we want to solve can be written in matrix form: + +@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} + +#### Grad-Div stablization #### + +Similar to step-57, we add $\gamma B^T M_p^{-1} B$ to the upper left block of the system, +thus the system becomes: + +@f{eqnarray*} + \left( + \begin{array}{cc} + \tilde{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} +where $\tilde{A} = A + \gamma B^T M_p^{-1} B$. + +Detailed explaination of the Grad-Div stablization can be found in [1]. + +#### Block preconditioner #### + +The block preconditioner is pretty much the same as in step-22, except for two additional terms, +namely the inertial term (mass matrix) and the Grad-Div term. + +The block preconditioner can be written as: +@f{eqnarray*} + P^{-1} + = + \left( + \begin{array}{cc} + {\tilde{A}}^{-1} & 0 \\ + {\tilde{S}}^{-1}B{\tilde{A}}^{-1} & -{\tilde{S}}^{-1} \\ + \end{array} + \right) +@f} +where ${\tilde{S}}$ is the Schur complement of ${\tilde{A}}$, which can be decomposed +into the Schur complements of the mass matrix, diffusion matrix, and the Grad-Div term: +@f{eqnarray*} + {\tilde{S}}^{-1} + \approx + {S_{mass}}^{-1} + {S_{diff}}^{-1} + {S_{Grad-Div}}^{-1} + \approx + [B^T (diag M)^{-1} B]^{-1} + \Delta{t}(\nu + \gamma)M_p^{-1} +@f} + +For more information about preconditioning incompressible Navier-Stokes equations, please refer +to [1] and [2]. + +#### Test case #### +We test the code with a classical benchmark case, flow past a cylinder, with Reynold's +number 100. The geometry setup of the case can be found on +[this webpage](http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html). + +### Acknowledgements ### +Thank Wolfgang Bangerth, Timo Heister and Martin Kronbichler for their helpful discussions +on my numerical formulation and implementation. + +------------------------------------------ +### References ### +[1] Timo Heister. A massively parallel finite element framework with application to incompressible flows. Doctoral dissertation, University of Gottingen, 2011. + +[2] M. Kronbichler, A. Diagne and H. Holmgren. A fast massively parallel two-phase flow solver for microfluidic chip simulation, International Journal of High Performance Computing Applications, 2016. diff --git a/time_dependent_navier_stokes/doc/Re100.mp4 b/time_dependent_navier_stokes/doc/Re100.mp4 new file mode 100644 index 0000000..b050e0f Binary files /dev/null and b/time_dependent_navier_stokes/doc/Re100.mp4 differ diff --git a/time_dependent_navier_stokes/doc/Re100.png b/time_dependent_navier_stokes/doc/Re100.png new file mode 100644 index 0000000..77027e3 Binary files /dev/null and b/time_dependent_navier_stokes/doc/Re100.png differ diff --git a/time_dependent_navier_stokes/doc/author b/time_dependent_navier_stokes/doc/author new file mode 100644 index 0000000..2995ade --- /dev/null +++ b/time_dependent_navier_stokes/doc/author @@ -0,0 +1 @@ +Jie Cheng diff --git a/time_dependent_navier_stokes/doc/builds-on b/time_dependent_navier_stokes/doc/builds-on new file mode 100644 index 0000000..9818e36 --- /dev/null +++ b/time_dependent_navier_stokes/doc/builds-on @@ -0,0 +1 @@ +step-20 step-22 step-57 diff --git a/time_dependent_navier_stokes/doc/dependencies b/time_dependent_navier_stokes/doc/dependencies new file mode 100644 index 0000000..3f4cba3 --- /dev/null +++ b/time_dependent_navier_stokes/doc/dependencies @@ -0,0 +1 @@ +DEAL_II_WITH_CXX11 diff --git a/time_dependent_navier_stokes/doc/entry-name b/time_dependent_navier_stokes/doc/entry-name new file mode 100644 index 0000000..03d7b99 --- /dev/null +++ b/time_dependent_navier_stokes/doc/entry-name @@ -0,0 +1 @@ +Time-dependent Navier-Stokes diff --git a/time_dependent_navier_stokes/doc/tooltip b/time_dependent_navier_stokes/doc/tooltip new file mode 100644 index 0000000..ca55f62 --- /dev/null +++ b/time_dependent_navier_stokes/doc/tooltip @@ -0,0 +1 @@ +Solving time-dependent incompressible Navier-Stokes problem with Grad-Div stablization using IMEX scheme. diff --git a/time_dependent_navier_stokes/time_dependent_navier_stokes.cc b/time_dependent_navier_stokes/time_dependent_navier_stokes.cc new file mode 100644 index 0000000..85f0976 --- /dev/null +++ b/time_dependent_navier_stokes/time_dependent_navier_stokes.cc @@ -0,0 +1,1008 @@ +#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 +#include +#include + +#include + +#include +#include +#include + +namespace fluid +{ + using namespace dealii; + + // @sect3{Create the triangulation} + + // The code to create triangulation is copied from Martin Kronbichler's code + // (https://github.com/kronbichler/adaflo/blob/master/tests/flow_past_cylinder.cc) + // with very few modifications. + // Helper function used in both 2d and 3d: + void create_triangulation_2d(Triangulation<2> &tria, bool compute_in_2d = true) + { + SphericalManifold<2> boundary(Point<2>(0.5, 0.2)); + Triangulation<2> left, middle, right, tmp, tmp2; + GridGenerator::subdivided_hyper_rectangle( + left, + std::vector({3U, 4U}), + Point<2>(), + Point<2>(0.3, 0.41), + false); + GridGenerator::subdivided_hyper_rectangle( + right, + std::vector({18U, 4U}), + Point<2>(0.7, 0), + Point<2>(2.5, 0.41), + false); + + // Create middle part first as a hyper shell. + GridGenerator::hyper_shell(middle, Point<2>(0.5, 0.2), 0.05, 0.2, 4, true); + middle.set_manifold(0, boundary); + middle.refine_global(1); + + // Then move the vertices to the points where we want them to be to create a + // slightly asymmetric cube with a hole + for (Triangulation<2>::cell_iterator cell = middle.begin(); + cell != middle.end(); + ++cell) + for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell; ++v) + { + Point<2> &vertex = cell->vertex(v); + if (std::abs(vertex[0] - 0.7) < 1e-10 && + std::abs(vertex[1] - 0.2) < 1e-10) + vertex = Point<2>(0.7, 0.205); + else if (std::abs(vertex[0] - 0.6) < 1e-10 && + std::abs(vertex[1] - 0.3) < 1e-10) + vertex = Point<2>(0.7, 0.41); + else if (std::abs(vertex[0] - 0.6) < 1e-10 && + std::abs(vertex[1] - 0.1) < 1e-10) + vertex = Point<2>(0.7, 0); + else if (std::abs(vertex[0] - 0.5) < 1e-10 && + std::abs(vertex[1] - 0.4) < 1e-10) + vertex = Point<2>(0.5, 0.41); + else if (std::abs(vertex[0] - 0.5) < 1e-10 && + std::abs(vertex[1] - 0.0) < 1e-10) + vertex = Point<2>(0.5, 0.0); + else if (std::abs(vertex[0] - 0.4) < 1e-10 && + std::abs(vertex[1] - 0.3) < 1e-10) + vertex = Point<2>(0.3, 0.41); + else if (std::abs(vertex[0] - 0.4) < 1e-10 && + std::abs(vertex[1] - 0.1) < 1e-10) + vertex = Point<2>(0.3, 0); + else if (std::abs(vertex[0] - 0.3) < 1e-10 && + std::abs(vertex[1] - 0.2) < 1e-10) + vertex = Point<2>(0.3, 0.205); + else if (std::abs(vertex[0] - 0.56379) < 1e-4 && + std::abs(vertex[1] - 0.13621) < 1e-4) + vertex = Point<2>(0.59, 0.11); + else if (std::abs(vertex[0] - 0.56379) < 1e-4 && + std::abs(vertex[1] - 0.26379) < 1e-4) + vertex = Point<2>(0.59, 0.29); + else if (std::abs(vertex[0] - 0.43621) < 1e-4 && + std::abs(vertex[1] - 0.13621) < 1e-4) + vertex = Point<2>(0.41, 0.11); + else if (std::abs(vertex[0] - 0.43621) < 1e-4 && + std::abs(vertex[1] - 0.26379) < 1e-4) + vertex = Point<2>(0.41, 0.29); + } + + // Refine once to create the same level of refinement as in the + // neighboring domains: + middle.refine_global(1); + + // Must copy the triangulation because we cannot merge triangulations with + // refinement: + GridGenerator::flatten_triangulation(middle, tmp2); + + // Left domain is requred in 3d only. + if (compute_in_2d) + { + GridGenerator::merge_triangulations(tmp2, right, tria); + } + else + { + GridGenerator::merge_triangulations(left, tmp2, tmp); + GridGenerator::merge_triangulations(tmp, right, tria); + } + } + + // Create 2D triangulation: + void create_triangulation(Triangulation<2> &tria) + { + create_triangulation_2d(tria); + // Set the cylinder boundary to 1, the right boundary (outflow) to 2, the rest to 0. + for (Triangulation<2>::active_cell_iterator cell = tria.begin(); + cell != tria.end(); + ++cell) + { + for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f) + { + if (cell->face(f)->at_boundary()) + { + if (std::abs(cell->face(f)->center()[0] - 2.5) < 1e-12) + { + cell->face(f)->set_all_boundary_ids(2); + } + else if (Point<2>(0.5, 0.2).distance(cell->face(f)->center()) <= 0.05) + { + cell->face(f)->set_all_manifold_ids(10); + cell->face(f)->set_all_boundary_ids(1); + } + else + { + cell->face(f)->set_all_boundary_ids(0); + } + } + } + } + } + + // Create 3D triangulation: + void create_triangulation(Triangulation<3> &tria) + { + Triangulation<2> tria_2d; + create_triangulation_2d(tria_2d, false); + GridGenerator::extrude_triangulation(tria_2d, 5, 0.41, tria); + // Set the cylinder boundary to 1, the right boundary (outflow) to 2, the rest to 0. + for (Triangulation<3>::active_cell_iterator cell = tria.begin(); + cell != tria.end(); ++cell) + { + for (unsigned int f = 0; f::faces_per_cell; ++f) + { + if (cell->face(f)->at_boundary()) + { + if (std::abs(cell->face(f)->center()[0]-2.5) < 1e-12) + { + cell->face(f)->set_all_boundary_ids(2); + } + else if (Point<3>(0.5, 0.2, cell->face(f)->center()[2]).distance + (cell->face(f)->center()) <= 0.05) + { + cell->face(f)->set_all_manifold_ids(10); + cell->face(f)->set_all_boundary_ids(1); + } + else + { + cell->face(f)->set_all_boundary_ids(0); + } + } + } + } + } + + // @sect3{Time stepping} + class Time + { + public: + Time(const double time_end, const double delta_t) + : timestep(0), time_current(0.0), time_end(time_end), delta_t(delta_t) + { + } + virtual ~Time() {} + double current() const { return time_current; } + double end() const { return time_end; } + double get_delta_t() const { return delta_t; } + unsigned int get_timestep() const { return timestep; } + void increment() + { + time_current += delta_t; + ++timestep; + } + + private: + unsigned int timestep; + double time_current; + const double time_end; + const double delta_t; + }; + + // @sect3{Boundary values} + + // Dirichlet boundary conditions for the velocity inlet and walls + template + class BoundaryValues : public Function + { + public: + BoundaryValues() : Function(dim + 1) {} + virtual double value(const Point &p, + const unsigned int component) const; + + virtual void vector_value(const Point &p, + Vector &values) const; + }; + + template + double BoundaryValues::value(const Point &p, + const unsigned int component) const + { + Assert(component < this->n_components, + ExcIndexRange(component, 0, this->n_components)); + if (component == 0 && std::abs(p[0] - 0.3) < 1e-10) + { + double U = 1.5; + double y = p[1]; + return 4 * U * y * (0.41 - y) / (0.41 * 0.41); + } + return 0; + } + + template + void BoundaryValues::vector_value(const Point &p, + Vector &values) const + { + for (unsigned int c = 0; c < this->n_components; ++c) + values(c) = BoundaryValues::value(p, c); + } + + // @sect3{Preconditioners} + + // The LHS of the system matrix is the same as Stokes equation for IMEX scheme. + // A block preconditioner as in step-22 is used here. + + // @sect4{Inner preconditioner} + + // Adapted from step-22, used to solve for ${\tilde{A}}^{-1}$ + template + struct InnerPreconditioner; + + template <> + struct InnerPreconditioner<2> + { + typedef SparseDirectUMFPACK type; + }; + + template <> + struct InnerPreconditioner<3> + { + typedef SparseILU type; + }; + + // @sect4{Inverse matrix} + + // This is used for ${\tilde{S}}^{-1}$ and ${\tilde{A}}^{-1}$, which are symmetric so we use CG + // solver inside + template + class InverseMatrix : public Subscriptor + { + public: + InverseMatrix(const MatrixType &m, + const PreconditionerType &preconditioner); + void vmult(Vector &dst, const Vector &src) const; + + private: + const SmartPointer matrix; + const SmartPointer preconditioner; + }; + + template + InverseMatrix::InverseMatrix( + const MatrixType &m, const PreconditionerType &preconditioner) + : matrix(&m), preconditioner(&preconditioner) + { + } + + template + void InverseMatrix::vmult( + Vector &dst, const Vector &src) const + { + SolverControl solver_control(src.size(), 1e-6 * src.l2_norm()); + SolverCG<> cg(solver_control); + dst = 0; + cg.solve(*matrix, dst, src, *preconditioner); + } + + // @sect4{Approximate Schur complement of mass matrix} + + // The Schur complement of mass matrix is written as $S_M = BM^{-1}B^T$ + // Similar to step-20, we use $B(diag(M))^{-1}B^T$ to approximate it. + class ApproximateMassSchur : public Subscriptor + { + public: + ApproximateMassSchur(const BlockSparseMatrix &M); + void vmult(Vector &dst, const Vector &src) const; + + private: + const SmartPointer> mass_matrix; + mutable Vector tmp1, tmp2; + }; + + ApproximateMassSchur::ApproximateMassSchur( + const BlockSparseMatrix &M) + : mass_matrix(&M), tmp1(M.block(0, 0).m()), tmp2(M.block(0, 0).m()) + { + } + + void ApproximateMassSchur::vmult(Vector &dst, + const Vector &src) const + { + mass_matrix->block(0, 1).vmult(tmp1, src); + mass_matrix->block(0, 0).precondition_Jacobi(tmp2, tmp1); + mass_matrix->block(1, 0).vmult(dst, tmp2); + } + + // @sect4{The inverse matrix of the system Schur complement} + + // The inverse of the total Schur complement is the sum of the inverse of + // diffusion, Grad-Div term, and mass Schur complements. Note that the first + // two components add up to $\Delta{t}(\nu + \gamma)M_p^{-1}$ as introduced in step-57, + // in which the additional $\Delta{t}$ comes from the time discretization, + // and the last component is obtained by wrapping a InverseMatrix<\code> + // around ApproximateMassSchur<\code>. + template + class SchurComplementInverse : public Subscriptor + { + public: + SchurComplementInverse( + double gamma, double viscosity, double dt, + const InverseMatrix &Sm_inv, + const InverseMatrix, PreconditionerMp> &Mp_inv); + void vmult(Vector &dst, const Vector &src) const; + private: + const double gamma; + const double viscosity; + const double dt; + const SmartPointer> Sm_inverse; + const SmartPointer, + PreconditionerMp>> Mp_inverse; + }; + + template + SchurComplementInverse::SchurComplementInverse( + double gamma, double viscosity, double dt, + const InverseMatrix &Sm_inv, + const InverseMatrix, PreconditionerMp> &Mp_inv) : + gamma(gamma), viscosity(viscosity), dt(dt), Sm_inverse(&Sm_inv), Mp_inverse(&Mp_inv) + { + } + + template + void SchurComplementInverse::vmult( + Vector &dst, const Vector &src) const + { + Vector tmp(src.size()); + Sm_inverse->vmult(dst, src); + Mp_inverse->vmult(tmp, src); + tmp *= (viscosity + gamma) * dt; + dst += tmp; + } + + // @sect4{The block Schur preconditioner} + + // The block Schur preconditioner has the same form as in step-22, which is written as + // $P^{-1} = [\tilde{A}}^{-1}, 0; {\tilde{S}}^{-1}B{\tilde{A}}^{-1}, -{\tilde{S}}^{-1}]$ + // Note that ${\tilde{A}}^{-1}$ has contributions from the diffusion, Grad-Div and mass terms. + // This class has three template arguments: PreconditionerA is needed for ${\tilde{A}}^{-1}$, + // PreconditionerSm and PreconditionerMp are used in the inverse of the Schur complement + // of $\tilde{A}$, namely ${\tilde{S}}^{-1}$. + template + class BlockSchurPreconditioner : public Subscriptor + { + public: + BlockSchurPreconditioner( + const BlockSparseMatrix &system_m, + const InverseMatrix, PreconditionerA> &A_inv, + const SchurComplementInverse &S_inv); + void vmult(BlockVector &dst, const BlockVector &src) const; + + private: + const SmartPointer> system_matrix; + const SmartPointer< + const InverseMatrix, PreconditionerA>> A_inverse; + const SmartPointer< + const SchurComplementInverse> S_inverse; + mutable Vector tmp; + }; + + template + BlockSchurPreconditioner:: + BlockSchurPreconditioner( + const BlockSparseMatrix &system_m, + const InverseMatrix, PreconditionerA> &A_inv, + const SchurComplementInverse &S_inv) + : system_matrix(&system_m), A_inverse(&A_inv), S_inverse(&S_inv), + tmp(system_matrix->block(1, 1).m()) + { + } + + template + void BlockSchurPreconditioner::vmult( + BlockVector &dst, const BlockVector &src) const + { + A_inverse->vmult(dst.block(0), src.block(0)); + system_matrix->block(1, 0).residual(tmp, dst.block(0), src.block(1)); + tmp *= -1; + S_inverse->vmult(dst.block(1), tmp); + } + + // @sect3{The time-dependent Navier-Stokes class template} + template + class NavierStokes + { + public: + NavierStokes(const unsigned int degree); + void run(); + + private: + void setup(); + void assemble(bool assemble_lhs); + + std::pair solve_linear_system(bool update_preconditioner); + void output_results(const unsigned int index) const; + void process_solution(std::ofstream& out) const; + const ConstraintMatrix &get_constraints() const; + + double viscosity; + double gamma; + const unsigned int degree; + std::vector dofs_per_block; + + Triangulation triangulation; + FESystem fe; + DoFHandler dof_handler; + QGauss quadrature_formula; + QGauss face_quadrature_formula; + + ConstraintMatrix zero_constraints; + ConstraintMatrix nonzero_constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + // We need both velocity mass and pressure mass, so we use a block sparse matrix to store it. + BlockSparseMatrix mass_matrix; + + BlockVector solution; + BlockVector solution_increment; + BlockVector system_rhs; + + Time time; + mutable TimerOutput timer; + + // We use shared pointers for all the preconditioning-related stuff + std::shared_ptr approximate_Sm; + std::shared_ptr preconditioner_Sm; + std::shared_ptr> Sm_inverse; + + std::shared_ptr> preconditioner_Mp; + std::shared_ptr, SparseILU>> Mp_inverse; + + std::shared_ptr>> S_inverse; + + std::shared_ptr::type> preconditioner_A; + std::shared_ptr, + typename InnerPreconditioner::type>> A_inverse; + + std::shared_ptr::type, PreconditionIdentity, SparseILU>> preconditioner; + }; + + // @sect4{NavierStokes::NavierStokes} + template + NavierStokes::NavierStokes(const unsigned int degree) + : viscosity(0.001), + gamma(1), + degree(degree), + triangulation(Triangulation::maximum_smoothing), + fe(FE_Q(degree + 1), dim, FE_Q(degree), 1), + dof_handler(triangulation), + quadrature_formula(degree+2), + face_quadrature_formula(degree+2), + time(1e-2, 1e-3), + timer(std::cout, TimerOutput::summary, TimerOutput::wall_times) + { + } + + // @sect4{NavierStokes::setup} + template + void NavierStokes::setup() + { + timer.enter_subsection("Setup system"); + dof_handler.distribute_dofs(fe); + DoFRenumbering::Cuthill_McKee(dof_handler); + + // We renumber the components to have all velocity DoFs come before + // the pressure DoFs to be able to split the solution vector in two blocks + // which are separately accessed + std::vector block_component(dim + 1, 0); + block_component[dim] = 1; + DoFRenumbering::component_wise(dof_handler, block_component); + + dofs_per_block.resize(2); + DoFTools::count_dofs_per_block( + dof_handler, dofs_per_block, block_component); + unsigned int dof_u = dofs_per_block[0]; + unsigned int dof_p = dofs_per_block[1]; + + // The Dirichlet boundary condition is applied to boundaries 0 and 1. + FEValuesExtractors::Vector velocities(0); + { + nonzero_constraints.clear(); + + DoFTools::make_hanging_node_constraints(dof_handler, nonzero_constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + BoundaryValues(), + nonzero_constraints, + fe.component_mask(velocities)); + VectorTools::interpolate_boundary_values(dof_handler, + 1, + BoundaryValues(), + nonzero_constraints, + fe.component_mask(velocities)); + } + nonzero_constraints.close(); + + { + zero_constraints.clear(); + + DoFTools::make_hanging_node_constraints(dof_handler, zero_constraints); + VectorTools::interpolate_boundary_values( + dof_handler, + 0, + Functions::ZeroFunction(dim + 1), + zero_constraints, + fe.component_mask(velocities)); + VectorTools::interpolate_boundary_values( + dof_handler, + 1, + Functions::ZeroFunction(dim + 1), + zero_constraints, + fe.component_mask(velocities)); + } + zero_constraints.close(); + + std::cout << " Number of active cells: " << triangulation.n_active_cells() + << std::endl + << " Number of vertices: " << triangulation.n_vertices() + << std::endl + << " Number of degrees of freedom: " << dof_handler.n_dofs() + << " (" << dof_u << '+' << dof_p << ')' << std::endl; + + BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block); + DoFTools::make_sparsity_pattern(dof_handler, dsp, nonzero_constraints); + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit(sparsity_pattern); + mass_matrix.reinit(sparsity_pattern); + + solution.reinit(dofs_per_block); + solution_increment.reinit(dofs_per_block); + system_rhs.reinit(dofs_per_block); + + timer.leave_subsection(); + } + + // @sect4{NavierStokes::setup} + + // A helper function to determine which constrint to use based on the current timestep + template + const ConstraintMatrix &NavierStokes::get_constraints() const + { + return time.get_timestep() == 0 ? nonzero_constraints : zero_constraints; + } + + // @sect4{NavierStokes::assemble} + + // Note that we only need to assemble the LHS for twice: once with the nonzero constraint + // and once for zero constraint. But we must assemble the RHS at every time step. + template + void NavierStokes::assemble(bool assemble_lhs) + { + timer.enter_subsection("Assemble system"); + if (assemble_lhs) + { + system_matrix = 0; + mass_matrix = 0; + } + + system_rhs = 0; + + 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(); + + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + FullMatrix local_mass_matrix(dofs_per_cell, dofs_per_cell); + Vector local_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + std::vector> current_velocity_values(n_q_points); + std::vector> current_velocity_gradients(n_q_points); + std::vector current_velocity_divergences(n_q_points); + std::vector current_pressure_values(n_q_points); + + std::vector div_phi_u(dofs_per_cell); + std::vector> phi_u(dofs_per_cell); + std::vector> grad_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; + local_mass_matrix = 0; + + fe_values[velocities].get_function_values(solution, + current_velocity_values); + + fe_values[velocities].get_function_gradients( + solution, current_velocity_gradients); + + fe_values[velocities].get_function_divergences( + solution, current_velocity_divergences); + + fe_values[pressure].get_function_values(solution, + current_pressure_values); + + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + div_phi_u[k] = fe_values[velocities].divergence(k, q); + grad_phi_u[k] = fe_values[velocities].gradient(k, q); + phi_u[k] = fe_values[velocities].value(k, q); + phi_p[k] = fe_values[pressure].value(k, q); + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + if (assemble_lhs) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + // $LHS = a((u, p), (v, q))*dt + m(u, v) + // = ((grad_v, nu*grad_u) - (div_v, p) - (q, div_u))*dt + + // m(u, v)$ plus Grad-Div term. + local_matrix(i, j) += + ((viscosity * + scalar_product(grad_phi_u[j], grad_phi_u[i]) - + div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j] + + gamma*div_phi_u[j]*div_phi_u[i]) * + time.get_delta_t() + + phi_u[i] * phi_u[j]) * + fe_values.JxW(q); + // Besides the velocity and pressure mass matrices, we also + // assemble $B^T$ and $B$ into the block mass matrix for convenience + // because we need to use them to compute the Schur complement. + // As a result $M = [M_u, B^T; B, M_p]$. + local_mass_matrix(i, j) += + (phi_u[i] * phi_u[j] + phi_p[i] * phi_p[j] - + div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) * + fe_values.JxW(q); + } + } + // $RHS = - dt*[ a((u_prev, p_prev), (v, q)) + c(u_prev; u_prev, v)]$ + // plus Grad-Div term. + local_rhs(i) -= + (viscosity * scalar_product(current_velocity_gradients[q], + grad_phi_u[i]) - + current_velocity_divergences[q] * phi_p[i] - + current_pressure_values[q] * div_phi_u[i] + + current_velocity_gradients[q] * current_velocity_values[q] * + phi_u[i] + + gamma * current_velocity_divergences[q] * div_phi_u[i]) * + fe_values.JxW(q) * time.get_delta_t(); + } + } + + cell->get_dof_indices(local_dof_indices); + + const ConstraintMatrix &constraints_used = get_constraints(); + + if (assemble_lhs) + { + constraints_used.distribute_local_to_global(local_matrix, + local_rhs, + local_dof_indices, + system_matrix, + system_rhs); + constraints_used.distribute_local_to_global(local_mass_matrix, + local_dof_indices, + mass_matrix); + } + else + { + constraints_used.distribute_local_to_global( + local_rhs, local_dof_indices, system_rhs); + } + } + timer.leave_subsection(); + } + + // @sect4{NavierStokes::solve_linear_system} + + // Only updates the preconditioners when we assemble the LHS of the system. + template + std::pair NavierStokes::solve_linear_system(bool update_precondition) + { + const ConstraintMatrix &constraints_used = get_constraints(); + + if (update_precondition) + { + timer.enter_subsection("Precondition linear system"); + + preconditioner.reset(); + A_inverse.reset(); + preconditioner_A.reset(); + S_inverse.reset(); + Mp_inverse.reset(); + preconditioner_Mp.reset(); + Sm_inverse.reset(); + preconditioner_Sm.reset(); + approximate_Sm.reset(); + + approximate_Sm.reset(new ApproximateMassSchur(mass_matrix)); + preconditioner_Sm.reset(new PreconditionIdentity()); + Sm_inverse.reset(new InverseMatrix + (*approximate_Sm, *preconditioner_Sm)); + preconditioner_Mp.reset(new SparseILU()); + preconditioner_Mp->initialize(mass_matrix.block(1,1)); + Mp_inverse.reset(new InverseMatrix, SparseILU> + (mass_matrix.block(1,1), *preconditioner_Mp)); + S_inverse.reset(new SchurComplementInverse>(gamma, viscosity, time.get_delta_t(), *Sm_inverse, *Mp_inverse)); + preconditioner_A.reset(new typename InnerPreconditioner::type()); + preconditioner_A->initialize(system_matrix.block(0,0), + typename InnerPreconditioner::type::AdditionalData()); + + A_inverse.reset(new InverseMatrix, + typename InnerPreconditioner::type>(system_matrix.block(0,0), *preconditioner_A)); + preconditioner.reset(new BlockSchurPreconditioner< + typename InnerPreconditioner::type, PreconditionIdentity, + SparseILU>(system_matrix, *A_inverse, *S_inverse)); + + timer.leave_subsection(); + } + + // Solve with GMRES solver. + timer.enter_subsection("Solve linear system"); + SolverControl solver_control(system_matrix.m(), + 1e-8 * 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_increment, system_rhs, *preconditioner); + + constraints_used.distribute(solution_increment); + timer.leave_subsection(); + + return {solver_control.last_step(), solver_control.last_value()}; + } + + // @sect4{NavierStokes::run} + + template + void NavierStokes::run() + { + create_triangulation(triangulation); + triangulation.refine_global(2); + setup(); + + std::ofstream out("grid.eps"); + GridOut grid_out; + grid_out.write_eps(triangulation, out); + + std::ofstream out2("force.txt"); + out2 << std::setw(13) << std::left << "Time/s" + << std::setw(13) << std::left << " Drag" << std::setw(13) + << std::left << " Lift" << std::endl; + + // In IMEX scheme we do not need to implement the Newton's method, what we need + // to do at every time step is simple: + // 1. Solve for the solution increment; 2. Update the solution. + output_results(time.get_timestep()); + while (time.current() <= time.end()) + { + std::cout << "*****************************************" << std::endl; + std::cout << "Time = " << time.current() << std::endl; + + assemble(time.get_timestep() < 2); + + auto state = solve_linear_system(time.get_timestep() < 2); + solution.add(1.0, solution_increment); + + // solution is distributed using nonzero_constraints all the time + nonzero_constraints.distribute(solution); + solution_increment = 0; + + std::cout << " FGMRES steps = " << state.first + << " residual = " << std::setw(6) << state.second << std::endl; + + time.increment(); + + if (time.get_timestep() % 1 == 0) + { + output_results(time.get_timestep()); + process_solution(out2); + } + } + + out2.close(); + } + + // @sect4{NavierStokes::output_result} + + template + void NavierStokes::output_results(const unsigned int output_index) const + { + timer.enter_subsection("Output"); + std::cout << " Writing results..." << std::endl; + 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 << "Re100-" + << Utilities::int_to_string(output_index, 6) << ".vtu"; + + std::ofstream output(filename.str().c_str()); + data_out.write_vtu(output); + timer.leave_subsection(); + } + + // @sect4{NavierStokes::process_solution} + + // This function is used to calculate the drag and lift coefficients on the cylinder. + // We first calculate the traction of the fluid, which is nothing but the product of the + // stress tensor and the normal of the cylindrical surface, and then integrate it along + // the cylindrical surface and negate it. + template + void NavierStokes::process_solution(std::ofstream& out) const + { + timer.enter_subsection("Process solution"); + + Tensor<1, dim> force; + + FEFaceValues fe_face_values(fe, + face_quadrature_formula, + update_values | update_quadrature_points | + update_JxW_values | update_normal_vectors | + update_gradients); + + const unsigned int n_q_points = face_quadrature_formula.size(); + + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + + std::vector p(n_q_points); + std::vector> grad_sym_v(n_q_points); + + for (auto cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell) + { + for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f) + { + if (cell->face(f)->at_boundary() && cell->face(f)->boundary_id() == 1) + { + fe_face_values.reinit(cell, f); + fe_face_values[pressure].get_function_values(solution, p); + fe_face_values[velocities].get_function_symmetric_gradients(solution, grad_sym_v); + for (unsigned int q = 0; q < n_q_points; ++q) + { + const Tensor<1, dim> &N = fe_face_values.normal_vector(q); + SymmetricTensor<2, dim> stress = -p[q]*Physics::Elasticity::StandardTensors::I + + viscosity*grad_sym_v[q]; + force -= stress*N*fe_face_values.JxW(q); + } + } + } + } + + double drag_coef = 2*force[0]/(0.1); + double lift_coef = 2*force[dim-1]/(0.1); + + out.precision(6); + out.width(12); + + out << std::scientific << std::left << + time.current() << " " << drag_coef << " " << lift_coef << std::endl; + + timer.leave_subsection(); + } +} + +// @sect3{main function} + +int main() +{ + try + { + using namespace dealii; + using namespace fluid; + + NavierStokes<2> flow(/* degree = */ 1); + flow.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; +}