{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Tutorial 8: Incompressible Navier-Stokes"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"In this tutorial, we will learn\n",
" - How to solve nonlinear multi-field PDEs in Gridap\n",
" - How to build FE spaces whose functions have zero mean value\n",
"\n",
"## Problem statement\n",
"\n",
"The goal of this last tutorial is to solve a nonlinear multi-field PDE. As a model problem, we consider a well known benchmark in computational fluid dynamics, the lid-driven cavity for the incompressible Navier-Stokes equations. Formally, the PDE we want to solve is: find the velocity vector $u$ and the pressure $p$ such that\n",
"\n",
"$$\n",
"\\left\\lbrace\n",
"\\begin{aligned}\n",
"-\\Delta u + \\mathit{Re}\\ (u\\cdot \\nabla)\\ u + \\nabla p = 0 &\\text{ in }\\Omega,\\\\\n",
"\\nabla\\cdot u = 0 &\\text{ in } \\Omega,\\\\\n",
"u = g &\\text{ on } \\partial\\Omega,\n",
"\\end{aligned}\n",
"\\right.\n",
"$$\n",
"\n",
"where the computational domain is the unit square $\\Omega \\doteq (0,1)^d$, $d=2$, $\\mathit{Re}$ is the Reynolds number (here, we take $\\mathit{Re}=10$), and $(w \\cdot \\nabla)\\ u = (\\nabla u)^t w$ is the well known convection operator. In this example, the driving force is the Dirichlet boundary velocity $g$, which is a non-zero horizontal velocity with a value of $g = (1,0)^t$ on the top side of the cavity, namely the boundary $(0,1)\\times\\{1\\}$, and $g=0$ elsewhere on $\\partial\\Omega$. Since we impose Dirichlet boundary conditions on the entire boundary $\\partial\\Omega$, the mean value of the pressure is constrained to zero in order have a well posed problem,\n",
"\n",
"$$\n",
"\\int_\\Omega q \\ {\\rm d}\\Omega = 0.\n",
"$$\n",
"\n",
"## Numerical Scheme\n",
"\n",
"In order to approximate this problem we chose a formulation based on inf-sub stable $Q_k/P_{k-1}$ elements with continuous velocities and discontinuous pressures (see, e.g., [1] for specific details). The interpolation spaces are defined as follows. The velocity interpolation space is\n",
"\n",
"$$\n",
"V \\doteq \\{ v \\in [C^0(\\Omega)]^d:\\ v|_T\\in [Q_k(T)]^d \\text{ for all } T\\in\\mathcal{T} \\},\n",
"$$\n",
"where $T$ denotes an arbitrary cell of the FE mesh $\\mathcal{T}$, and $Q_k(T)$ is the local polynomial space in cell $T$ defined as the multi-variate polynomials in $T$ of order less or equal to $k$ in each spatial coordinate. Note that, this is the usual continuous vector-valued Lagrangian FE space of order $k$ defined on a mesh of quadrilaterals or hexahedra. On the other hand, the space for the pressure is\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"Q_0 &\\doteq \\{ q \\in Q: \\ \\int_\\Omega q \\ {\\rm d}\\Omega = 0\\}, \\text{ with}\\\\\n",
"Q &\\doteq \\{ q \\in L^2(\\Omega):\\ q|_T\\in P_{k-1}(T) \\text{ for all } T\\in\\mathcal{T}\\},\n",
"\\end{aligned}\n",
"$$\n",
"where $P_{k-1}(T)$ is the polynomial space of multi-variate polynomials in $T$ of degree less or equal to $k-1$. Note that functions in $Q_0$ are strongly constrained to have zero mean value. This is achieved in the code by removing one degree of freedom from the (unconstrained) interpolation space $Q$ and adding a constant to the computed pressure so that the resulting function has zero mean value.\n",
"\n",
"The weak form associated to these interpolation spaces reads: find $(u,p)\\in U_g \\times Q_0$ such that $[r(u,p)](v,q)=0$ for all $(v,q)\\in V_0 \\times Q_0$\n",
"where $U_g$ and $V_0$ are the set of functions in $V$ fulfilling the Dirichlet boundary condition $g$ and $0$ on $\\partial\\Omega$ respectively. The weak residual $r$ evaluated at a given pair $(u,p)$ is the linear form defined as\n",
"\n",
"$$\n",
"[r(u,p)](v,q) \\doteq a((u,p),(v,q))+ [c(u)](v),\n",
"$$\n",
"with\n",
"$$\n",
"\\begin{aligned}\n",
"a((u,p),(v,q)) &\\doteq \\int_{\\Omega} \\nabla v \\cdot \\nabla u \\ {\\rm d}\\Omega - \\int_{\\Omega} (\\nabla\\cdot v) \\ p \\ {\\rm d}\\Omega + \\int_{\\Omega} q \\ (\\nabla \\cdot u) \\ {\\rm d}\\Omega,\\\\\n",
"[c(u)](v) &\\doteq \\int_{\\Omega} v \t\\cdot \\left( (u\\cdot\\nabla)\\ u \\right)\\ {\\rm d}\\Omega.\\\\\n",
"\\end{aligned}\n",
"$$\n",
"Note that the bilinear form $a$ is associated with the linear part of the PDE, whereas $c$ is the contribution to the residual resulting from the convective term.\n",
"\n",
"In order to solve this nonlinear weak equation with a Newton-Raphson method, one needs to compute the Jacobian associated with the residual $r$. In this case, the Jacobian $j$ evaluated at a pair $(u,p)$ is the bilinear form defined as\n",
"\n",
"$$\n",
"[j(u,p)]((\\delta u, \\delta p),(v,q)) \\doteq a((\\delta u,\\delta p),(v,q)) + [{\\rm d}c(u)](\\delta u,v),\n",
"$$\n",
"where ${\\rm d}c$ results from the linearization of the convective term, namely\n",
"$$\n",
"[{\\rm d}c(u)](\\delta u,v) \\doteq \\int_{\\Omega} v \\cdot \\left( (u\\cdot\\nabla)\\ \\delta u \\right) \\ {\\rm d}\\Omega + \\int_{\\Omega} v \\cdot \\left( (\\delta u\\cdot\\nabla)\\ u \\right) \\ {\\rm d}\\Omega.\n",
"$$\n",
"The implementation of this numerical scheme is done in Gridap by combining the concepts previously seen for single-field nonlinear PDEs and linear multi-field problems.\n",
"\n",
"## Discrete model\n",
"\n",
"We start with the discretization of the computational domain. We consider a $100\\times100$ Cartesian mesh of the unit square."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using Gridap\n",
"n = 100\n",
"domain = (0,1,0,1)\n",
"partition = (n,n)\n",
"model = CartesianDiscreteModel(domain,partition)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"For convenience, we create two new boundary tags, namely `\"diri1\"` and `\"diri0\"`, one for the top side of the square (where the velocity is non-zero), and another for the rest of the boundary (where the velocity is zero)."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"labels = get_face_labeling(model)\n",
"add_tag_from_tags!(labels,\"diri1\",[6,])\n",
"add_tag_from_tags!(labels,\"diri0\",[1,2,3,4,5,7,8])"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"## FE spaces\n",
"\n",
"For the velocities, we need to create a conventional vector-valued continuous Lagrangian FE space. In this example, we select a second order interpolation."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"D = 2\n",
"order = 2\n",
"reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)\n",
"V = TestFESpace(model,reffeᵤ,conformity=:H1,labels=labels,dirichlet_tags=[\"diri0\",\"diri1\"])"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"The interpolation space for the pressure is built as follows"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"reffeₚ = ReferenceFE(lagrangian,Float64,order-1;space=:P)\n",
"Q = TestFESpace(model,reffeₚ,conformity=:L2,constraint=:zeromean)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"With the options `:Lagrangian`, `space=:P`, `valuetype=Float64`, and `order=order-1`, we select the local polynomial space $P_{k-1}(T)$ on the cells $T\\in\\mathcal{T}$. With the symbol `space=:P` we specifically chose a local Lagrangian interpolation of type \"P\". Without using `space=:P`, would lead to a local Lagrangian of type \"Q\" since this is the default for quadrilateral or hexahedral elements. On the other hand, `constraint=:zeromean` leads to a FE space, whose functions are constrained to have mean value equal to zero, which is just what we need for the pressure space. With these objects, we build the trial multi-field FE spaces"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"uD0 = VectorValue(0,0)\n",
"uD1 = VectorValue(1,0)\n",
"U = TrialFESpace(V,[uD0,uD1])\n",
"P = TrialFESpace(Q)\n",
"\n",
"Y = MultiFieldFESpace([V, Q])\n",
"X = MultiFieldFESpace([U, P])"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"## Triangulation and integration quadrature\n",
"\n",
"From the discrete model we can define the triangulation and integration measure"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"degree = order\n",
"Ωₕ = Triangulation(model)\n",
"dΩ = Measure(Ωₕ,degree)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"## Nonlinear weak form\n",
"\n",
"The different terms of the nonlinear weak form for this example are defined following an approach similar to the one discussed for the $p$-Laplacian equation, but this time using the notation for multi-field problems."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"const Re = 10.0\n",
"conv(u,∇u) = Re*(∇u')⋅u\n",
"dconv(du,∇du,u,∇u) = conv(u,∇du)+conv(du,∇u)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"The bilinear form reads"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"a((u,p),(v,q)) = ∫( ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u) )dΩ"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"The nonlinear term and its Jacobian are given by"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"c(u,v) = ∫( v⊙(conv∘(u,∇(u))) )dΩ\n",
"dc(u,du,v) = ∫( v⊙(dconv∘(du,∇(du),u,∇(u))) )dΩ"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"Finally, the Navier-Stokes weak form residual and Jacobian can be defined as"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"res((u,p),(v,q)) = a((u,p),(v,q)) + c(u,v)\n",
"jac((u,p),(du,dp),(v,q)) = a((du,dp),(v,q)) + dc(u,du,v)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"With the functions `res`, and `jac` representing the weak residual and the Jacobian, we build the nonlinear FE problem:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"op = FEOperator(res,jac,X,Y)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"## Nonlinear solver phase\n",
"\n",
"To finally solve the problem, we consider the same nonlinear solver as previously considered for the $p$-Laplacian equation."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using LineSearches: BackTracking\n",
"nls = NLSolver(\n",
" show_trace=true, method=:newton, linesearch=BackTracking())\n",
"solver = FESolver(nls)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"In this example, we solve the problem without providing an initial guess (a default one equal to zero will be generated internally)"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"uh, ph = solve(solver,op)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"Finally, we write the results for visualization (see next figure)."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"writevtk(Ωₕ,\"ins-results\",cellfields=[\"uh\"=>uh,\"ph\"=>ph])"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"\n",
"\n",
" ## References\n",
"\n",
" [1] H. C. Elman, D. J. Silvester, and A. J. Wathen. *Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics*. Oxford University Press, 2005."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"---\n",
"\n",
"*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
],
"metadata": {}
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.10.8"
},
"kernelspec": {
"name": "julia-1.10",
"display_name": "Julia 1.10.8",
"language": "julia"
}
},
"nbformat": 4
}