{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Curvature quantities on Riemannian manifolds\n", "\n", "In this notebook, we discuss the definition and computation of different curvature quantities of Riemannian manifolds and their relations, see [this notebook](riemannian_manifolds.ipynb) for a brief introduction to Riemannian manifolds. We focus on the Riemann curvature tensor, curvature operator, Ricci curvature tensor, scalar curvature, and (in 3D) the Einstein tensor.\n", "\n", "We use two- and three-dimensional numerical examples to verify the identities. We consider the Poincare disc in 2D and the three-sphere in 3D. For the numerical computations, we interpolate the exact metric into the [Regge finite elements](regge_metric.ipynb)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from ngsolve import *\n", "from ngsolve.webgui import Draw\n", "from netgen.occ import *\n", "import ngsdiffgeo as dg\n", "from ngsolve.fem import Einsum, LeviCivitaSymbol\n", "\n", "mesh2d = Mesh(\n", " OCCGeometry(Circle((0, 0), 0.7).Face(), dim=2).GenerateMesh(maxh=0.05)\n", ").Curve(3)\n", "pcd = dg.PoincareDisk()\n", "Draw(Norm(pcd.metric), mesh2d)\n", "\n", "eps = 2e-1\n", "mesh3d = Mesh(\n", " OCCGeometry(Box((eps, eps, eps), (pi - eps, pi - eps, 2 * pi - eps))).GenerateMesh(\n", " maxh=0.3\n", " )\n", ")\n", "sphere3 = dg.Sphere3()\n", "Draw(Norm(sphere3.metric), mesh3d, order=3)\n", "\n", "\n", "# Euclidean Levi-Civita permuting symbols in 2D and 3D\n", "Eps2d = LeviCivitaSymbol(2)\n", "Eps3d = LeviCivitaSymbol(3)\n", "\n", "# Interpolate exact metric to Regge space\n", "g2d = GridFunction(HCurlCurl(mesh2d, order=5))\n", "g3d = GridFunction(HCurlCurl(mesh3d, order=4))\n", "with TaskManager():\n", " g2d.Set(pcd.metric, dual=True)\n", " g3d.Set(sphere3.metric, dual=True)\n", "\n", "# Define Christoffel symbols of first and second kind\n", "chr1_2d = g2d.Operator(\"christoffel\")\n", "chr2_2d = g2d.Operator(\"christoffel2\")\n", "chr1_3d = g3d.Operator(\"christoffel\")\n", "chr2_3d = g3d.Operator(\"christoffel2\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Riemann curvature tensor\n", "We start with the fourth-order Riemann curvature tensor $\\mathfrak{R}$. Its coordinate expression is\n", "\\begin{align*}\n", "\\mathfrak{R}_{ijkl} = \\partial_i \\Gamma_{jkl}-\\partial_j \\Gamma_{ikl} + \\Gamma_{ik}^p\\Gamma_{jlp}- \\Gamma_{jk}^p\\Gamma_{ilp}.\n", "\\end{align*}\n", "\n", "First, we note that the first two terms can be rewritten using the incompatibility operator $\\mathrm{inc}(g)= \\varepsilon^{ij}\\varepsilon^{kl}\\partial_i\\partial_k g_{jl}$ in 2D and $\\mathrm{inc}(g)^{ij}= \\varepsilon^{ikl}\\varepsilon^{jmn}\\partial_k\\partial_m g_{ln}$ in 3D. There holds\n", "\\begin{align*}\n", "\\partial_i \\Gamma_{jkl}-\\partial_j \\Gamma_{ikl} = \\begin{cases}\\frac{1}{2}\\varepsilon_{ij}\\varepsilon_{kl}\\mathrm{inc}(g), & \\text{in 2D},\\\\\n", "\\frac{1}{2}\\varepsilon_{ija}\\varepsilon_{klb}\\mathrm{inc}(g)^{ab}, & \\text{in 3D}.\n", "\\end{cases} \n", "\\end{align*} " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We compare the Riemann differential operator from the Regge finite element with its coordinate expression from the metric $g$ and the exact Riemann curvature tensor. Note that the numerical computations fit up to numerical rounding error precision. Increasing the resolution of the mesh reduces the relative error between the numerical and exact curvature. We will learn later how to approximate the curvature more accurately from the metric tensor." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# NGSolve Riemann curvature tensor\n", "Riem1_2d = g2d.Operator(\"Riemann\")\n", "# manually compute Riemann curvature tensor\n", "Riem2_2d = (\n", " 0.5 * g2d.Operator(\"inc\") * Einsum(\"ij,kl->ijkl\", Eps2d, Eps2d)\n", " + Einsum(\"ikp,jlp->ijkl\", chr2_2d, chr1_2d)\n", " - Einsum(\"jkp,ilp->ijkl\", chr2_2d, chr1_2d)\n", ")\n", "# exact Riemann curvature tensor\n", "Riem_ex_2d = pcd.Riemann\n", "\n", "# The same in 3D\n", "Riem1_3d = g3d.Operator(\"Riemann\")\n", "Riem2_3d = (\n", " 0.5 * Einsum(\"ija,klb,ab->ijkl\", Eps3d, Eps3d, g3d.Operator(\"inc\"))\n", " + Einsum(\"ikp,jlp->ijkl\", chr2_3d, chr1_3d)\n", " - Einsum(\"jkp,ilp->ijkl\", chr2_3d, chr1_3d)\n", ")\n", "Riem_ex_3d = sphere3.Riemann\n", "\n", "# check that the expressions are equal\n", "with TaskManager():\n", " print(\n", " f\"err Riemann curvature tensor 2d = {sqrt(Integrate(InnerProduct(Riem1_2d - Riem2_2d,Riem1_2d - Riem2_2d),mesh2d)):.8f}\"\n", " )\n", " print(\n", " f\"rel err Riemann curvature tensor 2d exact = {sqrt(Integrate(InnerProduct(Riem1_2d - Riem_ex_2d,Riem1_2d - Riem_ex_2d),mesh2d))/sqrt(Integrate(InnerProduct(Riem_ex_2d, Riem_ex_2d),mesh2d)):.8f}\"\n", " )\n", " print(\n", " f\"err Riemann curvature tensor 3d = {sqrt(Integrate(InnerProduct(Riem1_3d - Riem2_3d,Riem1_3d - Riem2_3d),mesh3d)):.8f}\"\n", " )\n", "\n", " print(\n", " f\"rel err Riemann curvature tensor 3d exact = {sqrt(Integrate(InnerProduct(Riem1_3d - Riem_ex_3d,Riem1_3d - Riem_ex_3d),mesh3d))/sqrt(Integrate(InnerProduct(Riem_ex_3d, Riem_ex_3d),mesh3d)):.8f}\"\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Curvature operator\n", "Next, we look at the curvature operator $\\mathcal{Q}$, which incorporates the skew-symmetry property of the Riemann curvature tensor. It is defined via the identity\n", "\\begin{align*}\n", "\\mathcal{Q}(X^\\flat\\wedge Y^\\flat, W^\\flat \\wedge Z^\\flat) = \\mathfrak{R}(X,Y,Z,W),\\qquad \\forall X,Y,Z,W\\in\\mathfrak{X}(M).\n", "\\end{align*}\n", "Note the convention that the $Z$ and $W$ arguments are swapped and that e.g. $X^\\flat\\wedge Y^\\flat \\in \\Lambda^2(M)$ is a bi-vector. In two dimensions the curvature operator can be identified with a scalar function as $\\mathcal{Q}\\,\\omega(X,Y)\\,\\omega(W,Z) \\simeq \\mathcal{Q}(X^\\flat\\wedge Y^\\flat, W^\\flat\\wedge Z^\\flat)$ because the bi-vectors and the volume form are linearly dependent. Thus, in coordinates, there holds\n", "\\begin{align*}\n", "\\mathcal{Q} = -\\frac{1}{4}\\hat{\\varepsilon}^{kl}\\hat{\\varepsilon}^{mn}\\mathfrak{R}_{klmn} = -\\frac{1}{4\\,\\det g}\\varepsilon^{kl}\\varepsilon^{mn}\\mathfrak{R}_{klmn},\n", "\\end{align*}\n", "where $\\hat{\\varepsilon}^{kl}= \\frac{1}{\\sqrt{\\det g}}\\varepsilon^{kl}$ denotes the Levi-Civita symbol with respect to the metric. In three dimensions, $X^\\flat\\wedge Y^\\flat$ can be identified with the cross product $(X\\times Y)^\\flat$. The curvature operator becomes a symmetric operator acting on two 1-forms and reads in coordinates\n", "\\begin{align*}\n", "\\mathcal{Q}^{ij} = -\\frac{1}{4}\\hat{\\varepsilon}^{ikl}\\hat{\\varepsilon}^{jmn}\\mathfrak{R}_{klmn} = -\\frac{1}{4\\,\\det g}\\varepsilon^{ikl}\\varepsilon^{jmn}\\mathfrak{R}_{klmn}.\n", "\\end{align*}\n", "Using the coordinate expression of the Riemann curvature tensor and using the identity $\\hat{\\varepsilon}^{ija}\\hat{\\varepsilon}_{kla}={\\varepsilon}^{ija}{\\varepsilon}_{kla}= \\delta^i_k\\delta^j_l-\\delta^i_l\\delta^j_k$ we have after a short calculation in 2D and 3D\n", "\\begin{align*}\n", "\\mathcal{Q}=-\\frac{1}{2}\\hat{\\varepsilon}^{kl}\\hat{\\varepsilon}^{mn}\\left(\\partial_k\\partial_m g_{ln}+\\Gamma^q_{km}\\Gamma_{lnq}\\right),\\qquad \\mathcal{Q}^{ij}=-\\frac{1}{2}\\hat{\\varepsilon}^{ikl}\\hat{\\varepsilon}^{jmn}\\left(\\partial_k\\partial_m g_{ln}+\\Gamma^q_{km}\\Gamma_{lnq}\\right).\n", "\\end{align*}\n", "\n", "On the opposite, we can express the Riemann curvature tensor in terms of the curvature operator\n", "\\begin{align*}\n", "\\mathfrak{R}_{ijkl}=\\begin{cases} \\hat{\\varepsilon}_{ij}\\hat{\\varepsilon}_{kl}\\mathcal{Q}, &\\text{in 2D}\\\\\n", "\\hat{\\varepsilon}_{ija}\\hat{\\varepsilon}_{klb}\\mathcal{Q}^{ab}&\\text{in 3D}.\\end{cases}\n", "\\end{align*}\n", "We recall that $\\hat{\\varepsilon}_{ij} = \\sqrt{\\det g}\\varepsilon_{ij}$.\n", "\n", "WARNING: The curvature differential operator .Operator(\"curvature\") of a Regge object returns $\\det g\\,\\mathcal{Q}$ instead of $\\mathcal{Q}$! Thus, we need to divide by the determinant." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Q1_2d = 1 / Det(g2d) * g2d.Operator(\"curvature\")\n", "Q2_2d = (\n", " -1 / (4 * Det(g2d)) * Einsum(\"ij,kl,ijkl\", Eps2d, Eps2d, g2d.Operator(\"Riemann\"))\n", ")\n", "Q3_2d = (\n", " -0.5\n", " / Det(g2d)\n", " * (g2d.Operator(\"inc\") + Einsum(\"kl,mn,kmq,lnq\", Eps2d, Eps2d, chr2_2d, chr1_2d))\n", ")\n", "Q_ex_2d = pcd.curvature\n", "\n", "Q1_3d = 1 / Det(g3d) * g3d.Operator(\"curvature\")\n", "Q2_3d = (\n", " -1\n", " / (4 * Det(g3d))\n", " * Einsum(\"ikl,jmn,klmn->ij\", Eps3d, Eps3d, g3d.Operator(\"Riemann\"))\n", ")\n", "Q3_3d = (\n", " -0.5\n", " / Det(g3d)\n", " * (\n", " g3d.Operator(\"inc\")\n", " + Einsum(\"ikl,jmn,kmq,lnq->ij\", Eps3d, Eps3d, chr2_3d, chr1_3d)\n", " )\n", ")\n", "Q_ex_3d = sphere3.curvature\n", "\n", "with TaskManager():\n", " print(\n", " f\"err curvature operator 2d 1-2 = {sqrt(Integrate((Q1_2d - Q2_2d)**2,mesh2d)):.8f}\"\n", " )\n", " print(\n", " f\"err curvature operator 2d 1-3 = {sqrt(Integrate((Q1_2d - Q3_2d)**2,mesh2d)):.8f}\"\n", " )\n", " print(\n", " f\"rel err curvature operator 2d exact = {sqrt(Integrate((Q1_2d - Q_ex_2d)**2,mesh2d))/sqrt(Integrate(Q_ex_2d**2,mesh2d)):.8f}\"\n", " )\n", "\n", " print(\n", " f\"err curvature operator 3d 1-2 = {sqrt(Integrate(InnerProduct(Q1_3d - Q2_3d,Q1_3d - Q2_3d),mesh3d)):.8f}\"\n", " )\n", " print(\n", " f\"err curvature operator 3d 1-3 = {sqrt(Integrate(InnerProduct(Q1_3d - Q3_3d,Q1_3d - Q3_3d),mesh3d)):.8f}\"\n", " )\n", " print(\n", " f\"rel err curvature operator 3d exact = {sqrt(Integrate(InnerProduct(Q1_3d - Q_ex_3d,Q1_3d - Q_ex_3d),mesh3d))/sqrt(Integrate(InnerProduct(Q_ex_3d, Q_ex_3d),mesh3d)):.8f}\"\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ricci curvature tensor\n", "The Ricci curvature tensor is defined by contracting two of the indices of the Riemann curvature tensor (i.e., taking a trace)\n", "\\begin{align*}\n", "\\mathrm{Ric}_{ij} = \\mathfrak{R}_{kijl}g^{kl},\n", "\\end{align*}\n", "In terms of the curvature operator, we have\n", "\\begin{align*}\n", "\\mathrm{Ric} = \\begin{cases} \\mathcal{Q}\\,g, &\\text{in 2D}\\\\ g^{-1}\\times \\mathcal{Q}, &\\text{in 3D}\\end{cases},\n", "\\end{align*}\n", "where $(A \\times B)_{ij}= \\hat{\\varepsilon}_{ikl}\\hat{\\varepsilon}_{jmn}A^{km}B^{ln}$ is the tensor cross product. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Ric1_2d = g2d.Operator(\"Ricci\")\n", "Ric2_2d = 1 / Det(g2d) * g2d.Operator(\"curvature\") * g2d\n", "Ric_ex_2d = pcd.Ricci\n", "\n", "Ric1_3d = g3d.Operator(\"Ricci\")\n", "Ric2_3d = Det(g3d) * Einsum(\n", " \"ikl,jmn,km,ln->ij\",\n", " Eps3d,\n", " Eps3d,\n", " Inv(g3d),\n", " 1 / Det(g3d) * g3d.Operator(\"curvature\"),\n", ")\n", "Ric_ex_3d = sphere3.Ricci\n", "\n", "with TaskManager():\n", " print(\n", " f\"err Ricci curvature tensor 2d = {sqrt(Integrate(InnerProduct(Ric1_2d - Ric2_2d,Ric1_2d - Ric2_2d),mesh2d)):.8f}\"\n", " )\n", " print(\n", " f\"rel err Ricci curvature tensor 2d exact = {sqrt(Integrate(InnerProduct(Ric1_2d - Ric_ex_2d,Ric1_2d - Ric_ex_2d),mesh2d))/sqrt(Integrate(InnerProduct(Ric_ex_2d, Ric_ex_2d),mesh2d)):.8f}\"\n", " )\n", " print(\n", " f\"err Ricci curvature tensor 3d = {sqrt(Integrate(InnerProduct(Ric1_3d - Ric2_3d,Ric1_3d - Ric2_3d),mesh3d)):.8f}\"\n", " )\n", " print(\n", " f\"rel err Ricci curvature tensor 3d exact = {sqrt(Integrate(InnerProduct(Ric1_3d - Ric_ex_3d,Ric1_3d - Ric_ex_3d),mesh3d))/sqrt(Integrate(InnerProduct(Ric_ex_3d, Ric_ex_3d),mesh3d)):.8f}\"\n", " )\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Scalar curvature\n", "The scalar curvature is given by contracting the Riemann curvature tensor twice (or equivalently taking the trace of the Ricci curvature tensor)\n", "\\begin{align*}\n", "S = \\mathrm{Ric}_{ij} g^{ij} = \\mathfrak{R}_{kijl}g^{kl}g^{ij}.\n", "\\end{align*}\n", "\n", "In two dimensions, the Gauss curvature is given for linear independent $X,Y$ by\n", "\\begin{align*}\n", "K = \\frac{\\mathcal{R}(X,Y,Y,X)}{g(X,X)g(Y,Y)-g(X,Y)^2},\n", "\\end{align*}\n", "which is independent of the choice of $X,Y$. Further, it is related to the scalar curvature by $S=2K$.\n", "\n", "In two dimensions, the Gauss curvature $K$ coincides with the curvature operator $\\mathcal{Q}$, $K=\\mathcal{Q}$. In 3D, the scalar curvature can be computed from the curvature operator by taking its trace\n", "\\begin{align*}\n", "S = 2 \\mathrm{tr}_g(\\mathcal{Q}) = 2 g_{ij}\\mathcal{Q}^{ij}.\n", "\\end{align*}" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "S1_2d = g2d.Operator(\"scalar\")\n", "S2_2d = Einsum(\"il,jk,ijkl\", Inv(g2d), Inv(g2d), g2d.Operator(\"Riemann\"))\n", "S3_2d = 2 / Det(g2d) * g2d.Operator(\"curvature\")\n", "S_ex_2d = pcd.scalar\n", "\n", "S1_3d = g3d.Operator(\"scalar\")\n", "S2_3d = Einsum(\"il,jk,ijkl\", Inv(g3d), Inv(g3d), g3d.Operator(\"Riemann\"))\n", "S3_3d = 2 / Det(g3d) * InnerProduct(g3d, g3d.Operator(\"curvature\"))\n", "S_ex_3d = sphere3.scalar\n", "\n", "with TaskManager():\n", " print(\n", " f\"err scalar curvature 2d 1-2 = {sqrt(Integrate((S1_2d - S2_2d)**2,mesh2d)):.8f}\"\n", " )\n", " print(\n", " f\"err scalar curvature 2d 1-3 = {sqrt(Integrate((S1_2d - S3_2d)**2,mesh2d)):.8f}\"\n", " )\n", " print(\n", " f\"rel err scalar curvature 2d exact = {sqrt(Integrate((S1_2d - S_ex_2d)**2,mesh2d))/sqrt(Integrate(S_ex_2d**2,mesh2d)):.8f}\"\n", " )\n", "\n", " print(\n", " f\"err scalar curvature 3d 1-2 = {sqrt(Integrate((S1_3d - S2_3d)**2,mesh3d)):.8f}\"\n", " )\n", " print(\n", " f\"err scalar curvature 3d 1-3 = {sqrt(Integrate((S1_3d - S3_3d)**2,mesh3d)):.8f}\"\n", " )\n", " print(\n", " f\"rel err scalar curvature 3d exact = {sqrt(Integrate((S1_3d - S_ex_3d)**2,mesh3d))/sqrt(Integrate(S_ex_3d**2,mesh3d)):.8f}\"\n", " )\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Einstein tensor\n", "The Einstein tensor is defined by\n", "\\begin{align*}\n", "G_{ij}=\\mathrm{Ric}_{ij}-\\frac{1}{2}S\\,g_{ij},\\qquad G = J\\mathrm{Ric},\n", "\\end{align*}\n", "where $J\\sigma = \\sigma-\\frac{1}{2}\\mathrm{tr}(\\sigma)\\,g$. In three dimensions, it is the negative of the double lowering of the curvature operator \n", "\\begin{align*}\n", "G = -\\mathcal{Q}^{\\flat\\flat},\\qquad G_{ij} = -g_{ik}\\mathcal{Q}^{kl}g_{lj}.\n", "\\end{align*}" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "G1_3d = g3d.Operator(\"Einstein\")\n", "G2_3d = g3d.Operator(\"Ricci\") - 0.5 * g3d.Operator(\"scalar\") * g3d\n", "G3_3d = -g3d * (1 / Det(g3d) * g3d.Operator(\"curvature\")) * g3d\n", "G_ex_3d = sphere3.Einstein\n", "\n", "with TaskManager():\n", " print(\n", " f\"err Einstein tensor 3d 1-2 = {sqrt(Integrate(InnerProduct(G1_3d - G2_3d,G1_3d - G2_3d),mesh3d)):.8f}\"\n", " )\n", " print(\n", " f\"err Einstein tensor 3d 1-3 = {sqrt(Integrate(InnerProduct(G1_3d - G3_3d,G1_3d - G3_3d),mesh3d)):.8f}\"\n", " )\n", " print(\n", " f\"rel err Einstein tensor 3d exact = {sqrt(Integrate(InnerProduct(G1_3d - G_ex_3d,G1_3d - G_ex_3d),mesh3d))/sqrt(Integrate(InnerProduct(G_ex_3d, G_ex_3d),mesh3d)):.8f}\"\n", " )\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.13.1" } }, "nbformat": 4, "nbformat_minor": 2 }