[Leopart, Geopart] HDG Stokes velocity interpolation

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