{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# Demo - Stokes equations\n", "\n", " \n", "**Mikael Mortensen** (email: `mikaem@math.uio.no`), Department of Mathematics, University of Oslo.\n", "\n", "Date: **January 23, 2019**\n", "\n", "**Summary.** The Stokes equations describe the flow of highly viscous fluids.\n", "This is a demonstration of how the Python module [shenfun](https://github.com/spectralDNS/shenfun) can be used to solve Stokes\n", "equations using a mixed (coupled) basis in a 3D tensor product domain.\n", "We assume homogeneous Dirichlet boundary conditions in one direction\n", "and periodicity in the remaining two. The solver described runs with MPI\n", "without any further considerations required from the user.\n", "The solver assembles a block matrix with sparsity pattern as shown below\n", "for the Legendre basis.\n", "\n", "\n", "\n", "\n", "\n", "

Coupled block matrix for Stokes equations.

\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Stokes' equations\n", "\n", "\n", "Stokes' equations are given in strong form as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\nabla^2 \\boldsymbol{u} - \\nabla p &= \\boldsymbol{f} \\quad \\text{in } \\Omega, \\\\ \n", "\\nabla \\cdot \\boldsymbol{u} &= h \\quad \\text{in } \\Omega \\\\ \n", "\\int_{\\Omega} p dx &= 0\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $\\boldsymbol{u}$ and $p$ are, respectively, the\n", "fluid velocity vector and pressure, and the domain\n", "$\\Omega = [0, 2\\pi]^2 \\times [-1, 1]$. The flow is assumed periodic\n", "in $x$ and $y$-directions, whereas there is a no-slip homogeneous Dirichlet\n", "boundary condition on $\\boldsymbol{u}$ on the boundaries of the $z$-direction, i.e.,\n", "$\\boldsymbol{u}(x, y, \\pm 1) = (0, 0, 0)$. (Note that we can configure shenfun with\n", "non-periodicity in any of the three directions. However, since we are to\n", "solve linear algebraic systems in the non-periodic direction, there is a speed\n", "benefit from having the nonperiodic direction last. This has to do with Numpy\n", "using a C-style row-major storage of arrays by default.)\n", "The right hand side vector $\\boldsymbol{f}(\\boldsymbol{x})$ is an external applied body force.\n", "The right hand side $h$ is usually zero in the regular Stokes equations. Here\n", "we include it because it will be nonzero in the verification, which is using the\n", "method of manufactured solutions. Note that the final $\\int_{\\Omega} p dx = 0$\n", "is there because there is no Dirichlet boundary condition on the pressure\n", "and the system of equations would otherwise be ill conditioned.\n", "\n", "To solve Stokes' equations with the Galerkin method we need basis\n", "functions for both velocity and pressure. A\n", "Dirichlet basis will be used for velocity, whereas there is no boundary restriction\n", "on the pressure basis. For both three-dimensional bases we will use one basis\n", "function for the $x$-direction,\n", "$\\mathcal{X}(x)$, one for the $y$-direction, $\\mathcal{Y}(y)$, and one for the\n", "$z$-direction, $\\mathcal{Z}(z)$. And\n", "then we create three-dimensional basis functions like" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "v(x, y, z) = \\mathcal{X}(x) \\mathcal{Y}(y) \\mathcal{Z} (z).\n", "\\label{_auto1} \\tag{1}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The basis functions $\\mathcal{X}(x)$ and $\\mathcal{Y}(y)$ are chosen as Fourier\n", "exponentials, since these functions are periodic:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\mathcal{X}_l(x) = e^{\\imath l x}, \\forall \\, l \\in \\boldsymbol{l}^{N_0}, \n", "\\label{_auto2} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "\\mathcal{Y}_m(y) = e^{\\imath m y}, \\forall \\, m \\in \\boldsymbol{m}^{N_1},\n", "\\label{_auto3} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $\\boldsymbol{l}^{N_0} = (-N_0/2, -N_0/2+1, \\ldots, N_0/2-1)$ and\n", "$\\boldsymbol{m}^{N_1} = (-N_1/2, -N_1/2+1, \\ldots, N_1/2-1)$.\n", "The size of the discretized problem in real physical space is\n", "$\\boldsymbol{N} = (N_0, N_1, N_2)$, i.e., there are $N_0 \\cdot N_1 \\cdot N_2$ quadrature points\n", "in total.\n", "\n", "The basis functions for $\\mathcal{Z}(z)$ remain to be decided.\n", "For the velocity we need homogeneous Dirichlet boundary conditions, and for this\n", "we use composite Legendre or Chebyshev polynomials" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\mathcal{Z}^0_n(z) = \\phi_n(z) - \\phi_{n+2}(z), \\forall \\, n \\in \\boldsymbol{n}^{N_2-2},\n", "\\label{_auto4} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $\\phi_n$ is the n'th Legendre or Chebyshev polynomial of the first kind.\n", "$\\boldsymbol{n}^{N_2-2} = (0, 1, \\ldots, N_2-3)$, and the zero on $\\mathcal{Z}^0$\n", "is there to indicate the zero value on the boundary.\n", "\n", "The pressure basis that comes with no restrictions for the boundary is a\n", "little trickier. The reason for this has to do with\n", "inf-sup stability. The obvious choice of basis is the regular Legendre or\n", "Chebyshev basis, which is denoted as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\mathcal{Z}_n(z) = \\phi_n(z), \\forall \\, n \\in \\boldsymbol{n}^{N_2}. \\label{eq:Zn} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The problem is that for the natural choice of $n \\in (0, 1, \\ldots, N_2-1)$\n", "there is a nullspace and one degree of freedom remains unresolved. It turns out\n", "that the proper choice for the pressure basis is simply ([5](#eq:Zn)) for\n", "$n \\in \\boldsymbol{n}^{N_2-2}$. (Also remember that we have to fix $\\int_{\\Omega} p dx = 0$.)\n", "\n", "With given basis functions we obtain the spaces" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "V^{N_0} = \\text{span}\\{ \\mathcal{X}_l \\}_{l\\in\\boldsymbol{l}^{N_0}}, \n", "\\label{_auto5} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "V^{N_1} = \\text{span}\\{ \\mathcal{Y}_m \\}_{m\\in\\boldsymbol{m}^{N_1}}, \n", "\\label{_auto6} \\tag{7}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "V^{N_2} = \\text{span}\\{ \\mathcal{Z}_n \\}_{n\\in\\boldsymbol{n}^{N_2-2}}, \n", "\\label{_auto7} \\tag{8}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "V_0^{N_2} = \\text{span}\\{ \\mathcal{Z}^0_n \\}_{n\\in\\boldsymbol{n}^{N_2-2}},\n", "\\label{_auto8} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and from these we create two different tensor product spaces" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "W_0^{\\boldsymbol{N}}(\\boldsymbol{x}) = V^{N_0}(x) \\otimes V^{N_1}(y) \\otimes V_0^{N_2}(z), \n", "\\label{_auto9} \\tag{10}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "W^{\\boldsymbol{N}}(\\boldsymbol{x}) = V^{N_0}(x) \\otimes V^{N_1}(y) \\otimes V^{N_2}(z).\n", "\\label{_auto10} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The velocity vector is using a mixed basis, such that we will look for\n", "solutions $\\boldsymbol{u} \\in [W_0^{\\boldsymbol{N}}]^3 \\, (=W_0^{\\boldsymbol{N}} \\times W_0^{\\boldsymbol{N}} \\times W_0^{\\boldsymbol{N}})$,\n", "whereas we look for the pressure\n", "$p \\in W^{\\boldsymbol{N}}$. We now formulate a variational problem using the Galerkin method: Find\n", "$\\boldsymbol{u} \\in [W_0^{\\boldsymbol{N}}]^3$ and $p \\in W^{\\boldsymbol{N}}$ such that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\int_{\\Omega} (\\nabla^2 \\boldsymbol{u} - \\nabla p ) \\cdot \\overline{\\boldsymbol{v}} \\, dx_w = \\int_{\\Omega} \\boldsymbol{f} \\cdot \\overline{\\boldsymbol{v}}\\, dx_w \\quad\\forall \\boldsymbol{v} \\, \\in \\, [W_0^{\\boldsymbol{N}}]^3, \\label{eq:varform} \\tag{12} \n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "\\int_{\\Omega} \\nabla \\cdot \\boldsymbol{u} \\, \\overline{q} \\, dx_w = \\int_{\\Omega} h \\overline{q} \\, dx_w \\quad\\forall q \\, \\in \\, W^{\\boldsymbol{N}}.\n", "\\label{_auto11} \\tag{13}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here $dx_w=w_xdxw_ydyw_zdz$ represents a weighted measure, with weights $w_x(x), w_y(y), w_z(z)$.\n", "Note that it is only Chebyshev polynomials that\n", "make use of a non-constant weight $w_x=1/\\sqrt{1-x^2}$. The Fourier weights are $w_y=w_z=1/(2\\pi)$\n", "and the Legendre weight is $w_x=1$.\n", "The overline in $\\boldsymbol{\\overline{v}}$ and $\\overline{q}$ represents a complex conjugate, which is needed here because\n", "the Fourier exponentials are complex functions.\n", "\n", "### Mixed variational form\n", "\n", "\n", "\n", "Since we are to solve for $\\boldsymbol{u}$ and $p$ at the same time, we formulate a\n", "mixed (coupled) problem: find $(\\boldsymbol{u}, p) \\in [W_0^{\\boldsymbol{N}}]^3 \\times W^{\\boldsymbol{N}}$\n", "such that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "a((\\boldsymbol{u}, p), (\\boldsymbol{v}, q)) = L((\\boldsymbol{v}, q)) \\quad \\forall (\\boldsymbol{v}, q) \\in [W_0^{\\boldsymbol{N}}]^3 \\times W^{\\boldsymbol{N}},\n", "\\label{_auto12} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where bilinear ($a$) and linear ($L$) forms are given as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " a((\\boldsymbol{u}, p), (\\boldsymbol{v}, q)) = \\int_{\\Omega} (\\nabla^2 \\boldsymbol{u} - \\nabla p) \\cdot \\overline{\\boldsymbol{v}} \\, dx_w + \\int_{\\Omega} \\nabla \\cdot \\boldsymbol{u} \\, \\overline{q} \\, dx_w, \n", "\\label{_auto13} \\tag{15}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", " L((\\boldsymbol{v}, q)) = \\int_{\\Omega} \\boldsymbol{f} \\cdot \\overline{\\boldsymbol{v}}\\, dx_w + \\int_{\\Omega} h \\overline{q} \\, dx_w.\n", "\\label{_auto14} \\tag{16}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the bilinear form will assemble to block matrices, whereas the right hand side\n", "linear form will assemble to block vectors.\n", "\n", "## Implementation\n", "\n", "### Preamble\n", "\n", "We will solve the Stokes equations using the [shenfun](https://github.com/spectralDNS/shenfun) Python module. The first thing needed\n", "is then to import some of this module's functionality\n", "plus some other helper modules, like [Numpy](https://numpy.org) and [Sympy](https://sympy.org):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "import sys\n", "import numpy as np\n", "from sympy import symbols, sin, cos\n", "from shenfun import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use `Sympy` for the manufactured solution and `Numpy` for testing.\n", "\n", "### Manufactured solution\n", "\n", "\n", "\n", "The exact solutions $\\boldsymbol{u}_e(\\boldsymbol{x})$ and $p(\\boldsymbol{x})$ are chosen to satisfy boundary\n", "conditions, and the right hand sides $\\boldsymbol{f}(\\boldsymbol{x})$ and $h(\\boldsymbol{x})$ are then\n", "computed exactly using `Sympy`. These exact right hand sides will then be used to\n", "compute a numerical solution that can be verified against the manufactured\n", "solution. The chosen solution with computed right hand sides are:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x, y, z = symbols('x,y,z')\n", "uex = sin(2*y)*(1-z**2)\n", "uey = sin(2*x)*(1-z**2)\n", "uez = sin(2*z)*(1-z**2)\n", "pe = -0.1*sin(2*x)*cos(4*y)\n", "fx = uex.diff(x, 2) + uex.diff(y, 2) + uex.diff(z, 2) - pe.diff(x, 1)\n", "fy = uey.diff(x, 2) + uey.diff(y, 2) + uey.diff(z, 2) - pe.diff(y, 1)\n", "fz = uez.diff(x, 2) + uez.diff(y, 2) + uez.diff(z, 2) - pe.diff(z, 1)\n", "h = uex.diff(x, 1) + uey.diff(y, 1) + uez.diff(z, 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Tensor product spaces\n", "\n", "One-dimensional spaces are created using the [FunctionSpace()](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.FunctionSpace) function. A choice of\n", "polynomials between Legendre or Chebyshev can be made, and the size\n", "of the domain is given" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = (20, 20, 20)\n", "family = 'Legendre'\n", "K0 = FunctionSpace(N[0], 'Fourier', dtype='D', domain=(0, 2*np.pi))\n", "K1 = FunctionSpace(N[1], 'Fourier', dtype='d', domain=(0, 2*np.pi))\n", "SD = FunctionSpace(N[2], family, bc=(0, 0))\n", "ST = FunctionSpace(N[2], family)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next the one-dimensional spaces are used to create two tensor product spaces Q = $W^{\\boldsymbol{N}}$\n", "and TD = $W_0^{\\boldsymbol{N}}$, one vector V = $[W_0^{\\boldsymbol{N}}]^3$ and one mixed\n", "space VQ = V $\\times$ Q." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "TD = TensorProductSpace(comm, (K0, K1, SD), axes=(2, 0, 1))\n", "Q = TensorProductSpace(comm, (K0, K1, ST), axes=(2, 0, 1))\n", "V = VectorSpace(TD)\n", "VQ = CompositeSpace([V, Q])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we choose to transform axes in the order $1, 0, 2$. This is to ensure\n", "that the fully transformed arrays are aligned in the non-periodic direction 2.\n", "And we need the arrays aligned in this direction, because this is the only\n", "direction where there are tensor product matrices that are non-diagonal. All\n", "Fourier matrices are, naturally, diagonal.\n", "\n", "Test- and trialfunctions are created much like in a regular, non-mixed,\n", "formulation. However, one has to create one test- and trialfunction for\n", "the mixed space, and then split them up afterwards" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "up = TrialFunction(VQ)\n", "vq = TestFunction(VQ)\n", "u, p = up\n", "v, q = vq" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the basisfunctions in place we may assemble the different blocks of the\n", "final coefficient matrix. Since Legendre is using a constant weight function,\n", "the equations may also be integrated by parts to obtain a symmetric system:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if family.lower() == 'chebyshev':\n", " A = inner(v, div(grad(u)))\n", " G = inner(v, -grad(p))\n", "else:\n", " A = inner(grad(v), -grad(u))\n", " G = inner(div(v), p)\n", "D = inner(q, div(u))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The assembled subsystems `A, G` and `D` are lists containg the different blocks of\n", "the complete, coupled matrix. `A` actually contains 6\n", "tensor product matrices of type [TPMatrix](https://shenfun.readthedocs.io/en/latest/shenfun.html#shenfun.matrixbase.TPMatrix). The first two\n", "matrices are for vector component zero of the test function `v[0]` and\n", "trial function `u[0]`, the\n", "matrices 2 and 3 are for components 1 and the last two are for components\n", "2. The first two matrices are as such for" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```Python\n", " A[0:2] = inner(v[0], div(grad(u[0])))\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Breaking it down this inner product is mathematically" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\label{eq:partialeq1} \\tag{17}\n", "\\int_{\\Omega} \\boldsymbol{v}[0] \\left(\\frac{\\partial^2 \\boldsymbol{u}[0]}{\\partial x^2} + \\frac{\\partial^2 \\boldsymbol{u}[0]}{\\partial y^2} + \\frac{\\partial^2 \\boldsymbol{u}[0]}{\\partial z^2}\\right) w_x dx w_y dy w_z dz.\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we now use test function $\\boldsymbol{v}[0]$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\boldsymbol{v}[0]_{lmn} = \\mathcal{X}_l \\mathcal{Y}_m \\mathcal{Z}_n,\n", "\\label{_auto15} \\tag{18}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and trialfunction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\boldsymbol{u}[0]_{pqr} = \\sum_{p} \\sum_{q} \\sum_{r} \\hat{\\boldsymbol{u}}[0]_{pqr} \\mathcal{X}_p \\mathcal{Y}_q \\mathcal{Z}_r,\n", "\\label{_auto16} \\tag{19}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $\\hat{\\boldsymbol{u}}$ are the unknown degrees of freedom, and then insert these functions\n", "into ([17](#eq:partialeq1)), then we obtain after\n", "performing some exact evaluations over the periodic directions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " \\Big( \\underbrace{-\\left(l^2 \\delta_{lp} + m^2 \\delta_{mq} \\right) \\int_{-1}^{1} \\mathcal{Z}_r(z) \\mathcal{Z}_n(z) w_z dz}_{A[0]} + \\underbrace{\\delta_{lp} \\delta_{mq} \\int_{-1}^{1} \\frac{\\partial^2 \\mathcal{Z}_r(z)}{\\partial z^2} \\mathcal{Z}_n(z) w_z dz}_{A[1]} \\Big) \\hat{\\boldsymbol{u}}[0]_{pqr},\n", "\\label{_auto17} \\tag{20}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly for components 1 and 2 of the test and trial vectors, leading to 6 tensor\n", "product matrices in total for `A`. Similarly, we get three components of `G`\n", "and three of `D`.\n", "\n", "Eliminating the Fourier diagonal matrices, we are left with block matrices like" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "H(l, m) =\n", " \\begin{bmatrix}\n", " A[0]+A[1] & 0 & 0 & G[0] \\\\ \n", " 0 & A[2]+A[3] & 0 & G[1] \\\\ \n", " 0 & 0 & A[4]+A[5] & G[2] \\\\ \n", " D[0] & D[1] & D[2] & 0\n", " \\end{bmatrix}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that there will be one large block matrix $H(l, m)$ for each Fourier\n", "wavenumber combination $(l, m)$. To solve the problem in the end we will need to\n", "loop over these wavenumbers and solve the assembled linear systems one by one.\n", "An example of the block matrix, for $l=m=5$ and $\\boldsymbol{N}=(20, 20, 20)$ is given\n", "in Fig. [fig:BlockMat](#fig:BlockMat).\n", "\n", "\n", "In the end we create a block matrix through" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "M = BlockMatrix(A+G+D)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The right hand side can easily be assembled since we have already\n", "defined the functions $\\boldsymbol{f}$ and $h$, see Sec. [Manufactured solution](#sec:mansol)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Get mesh (quadrature points)\n", "X = TD.local_mesh(True)\n", "\n", "# Get f and h on quad points\n", "fh = Array(VQ, buffer=(fx, fy, fz, h))\n", "f_, h_ = fh\n", "\n", "# Compute inner products\n", "fh_hat = Function(VQ)\n", "f_hat, h_hat = fh_hat\n", "f_hat = inner(v, f_, output_array=f_hat)\n", "h_hat = inner(q, h_, output_array=h_hat)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the end all that is left is to solve and compare with\n", "the exact solution." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Solve problem\n", "up_hat = M.solve(fh_hat, constraints=((3, 0, 0), (3, N[2]-1, 0)))\n", "up = up_hat.backward()\n", "u_, p_ = up\n", "\n", "# Exact solution\n", "ux, uy, uz = Array(V, buffer=(uex, uey, uez))\n", "pe = Array(Q, buffer=pe)\n", "\n", "error = [comm.reduce(np.linalg.norm(ux-u_[0])),\n", " comm.reduce(np.linalg.norm(uy-u_[1])),\n", " comm.reduce(np.linalg.norm(uz-u_[2])),\n", " comm.reduce(np.linalg.norm(pe-p_))]\n", "print(error)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that solve has a keyword argument\n", "`constraints=((3, 0, 0), (3, N[2]-1), 0)` that takes care of the restriction\n", "$\\int_{\\Omega} p dx = 0$ by indenting the rows in M corresponding to the\n", "first and last degree of freedom for the pressure. The value $(3, 0, 0)$\n", "indicates that pressure is\n", "in block 3 of the block vector solution (the velocity vector holds\n", "positions 0, 1 and 2), whereas the two zeros ensures that the first dof\n", "(dof 0) should obtain value 0. The constraint on the highest\n", "wavenumber `(3, N[2]-1, 0)` is required to get a non-singular\n", "matrix.\n", "\n", "## Complete solver\n", "\n", "\n", "A complete solver can be found in demo [Stokes3D.py](https://github.com/spectralDNS/shenfun/blob/master/demo/Stokes3D.py)." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.6" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": false, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }