Exterior covariant derivatives of double forms

This notebook continues Double forms. We introduce the exterior covariant derivative and coderivative in either slot, and verify their basic algebraic and integration-by-parts identities.

\[\def\owedge{\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}}} \def\Tr{\mathrm{Tr}}\]

[1]:

import ngsolve
from ngsolve import (
    BND,
    VOL,
    CF,
    Id,
    InnerProduct,
    Integrate,
    Mesh,
    TaskManager,
    sqrt,
    x,
    y,
    z,
)
from netgen.occ import unit_cube
import ngsdiffgeo as dg

TOL = 1e-7
mesh = Mesh(unit_cube.GenerateMesh(maxh=2))
mf = dg.RiemannianManifold(dg.Heisenberg().metric, normal_sign=-1)


def l2_error(a, b, vb=VOL, bonus_intorder=4):
    dX = (
        ngsolve.dx(element_boundary=True, bonus_intorder=bonus_intorder)
        if vb == BND
        else ngsolve.dx(bonus_intorder=bonus_intorder)
    )
    return sqrt(Integrate(InnerProduct(a - b, a - b) * dX, mesh))

We use smooth double forms on a three-dimensional non-Euclidean manifold. A double form of bidegree \((p,q)\) is represented by dg.DoubleForm(..., p=p, q=q, dim=3).

[2]:

alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))
beta = dg.OneForm(CF((0.3 * z * y, x * z**2, y**2)))
gamma = dg.OneForm(CF((y * z - x, x**2 + z, x * y)))
eta = dg.OneForm(CF((z + x, y - z, x * y)))

beta_gamma = dg.Wedge(beta, gamma)
alpha_beta = dg.Wedge(alpha, beta)

A11 = dg.DoubleForm(dg.Einsum("i,j->ij", alpha, beta), p=1, q=1, dim=3)
B11 = dg.DoubleForm(dg.Einsum("i,j->ij", gamma, eta), p=1, q=1, dim=3)
A12 = dg.DoubleForm(dg.Einsum("i,jk->ijk", alpha, beta_gamma), p=1, q=2, dim=3)
A21 = dg.DoubleForm(dg.Einsum("ij,k->ijk", alpha_beta, gamma), p=2, q=1, dim=3)
A22 = dg.Wedge(A11, B11)
X = dg.VectorField(CF((x - 0.1 * y * z, y + 0.2 * z * x, z - 0.3 * x * y)))

Left and right slots

For \(\varphi\in\Lambda^{p,q}(\Omega)\), the operators mf.d_cov and mf.delta_cov act on the left slot by default. Set slot="right" to act on the right slot. Thus,

\[d^\nabla:\Lambda^{p,q}\to\Lambda^{p+1,q},\qquad d^{\nabla,\prime}:\Lambda^{p,q}\to\Lambda^{p,q+1},\]

and the two coderivatives lower the corresponding degree. Transposition exchanges the slots:

\[d^{\nabla,\prime}\varphi=\bigl(d^\nabla(\varphi^T)\bigr)^T,\qquad \delta^{\nabla,\prime}\varphi=\bigl(\delta^\nabla(\varphi^T)\bigr)^T.\]

[3]:

with TaskManager():
    for phi in (A11, A12, A21, A22):
        d_right = mf.d_cov(phi, slot="right")
        d_right_from_transpose = mf.d_cov(phi.trans, slot="left").trans
        delta_right = mf.delta_cov(phi, slot="right")
        delta_right_from_transpose = mf.delta_cov(phi.trans, slot="left").trans

        assert l2_error(d_right, d_right_from_transpose) < TOL
        assert l2_error(delta_right, delta_right_from_transpose) < TOL

Covariant divergence and graph compilation

mf.CovDiv contracts the covariant derivative in the selected slot. On a double form, the corresponding covariant coderivative is its negative:

\[\delta^\nabla\varphi=-\operatorname{CovDiv}_{L}\varphi,\qquad \delta^{\nabla,\prime}\varphi=-\operatorname{CovDiv}_{R}\varphi.\]

For a complicated coefficient expression that is evaluated repeatedly, compile_inner="graph" can compile the input graph shared by the gradient and connection terms of d_cov and delta_cov. This is an opt-in optimization for coefficient expressions. Trial and test proxy expressions should use the default mode.

[4]:

with TaskManager():
    for slot in ("left", "right"):
        assert l2_error(mf.CovDiv(A11, slot=slot), -mf.delta_cov(A11, slot=slot)) < TOL
        assert (
            l2_error(
                mf.d_cov(A11, slot=slot, compile_inner="graph"),
                mf.d_cov(A11, slot=slot),
            )
            < TOL
        )
        assert (
            l2_error(
                mf.delta_cov(A11, slot=slot, compile_inner="graph"),
                mf.delta_cov(A11, slot=slot),
            )
            < TOL
        )

