[LeoPart,GeoPart] HDG Stokes has large error in anisotropic mesh

I’m using geopart and leopart to solve Stokes flow for a while.

But I just found that the HDG scheme is very sensitive to mesh quality compare to the classic Taylor-Hood scheme. Any ideas?

Here is the testing case of simple square duct flow using several anisotropic meshes:

HDG Stokes velocity field (NZ/NX=5,4,3):

Taylor-Hood Stokes velocity field (NZ/NX=5,4,3):

Package:

import numpy as np
import dolfin
from dolfin import *

from mpi4py import MPI as pyMPI

from leopart import StokesStaticCondensation, FormsStokes
import geopart.stokes.incompressible
import dolfin_dg as dg

comm = pyMPI.COMM_WORLD
mpi_comm = MPI.comm_world


#mark whole boundary, inflow and outflow will overwrite)
class Noslip(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1.0)

#Create a unit box mesh
n_ele = 6
aspect_ratio = 4

mesh = BoxMesh(comm, Point(0.0, 0.0,0.0), Point(1.0, 1.0, 1.0),
                     n_ele, n_ele, n_ele*aspect_ratio)

#read mesh and boundaries from file
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)

mark = {"Internal":0, "wall": 1,"inlet": 2,"outlet": 3 }

boundaries.set_all(mark["Internal"])
wall=Noslip()
wall.mark(boundaries, mark["wall"])
left = Left()
left.mark(boundaries, mark["inlet"])
right = Right()
right.mark(boundaries, mark["outlet"])

#read viscosity coefficient from file
mu = Constant(0.001)

#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.create_solution_variable(W)

p_in = dolfin.Constant(1.0)         # 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.StokesModel(eta=mu,f=f)

#Form and solve 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)

#Output solution p,u to paraview
dolfin.XDMFFile("pressure.xdmf").write_checkpoint(ph, "p")
dolfin.XDMFFile("velocity.xdmf").write_checkpoint(uh, "u")

flux = [dolfin.assemble(dolfin.dot(uh, n)*ds(i)) for i in range(len(mark))]

if comm.Get_rank() == 0:
    for key, value in mark.items():
        print("Flux_%s= %.15lf"%(key,flux[value]))

You can find example code and mesh file here:
https://github.com/BinWang0213/TemporaryProject/tree/master/hdg_test/anistropic

After more testings, I found that Laplacian formulation with non-physical traction Neumann BC will fix this issue.

So I guess the implementation of weak form of HDG in LeoPart has some problems in it. Do you have any ideas?

The Laplacian formulation of the Navier-Stokes equations,

-\mu \nabla ^2\mathbf{u}+\nabla p=f

The above strong formulation is equivalent to the following divergence formulation arising from continuum mechanics (used in LeoPart):

-\nabla \cdot \left( \mathbf{\sigma } \right) =\mathbf{f} \\ \mathbf{\sigma }=-p\mathbf{I}+2\mu \mathbf{\varepsilon } \\ \mathbf{\varepsilon }=\frac{1}{2}\left( \nabla \mathbf{u}+\nabla \mathbf{u}^T \right)

Reference:

  1. Some thoughts on the Navier-Stokes equations — Harish Narayanan
  2. Page 26 https://www.math.colostate.edu/~bangerth/videos/676/slides.21.5.pdf

Here are my corrected ufl_forms() for the HDG Stokes to resolve this issue. I simply replace all sym(grad(w)) by grad(w) and 2*nu by nu.

