[Leopart, Geopart] HDG Stokes velocity interpolation

1

Domain
Domain938×640 249 KB

Green: Neumann BC Value=0.1
Red: Neumann BC Value=0.0
Cyan: no-slip Dirichelt BC

I’m using geopart and leopart to learn the HDG for Stokes problem. In their paper, one of important properties is that HDG will give us the pointwise divergence-free velocity solution.

But when I use fenics built-in coefficient interpolation uh([x,y]), it returns non-zero velocity at the no-slip wall.

In the Leopart, the cell-wise velocity is interpolated by

#https://bitbucket.org/jakob_maljaars/leopart/src/master/source/cpp/advect_particles.cpp#lines-385
Utils::return_expansion_coeffs(coeffs, *ci, &uh_step);
Utils::return_basis_matrix(basis_mat.data(), _P->x(ci->index(), i), *ci, _element);
Eigen::VectorXd u_p = basis_mat * exp_coeffs;

How can I correctly or manully interpolate the pointwise divergence-free velocity in each cell?

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

code:

import numpy as np
import geopart.stokes
import dolfin
import dolfin_dg as dg
import ufl
from ufl import sym, grad, Identity, div, dx
import matplotlib.pyplot as plt

#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()

#Define HDG element and function space
element_cls = geopart.stokes.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.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)


#visulize solution in Paraview
ufile_pvd = dolfin.File("velocity.pvd")
ufile_pvd << uh
pfile_pvd = dolfin.File("pressure.pvd")
pfile_pvd << ph #pressure is flipped for symmetric matrice

bd_pvd = dolfin.File("boundaries_mesh.pvd")
bd_pvd << boundaries


#Check boundary wall velocity

#Plot velocity u along bottom edge (0,0)->(12,0)
x=np.linspace(0,12,10)
vel_interp=np.array([uh([i,0.0]) for i in x])
print('BoundaryVelocity @ Line (0,0)->(12,0)=',np.sum(np.abs(vel_interp)))
# Non-zero velocity at non-slip wall 
#BoundaryVelocity @ Line (0,0)->(12,0)= 0.0534387418074

plt.plot(x,vel_interp[:,0],label='Velocity x')
plt.plot(x,vel_interp[:,1],label='Velocity y')
plt.legend()
plt.show()
1 Like
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
3

Hi Nate.

Thanks for your detailed analysis. I don’t have a strong background in math. It looks like we have to refine the mesh to avoid the “fluid leaking” issue in HDG.

As you can see in my original coarse mesh (figure below), the fuid is even “leaking” (velocity vector point out of our domain) on the bottom flat wall.

image
image952×690 199 KB

From the engineering perspective, it will be great if we have divergence-free velocity field (like finite volume, DG, etc) and strongly imposed boundary condition (like in Taylor-Hood scheme).

  1. Is it possible to strongly prescribe the no-slip boundary condition by using geopart/leopart or HDG, DG etc?

  2. Also, how can I manully interpolate the pointwise divergence-free velocity in each cell in python environment?

Thanks a lot for your comments!
Bin

4

Hi there!

It looks like the problem lies in your application of the no-slip boundary condition. In the lines

Dirichlet_wall=dg.DGDirichletBC(ds(mark["wall"]), noslip)
weak_bcs = [Dirichlet_wall,Neumann_inlet,Neumann_outlet]

you are weakly imposing the no-slip condition on the walls. Instead, it should be enforced directly in the function space for the facet velocity, \bar{u}, by simply setting the degrees of freedom on the walls to zero.

Hope this helps,

Joe

2 Likes
5

GeoPart is written poorly in the syntax for imposing BCs. It was initially designed to quickly prototype several Stokes FE pairs as supporting information for a paper. BCs are contained in DGDirichletBC/DGNeumnnBCs and then extracted for strong imposition or weak imposition depending on whether the underlying FEM is conforming or nonconforming, respectively.

So in the case of HDG, GeoPart does indeed strongly impose the BCs on the space \bar{V} which is reflected in the noslip measure \Vert\bar{\mathbf{u}}_h\Vert_\text{noslip} = 0. The weak coupling of this facet data to the cell data appears to be causing the issue though, as I noted that \Vert \mathbf{u}_h \Vert_\text{noslip} converges at suboptimal rates \backsim h.

If I get time I’ll design a MWE.

1 Like
6

Just in addition to the answers of @jpdean and @nate:

  • to get a better picture of what’s precisely happening at the wall, plot the discontinuous velocity field \mathbf{u}_h rather than an interpolated version - which in fact is done when writing output to pvd format. So write the function to an xdmf file. See Nate’s snippet for an example.
  • you can evaluate the pointwise div-free field \mathbf{u}_h in various ways. The easiest - but probably also the slowest - is just evaluating your field uh at a Point. If I remember well, you just can do something like uh(Point(x,y,z)). In case you have an array of coordinates on which you want to evaluate uh, you could use the particles.interpolate function offered by LEoPart as it will be much faster. As an example, see https://bitbucket.org/jakob_maljaars/leopart/src/70124646f539397c470adc39417009f35fc4eae7/unit_tests/test_2d_project.py?at=master#lines-167.
2 Likes
8

Hi Jakob,

Thanks for your detailed instruction. As you can see in the figure below, the velocity field looks great!

But we still have an error of 1e-2 on the non-slip wall. And the error makes the noslip wall behaves like a slip wall.

  • Can we further “fix” this error somehow? This is non-physical solution in term of engineering usage even it is divergence-free.

  • BTW, which method was actually implemented in the Geopart/Leopart now? EDG-HDG (Rhebergen & Wells, 2020) or HDG (Rhebergen & Wells, 2018) ? I remember that in the Leopart HDG was implemented.

[1] Rhebergen, S., & Wells, G. N. (2018). A hybridizable discontinuous Galerkin method for
the Navier–Stokes equations with pointwise divergence-free velocity field. Journal of
Scientific Computing, 76 (3), 1484–1501.

[2] Rhebergen, S., & Wells, G. N. (2020). An embedded–hybridized discontinuous Galerkin
finite element method for the Stokes equations. Computer Methods in Applied
Mechanics and Engineering, 358 , 112619.

image
image908×701 232 KB

image

Here is the updated code:

Post-processing (Velocity interpolation)

#Assemble and solve the problem is the same as in the above code

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

#Divergence-free velocity interpolation
#We plot the velocity profile on the bottom edge
# (0, 0)->(12, 0)
NumPts=50
x,y=np.linspace(0,12,NumPts),np.ones(NumPts)*0.0
probe_pts=np.vstack([x,y]).T

#Method1 (slow)
u_interp=np.array([uh( dolfin.Point(i,j) ) for i,j in probe_pts])

#Plot velocity profile
import matplotlib.pyplot as plt
plt.title("uh along bottom edge")
plt.plot(probe_pts[:,0],u_interp[:,0],'ro',label='Velocity x')
plt.plot(probe_pts[:,0],u_interp[:,1],'bo',label='Velocity y')
plt.legend()
plt.show()

#Method2 (Fast)
u_interp=np.zeros([NumPts,2])
probe_particles=particles(probe_pts, [u_interp], mesh)
probe_particles.interpolate(uh,1)
u_interp=probe_particles.return_property(mesh,1)
probe_pts=probe_particles.positions() #noted that pts is reordered in Leopart
9

As Nate mentioned, some slip is to be expected due to the weak coupling between u and \bar{u}. However, since convergence is suboptimal and (assuming the solution is being interpolated correctly) the slip you are seeing on the boundary is quite large compared to the solution, it is possible something still might be wrong.

What value are you using for the penalty parameter? If it is too small, it might explain what we are seeing here.

10

Hi Joe,

I think they used \alpha=6k^2, in my testing code \alpha=6*2^2=24.
See code in Geopart here:

Does this value is too small?

11

I see. That value looks fine to me. In that case, I’m not sure and will have to have a think!

Thanks,

Joe