Leibniz rules

The exterior covariant derivative is a graded derivation in the selected slot. For \(\varphi\in\Lambda^{p,q}\) and \(\psi\in\Lambda^{k,l}\),

\[\begin{split}\begin{aligned} d^\nabla(\varphi\owedge\psi)&=d^\nabla\varphi\owedge\psi+(-1)^p\varphi\owedge d^\nabla\psi,\\ d^{\nabla,\prime}(\varphi\owedge\psi)&=d^{\nabla,\prime}\varphi\owedge\psi+(-1)^q\varphi\owedge d^{\nabla,\prime}\psi. \end{aligned}\end{split}\]

[5]:

def d_leibniz(phi, psi, slot):
    degree = phi.degree_left if slot == "left" else phi.degree_right
    return dg.Wedge(mf.d_cov(phi, slot=slot), psi) + (-1) ** degree * dg.Wedge(
        phi, mf.d_cov(psi, slot=slot)
    )


with TaskManager():
    for phi, psi in ((A11, B11), (A12, B11), (A21, B11)):
        for slot in ("left", "right"):
            left = mf.d_cov(dg.Wedge(phi, psi), slot=slot)
            right = d_leibniz(phi, psi, slot)
            assert l2_error(left, right) < TOL

The trace relates the two coderivatives to the derivative in the opposite slot:

\[\delta^{\nabla,\prime}=-\Tr d^\nabla-d^\nabla\Tr,\qquad \delta^\nabla=-\Tr d^{\nabla,\prime}-d^{\nabla,\prime}\Tr.\]

[6]:

with TaskManager():
    delta_right = mf.delta_cov(A22, slot="right")
    trace_formula_right = -mf.Trace(mf.d_cov(A22, slot="left")) - mf.d_cov(
        mf.Trace(A22), slot="left"
    )
    delta_left = mf.delta_cov(A22, slot="left")
    trace_formula_left = -mf.Trace(mf.d_cov(A22, slot="right")) - mf.d_cov(
        mf.Trace(A22), slot="right"
    )

    assert l2_error(delta_right, trace_formula_right) < TOL
    assert l2_error(delta_left, trace_formula_left) < TOL

The interior product is also a graded derivation. The method ContractSlot preserves the unselected slot and contracts the selected slot with a vector field.

[7]:

def contraction_leibniz(phi, psi, slot):
    degree = phi.degree_left if slot == "left" else phi.degree_right
    return dg.Wedge(mf.ContractSlot(phi, X, slot=slot), psi) + (
        -1
    ) ** degree * dg.Wedge(phi, mf.ContractSlot(psi, X, slot=slot))


with TaskManager():
    for slot in ("left", "right"):
        left = mf.ContractSlot(dg.Wedge(A11, B11), X, slot=slot)
        right = contraction_leibniz(A11, B11, slot)
        assert l2_error(left, right) < TOL

Integration by parts

For \(\chi\in\Lambda^{p,q}(\Omega)\) and \(\psi\in\Lambda^{p+1,q}(\Omega)\), the left-slot operators satisfy

\[\int_\Omega \langle d^\nabla\chi,\psi\rangle\,\omega =\int_\Omega \langle\chi,\delta^\nabla\psi\rangle\,\omega -\int_{\partial\Omega}\langle P_F\chi,P_n\psi\rangle\,\omega_{\partial\Omega}.\]

Here ProjectDoubleForm(..., left="F") projects the left slot tangentially and left="n" contracts it with the unit normal. We chose normal_sign=-1, so the normal is inward and the boundary term has the displayed minus sign. The form inner product includes the factorial normalization in each alternating slot.

[8]:

chi = A11
psi = A21

with TaskManager():
    volume_left = Integrate(
        mf.InnerProduct(mf.d_cov(chi), psi, forms=True)
        * mf.VolumeForm(VOL)
        * ngsolve.dx(bonus_intorder=6),
        mesh,
    )
    volume_right = Integrate(
        mf.InnerProduct(chi, mf.delta_cov(psi), forms=True)
        * mf.VolumeForm(VOL)
        * ngsolve.dx(bonus_intorder=6),
        mesh,
    )
    boundary = Integrate(
        mf.InnerProduct(
            mf.ProjectDoubleForm(chi, left="F"),
            mf.ProjectDoubleForm(psi, left="n"),
            forms=True,
        )
        * mf.VolumeForm(BND)
        * ngsolve.dx(element_boundary=True, bonus_intorder=6),
        mesh,
    )

assert abs(volume_left - volume_right + boundary) < 10 * TOL
print(f"integration-by-parts error: {abs(volume_left - volume_right + boundary):.3e}")

integration-by-parts error: 2.367e-08