[Leopart, Geopart] HDG Stokes velocity interpolation

2

It looks like you’re observing suboptimal convergence of the Stokes system in a non-convex domain. This is a pretty well known result and is typically noticed by singularities in the pressure approximation.

Consider my modification of your example. Let

\Vert \mathbf{v} \Vert_\text{noslip} = \sqrt{\int_{\partial\Omega_\text{wall}} \mathbf{v}^2 \mathrm{d} s}

such that we measure \Vert \mathbf{u}_h \Vert_\text{noslip} and \Vert \bar{\mathbf{u}}_h \Vert_\text{noslip}.

You should see that \Vert \bar{\mathbf{u}}_h \Vert_\text{noslip} = 0 by virtue of the strong imposition of boundary data, however \Vert \mathbf{u}_h \Vert_\text{noslip} \leq C h converges at the suboptimal rate. I’d expect to see this in typical FE schemes that weakly impose boundary data in non-convex computational domains, e.g: Nitsche’s method, DG methods, Langrange multiplier methods etc.

import numpy as np
import geopart.stokes.incompressible
import dolfin
import dolfin_dg as dg

dolfin.parameters["refinement_algorithm"] = "plaza_with_parent_facets"

#Load mesh and boundary
mesh = dolfin.Mesh()
hdf = dolfin.HDF5File(mesh.mpi_comm(), "mesh_boundaries.h5", "r")
hdf.read(mesh, "/mesh", False)

boundaries = dolfin.MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
mark = {"Internal":0, "wall": 1,"inlet": 2,"outlet": 3 }
hdf.read(boundaries, "/boundaries")
hdf.close()

m_max = 3
u_noslip_l2 = np.zeros(m_max, dtype=np.double)
ubar_noslip_l2 = np.zeros(m_max, dtype=np.double)
hmaxs = np.zeros(m_max, dtype=np.double)
for m in range(m_max):
    for j in range(m):
        cf = dolfin.MeshFunction("bool", mesh, mesh.topology().dim(), True)
        mesh2 = dolfin.refine(mesh, cf, False)
        boundaries = dolfin.adapt(boundaries, mesh2)
        mesh = mesh2

    #Define HDG element and function space
    element_cls = geopart.stokes.incompressible.HDG2()
    W = element_cls.function_space(mesh)

    ds = dolfin.Measure('ds',domain=mesh,subdomain_data=boundaries)
    n = dolfin.FacetNormal(mesh)

    #Define boundary condition
    U = element_cls.solution_variable(W)

    mu = dolfin.Constant(0.001)         # dynamic viscosity
    p_in = dolfin.Constant(0.1)         # pressure inlet
    p_out = dolfin.Constant(0.0)        # pressure outlet
    noslip = dolfin.Constant([0.0]*mesh.geometry().dim())  # no-slip wall

    #Boundary conditions
    gN1 = (- p_out*dolfin.Identity(mesh.geometry().dim())) * n
    Neumann_outlet=dg.DGNeumannBC(ds(mark["outlet"]), gN1)
    gN2 = (- p_in*dolfin.Identity(mesh.geometry().dim())) * n
    Neumann_inlet=dg.DGNeumannBC(ds(mark["inlet"]), gN2)
    Dirichlet_wall=dg.DGDirichletBC(ds(mark["wall"]), noslip)

    weak_bcs = [Dirichlet_wall,Neumann_inlet,Neumann_outlet]

    #Body force term
    f = dolfin.Constant([0.0]*mesh.geometry().dim())
    model = geopart.stokes.incompressible.StokesModel(eta=mu,f=f)

    #Form Stokes
    A, b = dolfin.PETScMatrix(), dolfin.PETScVector()
    element_cls.solve_stokes(W, U, (A, b), weak_bcs, model)
    uh, ph = element_cls.get_velocity(U), element_cls.get_pressure(U)
    uhbar = U[1].sub(0)

    u_noslip_l2[m] = dolfin.assemble(uh ** 2 * ds(mark["wall"])) ** 0.5
    ubar_noslip_l2[m] = dolfin.assemble(uhbar ** 2 * ds(mark["wall"])) ** 0.5
    hmaxs[m] = dolfin.MPI.max(mesh.mpi_comm(), mesh.hmax())

dolfin.XDMFFile("ph.xdmf").write_checkpoint(ph, "p")

print("||u||_noslip", u_noslip_l2)
print("||u||_noslip rates",
      np.log(u_noslip_l2[1:] / u_noslip_l2[:-1])
      / np.log(hmaxs[1:] / hmaxs[:-1]))
print("||ubar||_noslip", ubar_noslip_l2)

The pressure solution p_h on cells. Note the singularities at the re-entrant corners.

image
image751×666 45.5 KB

I can’t find a reference in the context of HDG problems, but perhaps this paper regarding the DG approximation of Stokes in a non convex domain is a good place to start.

That all being said. A key result of HDG methods is the error in the velocity approximation is independent of the pressure approximation. So perhaps there’s something I’m missing about EDG-HDG vs HDG.

3 Likes