def ufl_forms(self, nu, f, BCtype='physicalBC'):
        (w, q, wbar, qbar) = self.test_functions()
        (u, p, ubar, pbar) = self.trial_functions()

        # Infer geometric dimension
        zero_vec = np.zeros(self.gdim)

        ds = self.ds
        n = self.n
        he = self.he
        alpha = self.alpha
        beta_stab = self.beta_stab
        facet_integral = self.facet_integral

        pI = p * Identity(self.mixedL.sub(1).ufl_cell().topological_dimension())
        pbI = pbar * Identity(self.mixedL.sub(1).ufl_cell().topological_dimension())

        if(BCtype=='physicalBC'): #Physical traction of Neumann BC

            # Upper left block
            # Contribution comes from local momentum balance
            AB = (
                inner(2 * nu * sym(grad(u)), grad(w)) * dx
                + facet_integral(dot(-2 * nu * sym(grad(u)) * n + (2 * nu * alpha / he) * u, w))
                + facet_integral(dot(-2 * nu * u, sym(grad(w)) * n))
                - inner(pI, grad(w)) * dx
            )
            # Contribution comes from local mass balance
            BtF = -dot(q, div(u)) * dx - facet_integral(beta_stab * he / (nu + 1) * dot(p, q))
            A_S = AB + BtF

            # Upper right block
            # Contribution from local momentum
            CD = (
                facet_integral(-alpha / he * 2 * nu * inner(ubar, w))
                + facet_integral(2 * nu * inner(ubar, sym(grad(w)) * n))
                + facet_integral(dot(pbI * n, w))
            )
            H = facet_integral(beta_stab * he / (nu + 1) * dot(pbar, q))
            G_S = CD + H

            # Transpose block
            CDT = (
                facet_integral(-alpha / he * 2 * nu * inner(wbar, u))
                + facet_integral(2 * nu * inner(wbar, sym(grad(u)) * n))
                + facet_integral(qbar * dot(u, n))
            )
            HT = facet_integral(beta_stab * he / (nu + 1) * dot(p, qbar))
            G_ST = CDT + HT

            # Lower right block, penalty on ds(98) approximates free-slip
            KL = (
                facet_integral(alpha / he * 2 * nu * dot(ubar, wbar))
                - facet_integral(dot(pbar * n, wbar))
                + Constant(1e12) / he * inner(outer(ubar, wbar), outer(n, n)) * ds(98)
            )
            LtP = -facet_integral(dot(ubar, n) * qbar) - facet_integral(
                beta_stab * he / (nu + 1) * pbar * qbar
            )
            B_S = KL + LtP
        
        else:#Normal stress neumann BC

            # Upper left block
            # Contribution comes from local momentum balance
            AB = (
                inner(nu * grad(u), grad(w)) * dx
                + facet_integral(dot(-nu * grad(u) * n + (nu * alpha / he) * u, w))
                + facet_integral(dot(-nu * u, grad(w) * n))
                - inner(pI, grad(w)) * dx
            )
            # Contribution comes from local mass balance
            BtF = -dot(q, div(u)) * dx - facet_integral(beta_stab * he / (nu + 1) * dot(p, q))
            A_S = AB + BtF

            # Upper right block
            # Contribution from local momentum
            CD = (
                facet_integral(-alpha / he * nu * inner(ubar, w))
                + facet_integral(nu * inner(ubar, grad(w) * n))
                + facet_integral(dot(pbI * n, w))
            )
            H = facet_integral(beta_stab * he / (nu + 1) * dot(pbar, q))
            G_S = CD + H

            # Transpose block
            CDT = (
                facet_integral(-alpha / he * nu * inner(wbar, u))
                + facet_integral(nu * inner(wbar, grad(u) * n))
                + facet_integral(qbar * dot(u, n))
            )
            HT = facet_integral(beta_stab * he / (nu + 1) * dot(p, qbar))
            G_ST = CDT + HT

            # Lower right block, penalty on ds(98) approximates free-slip
            KL = (
                facet_integral(alpha / he * nu * dot(ubar, wbar))
                - facet_integral(dot(pbar * n, wbar))
                + Constant(1e12) / he * inner(outer(ubar, wbar), outer(n, n)) * ds(98)
            )
            LtP = -facet_integral(dot(ubar, n) * qbar) - facet_integral(
                beta_stab * he / (nu + 1) * pbar * qbar
            )
            B_S = KL + LtP

        # Righthandside
        Q_S = dot(f, w) * dx
        S_S = facet_integral(dot(Constant(zero_vec), wbar))
        return {"A_S": A_S, "G_S": G_S, "G_ST": G_ST, "B_S": B_S, "Q_S": Q_S, "S_S": S_S}

This correction posted didn’t really fix the problem. We can observe large velocity fluctuation when solving a problem with a poor mesh quality (aspect ratio > 2).

The problem is resolved. The HDG for 3D with high aspect ratio cells requires a large interior penalty parameter (alpha). For the same problem above,

In the extreme case above, when all of the mesh cells have a large aspect ratio, a small interior penalty parameter will result in an unstable solution shown in the first post.

As the stabilizer is calculated as \frac{\alpha}{h_e}, I wonder if the way of calculating mesh size h_e for a large aspect ratio cell affects the stability here? @nate @jpdean

1 Like

As Bin and I have discussed, we wonder if @jpdean has experienced HDG-IP methods being more sensitive to the penalty parameter.

What does worry me here is that the IP parameter should be independent of the mesh, including the aspect ratio. Am I right here?

EDIT: Don’t forget that the penalisation term typically chosen for interior penalty methods incorporates the polynomial degree \mathsf{k} whereby \alpha \frac{\max(1, \mathsf{k}^2)}{h_e}.

Hi @Wang_Bin and @nate,

Sorry for my slow reply, I have been away for a week.

As far as I’m aware, the minimum value of the penalty parameter required to ensure coercivity typically depends on the constant appearing in the discrete trace inequality, which in turn depends on mesh regularity parameters, and therefore on the aspect ratio of the cells.

This isn’t ideal, as I realise that high aspect ratio cells are commonly encountered in meshes of real engineering geometry. However, there are other reasons to avoid poor quality meshes. For instance, it has been shown that when working with properly configured preconditioned iterative solvers, it is more efficient to use high quality meshes at the expense of increased cell count; see for instance https://doi.org/10.1016/j.finel.2018.11.002

Thanks,

Joe

3 Likes

Excellent insight, thanks @jpdean.

1 Like