![](\"https://raw.githack.com/spectralDNS/spectralutilities/master/figures/Helmholtzdisc.png\")
Figure 1: Helmholtz on the unit disc.
\n", "" ] }, { "cell_type": "markdown", "id": "a8436390", "metadata": { "editable": true }, "source": [ "## Helmholtz equation\n", "\n", "\n", "The Helmholtz equation is given as" ] }, { "cell_type": "markdown", "id": "5432fbb6", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "-\\nabla^2 u(\\boldsymbol{x}) + \\alpha u(\\boldsymbol{x}) = f(\\boldsymbol{x}) \\quad \\text{for }\\, \\boldsymbol{x}=(x, y) \\in \\Omega, \\label{eq:helmholtz} \\tag{1}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "74c368e2", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "u =0 \\text{ on } \\partial \\Omega,\n", "\\label{_auto1} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "28ce39e0", "metadata": { "editable": true }, "source": [ "where $u(\\boldsymbol{x})$ is the solution, $f(\\boldsymbol{x})$ is a function and $\\alpha$ a constant.\n", "The domain is a circular disc $\\Omega = \\{(x, y): x^2+y^2 < a^2\\}$ with radius $a$.\n", "We use polar coordinates $(\\theta, r)$, defined as" ] }, { "cell_type": "markdown", "id": "abcc039f", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " x = r \\cos \\theta, \n", "\\label{_auto2} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "efcb355f", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", " y = r \\sin \\theta,\n", "\\label{_auto3} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "ec76706d", "metadata": { "editable": true }, "source": [ "which leads to a Cartesian product mesh $(\\theta, r) \\in [0, 2\\pi) \\times (0, a)$\n", "suitable for numerical implementations. Note that the\n", "two directions are ordered with $\\theta$ first and then $r$, which is less common\n", "than $(r, \\theta)$. This has to do with the fact that we will need to\n", "solve linear equation systems along the radial direction, but not\n", "the $\\theta$-direction, since Fourier matrices are diagonal. When\n", "the radial direction is placed last, the data in the radial direction\n", "will be contigeous in a row-major C memory, leading to faster memory\n", "access where it is needed the most. Note that it takes very few\n", "changes in `shenfun` to switch the directions to $(r, \\theta)$ if this\n", "is still desired.\n", "\n", "We will use Chebyshev\n", "or Legendre basis functions $\\psi_j(r)$ for the radial direction and\n", "a periodic Fourier expansion in $\\exp(\\imath k \\theta)$ for the\n", "azimuthal direction. The polar basis functions are as such" ] }, { "cell_type": "markdown", "id": "f1e0f267", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "v_{kj}(\\theta, r) = \\exp(\\imath k \\theta) \\psi_j(r),\n", "\\label{_auto4} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "c73ec530", "metadata": { "editable": true }, "source": [ "and we look for solutions" ] }, { "cell_type": "markdown", "id": "9af694a9", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u(\\mathbf{x}) = \\tilde{u}(\\theta, r) = \\sum_{k} \\sum_{j} \\hat{u}_{kj} v_{kj}(\\theta, r).\n", "\\label{_auto5} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "7dc44a76", "metadata": { "editable": true }, "source": [ "Note that $\\tilde{u}$ is the function $u$ mapped to computational space.\n", "From now on we will simply use $u(\\theta, r)$ without the tilde, and assume that\n", "the proper version of the function is understood from its arguments.\n", "\n", "A discrete Fourier approximation space with $N$ basis functions is then" ] }, { "cell_type": "markdown", "id": "e2c0c9b8", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "V_F^N = \\text{span} \\{\\exp(\\imath k \\theta)| \\text{ for } k \\in K\\},\n", "\\label{_auto6} \\tag{7}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "b675c87a", "metadata": { "editable": true }, "source": [ "where the index set $K = \\{-N/2, -N/2+1, \\ldots, N/2-1\\}$. Since the solution $u(\\theta, r)$\n", "is real, there is Hermitian symmetry and $\\hat{u}_{k,j} = \\hat{u}_{k,-j}^*$\n", "(with $*$ denoting a complex conjugate).\n", "For this reason we use only $k \\in K=\\{0, 1, \\ldots, N/2\\}$ in solving for\n", "$\\hat{u}_{kj}$, and then use Hermitian symmetry to get the remaining\n", "unknowns. This is handled under the hood by fast Fourier transforms.\n", "\n", "The radial basis is more tricky, because there is a nontrivial 'boundary'\n", "condition (pole condition) that needs to be applied at the center of the disc $(r=0)$" ] }, { "cell_type": "markdown", "id": "a2a580a7", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial u(\\theta, 0)}{\\partial \\theta} = 0.\n", "\\label{_auto7} \\tag{8}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "29ba05d7", "metadata": { "editable": true }, "source": [ "To apply this condition we split the solution into Fourier\n", "coefficients with wavenumber 0 and $K\\backslash \\{0\\}$,\n", "remembering that the Fourier basis function with $k=0$ is\n", "simply 1" ] }, { "cell_type": "markdown", "id": "8acc3dc5", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u(\\theta, r) = \\sum_{j} \\left( \\hat{u}_{0j} \\psi_{j}(r) + \\sum_{k=1}^{N/2} \\hat{u}_{kj} \\exp(\\imath k \\theta) \\psi_j(r) \\right).\n", "\\label{_auto8} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "3712274c", "metadata": { "editable": true }, "source": [ "We then apply a different radial basis for the two $\\psi$'s in\n", "the above equation (renaming the first $\\overline{\\psi}$)" ] }, { "cell_type": "markdown", "id": "38afaa1e", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u(\\theta, r) = \\sum_{j} \\left( \\hat{u}_{0j} \\overline{\\psi}_{j}(r) + \\sum_{k=1}^{N/2} \\hat{u}_{kj} \\exp(\\imath k \\theta) \\psi_j(r) \\right).\n", "\\label{_auto9} \\tag{10}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "f444ab42", "metadata": { "editable": true }, "source": [ "Note that the first term $\\sum_{j} \\hat{u}_{0j} \\overline{\\psi}_{j}(r)$ is independent\n", "of $\\theta$. Now, to enforce conditions" ] }, { "cell_type": "markdown", "id": "eec46081", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u(\\theta, a) = 0, \n", "\\label{_auto10} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "6a63be68", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{\\partial u(\\theta, 0)}{\\partial \\theta} = 0,\n", "\\label{_auto11} \\tag{12}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "f41e1913", "metadata": { "editable": true }, "source": [ "it is sufficient for the two bases ($\\overline{\\psi}$ and $\\psi$) to\n", "satisfy" ] }, { "cell_type": "markdown", "id": "06d79350", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\overline{\\psi}_j(a) = 0, \n", "\\label{_auto12} \\tag{13}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "869ecbbc", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "\\psi_j(a) = 0,\n", "\\label{_auto13} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "21ee8cae", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "\\psi_j(0) = 0.\n", "\\label{_auto14} \\tag{15}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "4319dd50", "metadata": { "editable": true }, "source": [ "Bases that satisfy these conditions can be found both with Legendre and\n", "Chebyshev polynomials.\n", "If $\\phi_j(x)$ is used for either the Legendre polynomial $L_j(x)$ or the\n", "Chebyshev polynomial of the first kind $T_j(x)$, we can have" ] }, { "cell_type": "markdown", "id": "8227d2da", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\overline{\\psi}_j(r) = \\phi_j(2r/a-1) - \\phi_{j+1}(2r/a-1), \\text{ for } j \\in 0, 1, \\ldots N-1, \n", "\\label{_auto15} \\tag{16}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "0aac6186", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "\\psi_j(r) = \\phi_j(2r/a-1) - \\phi_{j+2}(2r/a-1), \\text{ for } j \\in 0, 1, \\ldots N-2.\n", "\\label{eq:psi} \\tag{17}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "d6161bd1", "metadata": { "editable": true }, "source": [ "Define the following approximation spaces for the radial direction" ] }, { "cell_type": "markdown", "id": "ed53c2c4", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "V_D^N = \\text{span} \\{\\psi_j\\}_{j=0}^{N-3} \n", "\\label{_auto16} \\tag{18}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "b6fdf5eb", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "V_U^N = \\text{span} \\{\\overline{\\psi}_j\\}_{j=0}^{N-2} \n", "\\label{_auto17} \\tag{19}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "976a21cf", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "\\label{_auto18} \\tag{20}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "86adfe7c", "metadata": { "editable": true }, "source": [ "and split the function space for the azimuthal direction into" ] }, { "cell_type": "markdown", "id": "6c997fcf", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "V_F^0 = \\text{span}\\{1\\}, \n", "\\label{_auto19} \\tag{21}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "e0a07851", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "V_F^{1} = \\text{span} \\{\\exp(\\imath k \\theta)\\}, \\text{ for } k \\in K \\backslash \\{0\\}.\n", "\\label{_auto20} \\tag{22}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "c171d086", "metadata": { "editable": true }, "source": [ "We then look for solutions" ] }, { "cell_type": "markdown", "id": "42be4adc", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u(\\theta, r) = u^0(r) + u^1(\\theta, r),\n", "\\label{_auto21} \\tag{23}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "ab4aeed8", "metadata": { "editable": true }, "source": [ "where" ] }, { "cell_type": "markdown", "id": "5d3b3e3d", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u^0(r) = \\sum_{j=0}^{N-2} \\hat{u}^0_j \\overline{\\psi}_j(r), \n", "\\label{_auto22} \\tag{24}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "5d6431c2", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "u^1(\\theta, r) = \\sum_{j=0}^{N-3}\\sum_{k=1}^{N/2} \\hat{u}^1_{kj} \\exp(\\imath k \\theta) \\psi_j(r) .\n", "\\label{_auto23} \\tag{25}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "6188851a", "metadata": { "editable": true }, "source": [ "As such the Helmholtz problem is split in two smaller problems.\n", "The two problems read with the spectral Galerkin method:\n", "\n", "Find $u^0 \\in V_F^0 \\otimes V_U^N$ such that" ] }, { "cell_type": "markdown", "id": "5d4eec5a", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " \\int_{\\Omega} (-\\nabla^2 u^0 + \\alpha u^0) v^0 w d\\sigma = \\int_{\\Omega} f v^0 w d\\sigma, \\quad \\forall \\, v^0 \\in V_F^0 \\otimes V_U^N.\n", "\\label{eq:u0} \\tag{26}\n", " \\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "1283888c", "metadata": { "editable": true }, "source": [ "Find $u^1 \\in V_F^1 \\otimes V_D^N$ such that" ] }, { "cell_type": "markdown", "id": "c629c307", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " \\int_{\\Omega} (-\\nabla^2 u^1 + \\alpha u^1) v^1 w d\\sigma = \\int_{\\Omega} f v^1 w d\\sigma, \\quad \\forall \\, v^1 \\in V_F^1 \\otimes V_D^N.\n", "\\label{eq:u1} \\tag{27}\n", " \\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "9d10377e", "metadata": { "editable": true }, "source": [ "Note that integration over the domain is done using\n", "polar coordinates with an integral measure of $d\\sigma=rdrd\\theta$.\n", "However, the integral in the radial direction needs to be mapped\n", "to $t=2r/a-1$, where $t \\in [-1, 1]$, which suits the basis functions used,\n", "see ([17](#eq:psi)). This leads to a measure of $0.5(t+1)adtd\\theta$.\n", "Furthermore, the weight $w(t)$ will be unity for the Legendre basis and\n", "$(1-t^2)^{-0.5}$ for the Chebyshev bases." ] }, { "cell_type": "markdown", "id": "be46f19c", "metadata": { "editable": true }, "source": [ "## Implementation\n", "\n", "\n", "A complete implementation is found in the file [unitdisc_helmholtz.py](https://github.com/spectralDNS/shenfun/blob/master/demo/unitdisc_helmholtz.py).\n", "Here we give a brief explanation for the implementation. Start by\n", "importing all functionality from [shenfun](https://github.com/spectralDNS/shenfun)\n", "and [sympy](https://sympy.org), where Sympy is required for handeling the\n", "polar coordinates. Also, we choose to work with covariant\n", "basis vectors." ] }, { "cell_type": "code", "execution_count": 1, "id": "1acc570e", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:55.315310Z", "iopub.status.busy": "2026-06-03T08:49:55.314973Z", "iopub.status.idle": "2026-06-03T08:49:55.703410Z", "shell.execute_reply": "2026-06-03T08:49:55.703142Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "from shenfun import *\n", "import sympy as sp\n", "config['basisvectors'] = 'covariant'\n", "\n", "# Define polar coordinates using angle along first axis and radius second\n", "theta, r = psi = sp.symbols('x,y', real=True, positive=True)\n", "rv = (r*sp.cos(theta), r*sp.sin(theta)) # Map to Cartesian (x, y)" ] }, { "cell_type": "markdown", "id": "355e26ca", "metadata": { "editable": true }, "source": [ "Note that Sympy symbols are both positive and real, $\\theta$ is\n", "chosen to be along the first axis and $r$ second. This has to agree with\n", "the next step, which is the creation of tensorproductspaces\n", "$V_F^0 \\otimes V_U^N$ and $V_F^1 \\otimes V_D^N$. We use\n", "`domain=(0, 1)` for the radial direction to get a unit disc, whereas\n", "the default domain for the Fourier bases is already the\n", "required $(0, 2\\pi)$." ] }, { "cell_type": "code", "execution_count": 2, "id": "3162ee49", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:55.705199Z", "iopub.status.busy": "2026-06-03T08:49:55.704975Z", "iopub.status.idle": "2026-06-03T08:49:55.873634Z", "shell.execute_reply": "2026-06-03T08:49:55.873333Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "N = 32\n", "F = FunctionSpace(N, 'F', dtype='d')\n", "F0 = FunctionSpace(1, 'F', dtype='d')\n", "L = FunctionSpace(N, 'L', bc=(0, 0), domain=(0, 1))\n", "L0 = FunctionSpace(N, 'L', bc=(None, 0), domain=(0, 1))\n", "T = TensorProductSpace(comm, (F, L), axes=(1, 0), coordinates=(psi, rv))\n", "T0 = TensorProductSpace(MPI.COMM_SELF, (F0, L0), axes=(1, 0), coordinates=(psi, rv))" ] }, { "cell_type": "markdown", "id": "6b19eba1", "metadata": { "editable": true }, "source": [ "Note that since `F0` only has one component we could actually use\n", "`L0` without creating `T0`. But the code turns out to be simpler\n", "if we use `T0`, much because the additional $\\theta$-direction is\n", "required for the polar coordinates to apply. Using one single basis\n", "function for the $\\theta$ direction is as such a generic way to handle\n", "polar 1D problems (i.e., problems that are only functions of the\n", "radial direction, but still using polar coordinates).\n", "Also note that `F` is created using the entire range of wavenumbers\n", "even though it should not include wavenumber 0.\n", "As such we need to make sure that the coefficient created for\n", "$k=0$ (i.e., $\\hat{u}^1_{0,j}$) will be exactly zero.\n", "Finally, note that\n", "`T0` is not distributed with MPI, which is accomplished using\n", "`MPI.COMM_SELF` instead of `comm` (which equals `MPI.COMM_WORLD`).\n", "The purely radial problem ([26](#eq:u0)) is only solved on the one\n", "processor with rank = 0.\n", "\n", "Polar coordinates are ensured by feeding `coordinates=(psi, rv)`\n", "to [TensorProductSpace](https://shenfun.readthedocs.io/en/latest/shenfun.html#shenfun.tensorproductspace.TensorProductSpace). Operators like [div()](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.operators.div)\n", "[grad()](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.operators.grad) and [curl()](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.operators.curl) will now work on\n", "items of [Function](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.Function), [TestFunction](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.TestFunction) and\n", "[TrialFunction](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.TrialFunction) using a polar coordinate system.\n", "\n", "To define the equations ([26](#eq:u0)) and ([27](#eq:u1)) we first declare\n", "these test- and trialfunctions, and then use code that\n", "is remarkably similar to the mathematics." ] }, { "cell_type": "code", "execution_count": 3, "id": "a9f683c4", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:55.875162Z", "iopub.status.busy": "2026-06-03T08:49:55.874995Z", "iopub.status.idle": "2026-06-03T08:49:55.953847Z", "shell.execute_reply": "2026-06-03T08:49:55.953565Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "v = TestFunction(T)\n", "u = TrialFunction(T)\n", "v0 = TestFunction(T0)\n", "u0 = TrialFunction(T0)\n", "alpha = 1\n", "\n", "mats = inner(v, -div(grad(u))+alpha*u)\n", "if comm.Get_rank() == 0:\n", " mats0 = inner(v0, -div(grad(u0))+alpha*u0)" ] }, { "cell_type": "markdown", "id": "ac2378ec", "metadata": { "editable": true }, "source": [ "Here `mats` and `mats0` will contain several tensor product\n", "matrices in the form of\n", "[TPMatrix](https://shenfun.readthedocs.io/en/latest/shenfun.html#shenfun.matrixbase.TPMatrix). Since there is only one non-periodic direction\n", "the matrices can be easily solved using [SolverGeneric1ND](https://shenfun.readthedocs.io/en/latest/shenfun.html#shenfun.la.SolverGeneric1ND).\n", "But first we need to define the function $f(\\theta, r)$.\n", "To this end we use sympy and the method of\n", "manufactured solution to define a possible solution `ue`,\n", "and then compute `f` exactly using exact differentiation" ] }, { "cell_type": "code", "execution_count": 4, "id": "56e8f67f", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:55.955165Z", "iopub.status.busy": "2026-06-03T08:49:55.955091Z", "iopub.status.idle": "2026-06-03T08:49:56.053234Z", "shell.execute_reply": "2026-06-03T08:49:56.052983Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "\n", "# Manufactured solution\n", "ue = (r*(1-r))**2*sp.cos(8*theta)-0.1*(r-1)\n", "#f = -ue.diff(r, 2) - (1/r)*ue.diff(r, 1) - (1/r**2)*ue.diff(theta, 2) + alpha*ue\n", "f = (-div(grad(u))+alpha*u).tosympy(basis=ue, psi=psi)\n", "\n", "# Compute the right hand side on the quadrature mesh\n", "fj = Array(T, buffer=f)\n", "\n", "# Take scalar product\n", "f_hat = Function(T)\n", "f_hat = inner(v, fj, output_array=f_hat)\n", "if T.local_slice(True)[0].start == 0: # The processor that owns k=0\n", " f_hat[0] = 0\n", "\n", "# For k=0 we solve only a 1D equation. Do the scalar product for Fourier\n", "# coefficient 0 by hand (or sympy).\n", "if comm.Get_rank() == 0:\n", " f0_hat = Function(T0)\n", " gt = sp.lambdify(r, sp.integrate(f, (theta, 0, 2*sp.pi))/2/sp.pi)(L0.mesh())\n", " f0_hat = T0.scalar_product(gt, f0_hat)" ] }, { "cell_type": "markdown", "id": "ebf9e1dc", "metadata": { "editable": true }, "source": [ "Note that for $u^0$ we perform the interal in the $\\theta$ direction\n", "exactly using sympy. This is necessary since one Fourier coefficient\n", "is not sufficient to do this integral numerically. For the $u^1$\n", "case we do the integral numerically as part of the [inner()](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.inner.inner) product.\n", "With the correct right hand side assembled we can solve the\n", "linear system of equations" ] }, { "cell_type": "code", "execution_count": 5, "id": "62ab606c", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:56.054642Z", "iopub.status.busy": "2026-06-03T08:49:56.054567Z", "iopub.status.idle": "2026-06-03T08:49:56.065487Z", "shell.execute_reply": "2026-06-03T08:49:56.065265Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "u_hat = Function(T)\n", "Sol1 = la.SolverGeneric1ND(mats)\n", "u_hat = Sol1(f_hat, u_hat)\n", "\n", "# case k = 0\n", "u0_hat = Function(T0)\n", "if comm.Get_rank() == 0:\n", " Sol0 = la.SolverGeneric1ND(mats0)\n", " u0_hat = Sol0(f0_hat, u0_hat)\n", "comm.Bcast(u0_hat, root=0)" ] }, { "cell_type": "markdown", "id": "8923490d", "metadata": { "editable": true }, "source": [ "Having found the solution in spectral space all that is\n", "left is to transform it back to real space." ] }, { "cell_type": "code", "execution_count": 6, "id": "410b965c", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:56.066759Z", "iopub.status.busy": "2026-06-03T08:49:56.066683Z", "iopub.status.idle": "2026-06-03T08:49:56.068855Z", "shell.execute_reply": "2026-06-03T08:49:56.068633Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "# Transform back to real space. Broadcast 1D solution\n", "sl = T.local_slice(False)\n", "uj = u_hat.backward() + u0_hat.backward()[:, sl[1]]" ] }, { "cell_type": "markdown", "id": "4c5f0cee", "metadata": { "editable": true }, "source": [ "## Postprocessing\n", "The solution can now be compared with the exact solution\n", "through" ] }, { "cell_type": "code", "execution_count": 7, "id": "e8feed8f", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:56.070091Z", "iopub.status.busy": "2026-06-03T08:49:56.070022Z", "iopub.status.idle": "2026-06-03T08:49:56.072540Z", "shell.execute_reply": "2026-06-03T08:49:56.072331Z" }, "tags": [ "thebe-init" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Error = 7.825801253188558e-15\n" ] } ], "source": [ "uq = Array(T, buffer=ue)\n", "print('Error =', np.linalg.norm(uj-uq))" ] }, { "cell_type": "markdown", "id": "fe9d2c70", "metadata": { "editable": true }, "source": [ "We can also get the gradient of the solution. For this we need\n", "a space without boundary conditions, and a vector space" ] }, { "cell_type": "code", "execution_count": 8, "id": "de99b93f", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:56.073655Z", "iopub.status.busy": "2026-06-03T08:49:56.073591Z", "iopub.status.idle": "2026-06-03T08:49:56.109553Z", "shell.execute_reply": "2026-06-03T08:49:56.109278Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "TT = T.get_orthogonal()\n", "V = VectorSpace(TT)" ] }, { "cell_type": "markdown", "id": "f090f9c5", "metadata": { "editable": true }, "source": [ "Notice that we do not have the solution in one single space\n", "in spectral space, since it is a combination of `u_hat` and\n", "`u0_hat`. For this reason we first transform the solution from\n", "real space `uj` to the new orthogonal space `TT`" ] }, { "cell_type": "code", "execution_count": 9, "id": "cf17f776", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:56.110813Z", "iopub.status.busy": "2026-06-03T08:49:56.110746Z", "iopub.status.idle": "2026-06-03T08:49:56.114709Z", "shell.execute_reply": "2026-06-03T08:49:56.114500Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "ua = Array(TT, buffer=uj)\n", "uh = ua.forward()" ] }, { "cell_type": "markdown", "id": "355a95e5", "metadata": { "editable": true }, "source": [ "With the solution as a [Function](https://shenfun.readthedocs.io/en/latest/shenfun.forms.html#shenfun.forms.arguments.Function) we can simply project\n", "the gradient to `V`" ] }, { "cell_type": "code", "execution_count": 10, "id": "79ac588b", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:56.115917Z", "iopub.status.busy": "2026-06-03T08:49:56.115843Z", "iopub.status.idle": "2026-06-03T08:49:56.153636Z", "shell.execute_reply": "2026-06-03T08:49:56.153392Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "dv = project(grad(uh), V)\n", "du = dv.backward()" ] }, { "cell_type": "markdown", "id": "24534131", "metadata": { "editable": true }, "source": [ "Note that the gradient `du` now contains the contravariant components\n", "of the covariant basis vector `b`. The basis vector `b` is not normalized\n", "(it's length is not unity), because we have set\n", "`config['basisvectors']='covariant'`. The basisvectors can\n", "be seen as" ] }, { "cell_type": "code", "execution_count": 11, "id": "9d570561", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:49:56.155005Z", "iopub.status.busy": "2026-06-03T08:49:56.154923Z", "iopub.status.idle": "2026-06-03T08:49:56.158110Z", "shell.execute_reply": "2026-06-03T08:49:56.157900Z" }, "tags": [ "thebe-init" ] }, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\mathbf{b}_{\\theta} =- r \\sin{\\left(\\theta \\right)}\\,\\mathbf{i}+r \\cos{\\left(\\theta \\right)}\\,\\mathbf{j} \\\\ \\mathbf{b}_{r} =\\cos{\\left(\\theta \\right)}\\,\\mathbf{i}+\\sin{\\left(\\theta \\right)}\\,\\mathbf{j} \\\\ $" ], "text/plain": [ "