Problem with solving Stokes with RT1-DG0?

General
,
1

Hi all,

I am testing solving a basic Stokes problem with some different elements, and I have problems getting the correct solution with RT1-DG0. It works fine for RT2-DG0 and P2-P1 Taylor-Hood elements. Is there something I’m missing to do for the RT1-DG0 pair?

Here’s the MWE based on the Stokes solve in demo_navier_stokes.py:

import numpy as np
from numpy import sin, cos, pi
from mpi4py import MPI
from petsc4py import PETSc
import ufl
from ufl import inner, outer, grad, div, dx, dS, ds, avg, jump
import dolfinx
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block


N = 32 * 4
mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, N, N, dolfinx.mesh.CellType.triangle
)

# # P2-P1 works fine
# V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim,)))
# Q = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))

# RT2-DG0 works fine
V = dolfinx.fem.functionspace(mesh, ("RT", 2))
Q = dolfinx.fem.functionspace(mesh, ("DG", 0))

# # RT1-DG0 doesn't work
# V = dolfinx.fem.functionspace(mesh, ("RT", 1))
# Q = dolfinx.fem.functionspace(mesh, ("DG", 0))

u, p = ufl.TrialFunction(V), ufl.TrialFunction(Q)
v, q = ufl.TestFunction(V), ufl.TestFunction(Q)
n = ufl.FacetNormal(mesh)
h = ufl.CellDiameter(mesh)


def tensor_jump(v, n):
    return outer(v("+"), n("+")) + outer(v("-"), n("-"))


beta = dolfinx.fem.Constant(mesh, 10.0)
a_interior = (
    inner(grad(u), grad(v)) * dx
    - inner(avg(grad(u)), tensor_jump(v, n)) * dS
    - inner(tensor_jump(u, n), avg(grad(v))) * dS
    + beta / avg(h) * inner(jump(u), jump(v)) * dS
)
a_bdry = (
    -inner(grad(u), outer(v, n)) * ds
    - inner(outer(u, n), grad(v)) * ds
    + beta / h * inner(u, v) * ds
)
a_h = a_interior + a_bdry

b_h = -inner(div(u), q) * dx
bT_h = -inner(p, div(v)) * dx

a = dolfinx.fem.form([[a_h, bT_h], [b_h, None]])


def f_expr(x):
    # fmt: off
    return np.stack(
        [
            2*pi*sin(2*pi*x[1])*(cos(2*pi*x[0]) - 2*pi*pi*cos(2*pi*x[0]) + pi*pi),
            2*pi*sin(2*pi*x[0])*(cos(2*pi*x[1]) + 2*pi*pi*cos(2*pi*x[1]) - pi*pi)
        ]
    )
    # fmt: on


f = dolfinx.fem.Function(V)
f.interpolate(f_expr)
L0 = inner(f, v) * dx
L1 = inner(dolfinx.fem.Constant(mesh, 0.0), q) * dx
L = dolfinx.fem.form([L0, L1])


def zero(x):
    return np.stack((np.zeros_like(x[0]), np.zeros_like(x[1])))


bdry_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
bdry_dofs = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim - 1, bdry_facets)

u_D = dolfinx.fem.Function(V)
u_D.interpolate(zero)
bc_u = dolfinx.fem.dirichletbc(u_D, bdry_dofs)
bcs = [bc_u]

A = dolfinx.fem.petsc.assemble_matrix_block(a, bcs=bcs)
A.assemble()
b = dolfinx.fem.petsc.assemble_vector_block(L, a, bcs=bcs)

ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
opts = PETSc.Options()
opts["mat_mumps_icntl_14"] = 80
# opts["mat_mumps_icntl_24"] = 1
# opts["mat_mumps_icntl_25"] = 0
ksp.setFromOptions()
x = A.createVecRight()
ksp.solve(b, x)

u_h, p_h = dolfinx.fem.Function(V), dolfinx.fem.Function(Q)
offset = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
u_h.x.array[:offset] = x.array[:offset]


# From demo_navier_stokes:
def domain_average(mesh, v):
    """Compute the average of a function over the domain"""
    vol = mesh.comm.allreduce(
        dolfinx.fem.assemble_scalar(
            dolfinx.fem.form(dolfinx.fem.Constant(mesh, 1.0) * dx)
        ),
        op=MPI.SUM,
    )
    return (1 / vol) * mesh.comm.allreduce(
        dolfinx.fem.assemble_scalar(dolfinx.fem.form(v * dx)), op=MPI.SUM
    )


p_h.x.array[: (len(x.array) - offset)] = x.array[offset:]
p_h.x.array[:] -= domain_average(mesh, p_h)

W = dolfinx.fem.functionspace(mesh, ("DG", 2, (mesh.geometry.dim,)))
u_viz = dolfinx.fem.Function(W)
u_viz.interpolate(u_h)

t = 0.0
u_file = dolfinx.io.VTXWriter(mesh.comm, "u.bp", u_viz)
p_file = dolfinx.io.VTXWriter(mesh.comm, "p.bp", p_h)
u_file.write(t)
p_file.write(t)