![](\"https://cdn.jsdelivr.net/gh/spectralDNS/spectralutilities@master/figures/moebius8_trans.png\"){kind=spacer}
Figure 1
\n", "" ] }, { "cell_type": "markdown", "id": "e10f07d5", "metadata": { "editable": true }, "source": [ "## The Möbius strip\n", "\n", "A Möbius strip is the simplest non-orientable surface embedded in\n", "$\\mathbb{R}^3$. There are several realizations possible, and we\n", "will here consider the one parametrized by [[Kalvoda2020]](#Kalvoda2020)" ] }, { "cell_type": "markdown", "id": "aace7270", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " x(\\theta, t) = \\left(R-t\\cos\\frac{\\theta}{2R}\\right) \\cos \\frac{\\theta}{R}, \n", "\\label{_auto1} \\tag{1}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "1bc99e64", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", " y(\\theta, t) = \\left(R-t\\cos\\frac{\\theta}{2R}\\right) \\sin \\frac{\\theta}{R}, \n", "\\label{_auto2} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "b12752a5", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", " z(\\theta, t) = -t \\sin \\frac{\\theta}{2 R},\n", "\\label{_auto3} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "eeaf8367", "metadata": { "editable": true }, "source": [ "where $R$ is the main radius of the strip. $\\theta$ is the parameter that determines\n", "the angle for moving around the strip, like the angle of a circle.\n", "For one trip around the strip move $\\theta$ from $\\theta_0$ to $\\theta_0+2\\pi R$.\n", "A function in Cartesian coordinates $u(\\mathbf{x})$ for $\\mathbf{x} \\in \\mathbb{R}^3$\n", "is mapped to computational coordinates as $u(\\mathbf{x}) = \\tilde{u}(\\theta, t)$.\n", "By moving once around the strip a function can be seen to be twisted periodic,\n", "such that" ] }, { "cell_type": "markdown", "id": "24159786", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " \\tilde{u}(\\theta_0, t) = \\tilde{u}(\\theta_0 + 2\\pi R, -t), \n", "\\label{_auto4} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "3d69bbdd", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", " \\frac{\\partial \\tilde{u}}{\\partial \\theta}(\\theta_0, t) = \\frac{\\partial \\tilde{u}}{\\partial \\theta}(\\theta_0 + 2\\pi R, -t).\n", "\\label{_auto5} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "cdaec92b", "metadata": { "editable": true }, "source": [ "The twisted condition does not lend itself easily to a regular tensor\n", "product basis. On the other hand, by moving twice around the strip,\n", "regular periodic boundary conditions" ] }, { "cell_type": "markdown", "id": "c17b95fe", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " \\tilde{u}(\\theta_0, t) = \\tilde{u}(\\theta_0 + 4\\pi R, t),\n", "\\label{_auto6} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "63e6eb39", "metadata": { "editable": true }, "source": [ "will apply [[Kalvoda2020]](#Kalvoda2020), and we can discretize using Fourier\n", "exponentials in the $\\theta$-direction. Since the reference domain of periodic Fourier\n", "exponentials is $[-\\pi, \\pi)$ we define $\\varphi = \\theta /(2 R)$, such that\n", "the computational domain is $(\\varphi, t) \\in \\mathbb{I}^2 = [-\\pi, \\pi) \\times (-1, 1)$,\n", "with the parametrization" ] }, { "cell_type": "markdown", "id": "851f9be7", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " x(\\varphi, t) = \\left(R-t\\cos {\\varphi}\\right) \\cos {2 \\varphi}, \n", "\\label{_auto7} \\tag{7}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "60def052", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", " y(\\varphi, t) = \\left(R-t\\cos {\\varphi}\\right) \\sin {2 \\varphi}, \n", "\\label{_auto8} \\tag{8}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "ac99c073", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", " z(\\varphi, t) = -t \\sin {\\varphi}.\n", "\\label{_auto9} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "86946617", "metadata": { "editable": true }, "source": [ "Note that the reference domain corresponds to $(\\theta, t) \\in = [-2\\pi R, 2\\pi R) \\times (-1, 1)$.\n", "It is also trivial to adjust the $t$-domain with an affine map for a more narrow or\n", "wider strip, but this added complexity is not discussed here. One can simply\n", "choose the width of the strip below when choosing function space for the\n", "$t$-direction.\n", "\n", "For the $t$-direction, a mixed Legendre,\n", "$\\psi_{i} = L_{i} - L_{i+2}$, or Chebyshev, $\\psi_{i} = T_{i} - T_{i+2}$,\n", "basis can be used and we obtain a regular tensor product basis.\n", "The Legendre basis leads to more sparse matrices, and will be chosen\n", "as default, but Chebyshev also works just fine.\n", "Note that the same problem is\n", "solved by Kalvoda et al. [[Kalvoda2020]](#Kalvoda2020) using a tensor product\n", "basis with Fourier exponentials for the $\\theta$-direction and a mix\n", "of cosines and sines for the $t$-direction. Kalvoda et al. cannot\n", "make use of tensor product matrices and integrates using a\n", "two-dimensional quadrature scheme over the entire domain, leading\n", "to a dense matrix. We will here only use 1D quadrature and\n", "tensor products to get the coefficient matrix." ] }, { "cell_type": "markdown", "id": "c77d91c1", "metadata": { "editable": true }, "source": [ "## The Laplace Beltrami operator\n", "\n", "We consider the eigenvalue problem of the\n", "Laplace Beltrami operator by solving" ] }, { "cell_type": "markdown", "id": "85cf3c0c", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " -\\nabla^2 u(\\mathbf{x}) = \\lambda u(\\mathbf{x}), \\quad \\text{ for } \\mathbf{x} \\in \\Omega,\n", "\\label{_auto10} \\tag{10}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "c079b32b", "metadata": { "editable": true }, "source": [ "where $u$ is the solution, $\\lambda$ the eigenvalues and $\\Omega$ the Möbius\n", "strip. We consider only homogeneous Dirichlet boundary conditions on the boundary.\n", "\n", "To solve this problem with the spectral Galerkin method we\n", "choose an appropriate space $V$ and find $u\\in V$ such that" ] }, { "cell_type": "markdown", "id": "2e00a064", "metadata": { "editable": true }, "source": [ "$$\n", "\\int_{\\Omega} -\\nabla^2 u \\, v^* \\omega d\\sigma = \\int_{\\Omega} \\lambda u v^* \\omega d\\sigma, \\quad \\forall v \\in V\n", "$$" ] }, { "cell_type": "markdown", "id": "659c1475", "metadata": { "editable": true }, "source": [ "where $v^*$ is the complex conjugate of the test function $v$.\n", "\n", "With shenfun it is enough to operate in general coordinates as above\n", "and let the curvilinear mathematics all happen under the hood.\n", "However, for this example it is interesting to see what the\n", "Laplace-Beltrami operator looks like in computational coordinates.\n", "This is quite a bit of work by hand, and the starting point\n", "is the position vector" ] }, { "cell_type": "markdown", "id": "a9152953", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\mathbf{r} = x(\\varphi, t) \\mathbf{i} + y(\\varphi, t) \\mathbf{j} + z(\\varphi, t) \\mathbf{k}.\n", "\\label{_auto11} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "e19d947e", "metadata": { "editable": true }, "source": [ "From the position vector we compute the he covariant\n", "basis vectors $\\mathbf{b}_i = \\partial \\mathbf{r} / \\partial X^{i}$,\n", "where $\\mathbf{X} = (X^{i})_{i\\in(1,2)} = (\\varphi, t)$,\n", "and get the covariant metric tensor $g_{ij}=\\mathbf{b}_i \\cdot \\mathbf{b}_j$\n", "and its determinant $g = \\text{det}(g_{ij})$. We get" ] }, { "cell_type": "markdown", "id": "88f9eed6", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " {g} = (2R-2t\\cos \\varphi)^2+t^2.\n", "\\label{_auto12} \\tag{12}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "d478e8e2", "metadata": { "editable": true }, "source": [ "and" ] }, { "cell_type": "markdown", "id": "0382744f", "metadata": { "editable": true }, "source": [ "$$\n", "(g_{ij}) =\n", "\\begin{pmatrix}\n", " g & 0 \\\\ \n", " 0 & 1\n", "\\end{pmatrix}\n", "$$" ] }, { "cell_type": "markdown", "id": "1edbc444", "metadata": { "editable": true }, "source": [ "Likewise, the contravariant metric tensor $g^{ij}$ is the\n", "inverse of the covariant" ] }, { "cell_type": "markdown", "id": "8e3ccac6", "metadata": { "editable": true }, "source": [ "$$\n", "(g^{ij}) =\n", "\\begin{pmatrix}\n", " \\frac{1}{g} & 0 \\\\ \n", " 0 & 1\n", "\\end{pmatrix}\n", "$$" ] }, { "cell_type": "markdown", "id": "a3d85733", "metadata": { "editable": true }, "source": [ "Please note that all these metrics and other terms are computed\n", "under the hood by shenfun, and a user does not normally have to worry\n", "about these.\n", "\n", "It can be shown that the Laplace-Beltrami operator in curvilinear coordinates\n", "is given as" ] }, { "cell_type": "markdown", "id": "d2caccd3", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\nabla^2 \\tilde{u} = \\frac{1}{\\sqrt{g}}\\frac{\\partial}{\\partial X^{i}}\\left( g^{ij}\\sqrt{g} \\frac{\\partial \\tilde{u}}{\\partial X^{j}}\\right),\n", "\\label{_auto13} \\tag{13}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "856b5034", "metadata": { "editable": true }, "source": [ "with summation on repeated indices.\n", "Using this and the surface element $d\\sigma=\\sqrt{g} d\\varphi dt$, the\n", "variational form in computational coordinates becomes" ] }, { "cell_type": "markdown", "id": "fbdc9605", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " -\\int_{\\mathbb{I}^2}\\frac{\\partial}{\\partial X^{i}}\\left( g^{ij}\\sqrt{g} \\frac{\\partial \\tilde{u}}{\\partial X^{j}}\\right) \\, \\tilde{v}^*\\, \\tilde{\\omega} d\\varphi dt = \\int_{\\mathbb{I}^2} \\lambda \\tilde{u} \\tilde{v}^*\\, \\tilde{\\omega}\\sqrt{g} d\\varphi dt,\n", "\\label{_auto14} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "346133c3", "metadata": { "editable": true }, "source": [ "Summing and expanding some derivatives, we get" ] }, { "cell_type": "markdown", "id": "2fa380d8", "metadata": { "editable": true }, "source": [ "$$\n", "\\begin{align*}\n", " -\\int_{\\mathbb{I}^2} \\Big\\{ \\frac{1}{\\sqrt{g}}\\frac{\\partial^2 \\tilde{u}}{\\partial \\varphi^2 } &+\\frac{4 t \\sin \\varphi \\left(t \\cos {\\varphi} -R \\right) }{g^{\\frac{3}{2}}}\\frac{\\partial \\tilde{u}}{\\partial \\varphi } \\notag \\\\ \n", " & +\\sqrt{g}\\frac{\\partial^2 \\tilde{u}}{\\partial t^2 }\n", " +\\frac{4 \\cos {\\varphi} (t \\cos \\varphi -R) + t}{ \\sqrt{g}}\\frac{\\partial \\tilde{u}}{\\partial t } \\Big\\} \\tilde{v}^* \\tilde{\\omega} d\\varphi dt \\notag \\\\ \n", " &= \\int_{\\mathbb{I}^2} \\lambda \\tilde{u} \\tilde{v}^* \\tilde{\\omega} \\sqrt{g} d\\varphi dt,\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "id": "8eb40114", "metadata": { "editable": true }, "source": [ "Note that the variational problem contains unseparable variable coefficients,\n", "like $1/\\sqrt{g}$ and $\\sqrt{g}$ and as such cannot\n", "easily be solved using efficient tensor product algebra.\n", "However, note that both $g$ and $g^2$ are separable (i.e., they can\n", "be written as products of simpler functions separated by the arguments\n", "$g(\\varphi, t) = \\sum_k g_1^k(\\varphi)g_2^k(t)$),\n", "and using the unconventional weight $\\tilde{\\omega} = g^{3/2}$ we get" ] }, { "cell_type": "markdown", "id": "c2f7c6b4", "metadata": { "editable": true }, "source": [ "$$\n", "\\begin{align*}\n", " -\\int_{\\mathbb{I}^2} \\Big\\{ g \\frac{\\partial^2 \\tilde{u}}{\\partial \\varphi^2 } &+ {4 t \\sin \\varphi \\left(t \\cos {\\varphi} -R \\right) } \\frac{\\partial \\tilde{u}}{\\partial \\varphi } \\notag \\\\ \n", " & + g^2 \\frac{\\partial^2 \\tilde{u}}{\\partial t^2 }\n", " +g \\left({4 \\cos {\\varphi} (t \\cos \\varphi -R) + t} \\right) \\frac{\\partial \\tilde{u}}{\\partial t } \\Big\\} \\tilde{v}^* d\\varphi dt \\notag \\\\ \n", " &= \\int_{\\mathbb{I}^2} \\lambda g^{2} \\tilde{u} \\tilde{v}^* d\\varphi dt,\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "id": "1c9b6780", "metadata": { "editable": true }, "source": [ "where the left hand side can be assembled using tensor product matrices,\n", "where no single 1D matrix has more than 9 diagonals.\n", "\n", "Fortunately we do not have to do all this by hand since we have a\n", "software that automatically assembles such matrices for us.\n", "Now let's see how this problem can be handled with shenfun." ] }, { "cell_type": "markdown", "id": "13451900", "metadata": { "editable": true }, "source": [ "## Implementation\n", "\n", "First we need to create function spaces for each direction $\\varphi$\n", "and $t$, and then a tensor product space from the two. We use the\n", "parametrization given above and shenfun will then automatically\n", "differentiate to create basis functions and metrics, like $g$.\n", "One important factor, though. Sympy's [simplify](https://docs.sympy.org/latest/modules/simplify/simplify.html)\n", "will sometimes have problems finding the best possible simplification\n", "of a term, like $g$. And in this particular case we need\n", "to discourage the use of powers, or else sympy will end up with a\n", "much more complicated $g$ than we get below." ] }, { "cell_type": "code", "execution_count": 1, "id": "e3a23ab5", "metadata": { "collapsed": false, "editable": true, "execution": { "iopub.execute_input": "2024-02-19T13:48:13.558879Z", "iopub.status.busy": "2024-02-19T13:48:13.558593Z", "iopub.status.idle": "2024-02-19T13:48:16.662340Z", "shell.execute_reply": "2024-02-19T13:48:16.661624Z" } }, "outputs": [], "source": [ "from shenfun import *\n", "import sympy as sp\n", "from IPython.display import Math, Latex, display\n", "from scipy.sparse.linalg import eigs\n", "config['basisvectors'] = 'covariant'\n", "\n", "phi, t = psi = sp.symbols('x,y', real=True)\n", "\n", "RR = sp.Rational(132, 20)/sp.pi # Same as Kalvoda et al\n", "#RR = 2\n", "# Use a symbolic R first, then later substitute for the value in RR\n", "R = sp.Symbol('R', real=True, positive=True)\n", "rv = ((R-t*sp.cos(phi))*sp.cos(2*phi),\n", " (R-t*sp.cos(phi))*sp.sin(2*phi),\n", " -t*sp.sin(phi))\n", "\n", "def discourage_powers(expr):\n", " POW = sp.Symbol('POW')\n", " count = sp.count_ops(expr, visual=True)\n", " count = count.replace(POW, 100)\n", " count = count.replace(sp.Symbol, type(sp.S.One))\n", " return count\n", "\n", "N = (80, 40)\n", "B0 = FunctionSpace(N[0], 'F', domain=(-np.pi, np.pi), dtype='D')\n", "B1 = FunctionSpace(N[1], 'L', bc=(0, 0), domain=(-0.75, 0.75)) # Use same domain as Kalvoda et al\n", "T = TensorProductSpace(comm, (B0, B1), coordinates=(psi, rv, True, (), discourage_powers), axes=(1, 0))\n", "\n", "u = TrialFunction(T)\n", "v = TestFunction(T)" ] }, { "cell_type": "markdown", "id": "0aeb1ce8", "metadata": { "editable": true }, "source": [ "Note the `discourage_powers` function. Try to turn it off (if\n", "you are watching this in an interactive setting) and see what happens\n", "to, e.g., `T.coors.sg`, which is $\\sqrt{g}$.\n", "\n", "We can now check to see what the Laplace-Beltrami operator\n", "looks like when computed with shenfun:" ] }, { "cell_type": "code", "execution_count": 2, "id": "64eca61f", "metadata": { "collapsed": false, "editable": true, "execution": { "iopub.execute_input": "2024-02-19T13:48:16.665623Z", "iopub.status.busy": "2024-02-19T13:48:16.665219Z", "iopub.status.idle": "2024-02-19T13:48:18.645317Z", "shell.execute_reply": "2024-02-19T13:48:18.644836Z" } }, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\frac{1}{g{\\left(\\varphi,t \\right)}}\\frac{\\partial^2 f}{\\partial \\varphi^2 }- \\frac{4 t \\left(R - t \\cos{\\left(\\varphi \\right)}\\right) \\sin{\\left(\\varphi \\right)}}{g^{2}{\\left(\\varphi,t \\right)}}\\frac{\\partial f}{\\partial \\varphi }+\\frac{\\partial^2 f}{\\partial t^2 }+\\frac{- 4 R \\cos{\\left(\\varphi \\right)} + 4 t \\cos^{2}{\\left(\\varphi \\right)} + t}{g{\\left(\\varphi,t \\right)}}\\frac{\\partial f}{\\partial t }$" ], "text/plain": [ "