Element-wise post-processing in HDG method

General
1

Hi everyone. I am new to NGSolve. I am trying to do element-wise post-processing in the HDG method for the Poisson equation. Can anyone suggest how to do that, or provide any examples? Thank you!

2

Look at examples around here:

Joachim

1 Like
3

Thank you, Joachim, those are really helpful. Now that I have successfully applied HDG to the 2D Poisson equation, I am developing the HDG formulation for the Stokes equation. Here I have attached the formulation I want to do. For some reason, the code below is not converging. I am giving the bilinear formulation because the rest of the code is a minor modification of the 2D Poisson solver, and it already worked. (In the formulation, in some places InnerProduct is used to check whether it has any effect on computation or not)

image
image790×387 38.2 KB

def setup_HDG_system(fes, dt, **opts):

    (q1,q2,u1,u2,u1hat,u2hat,p),(r1,r2,v1,v2,v1hat,v2hat,q) = fes.TnT()
    eps=1e-6
    nu=1.0
    tau = 2.0
    dS = dx(element_boundary=True)
    n = specialcf.normal(2)

    a = BilinearForm(fes, **opts)

    # q + grad(u) = 0
    a += q1*r1*dx
    a += -u1*div(r1)*dx
    a += u1hat*(r1*n)*dS
    a += q2*r2*dx
    a += -u2*div(r2)*dx
    a += u2hat*(r2*n)*dS

    # backward Euler mass term
    a += (1/dt)*u1*v1*dx
    a += (1/dt)*u2*v2*dx

    # momentum
    a += -nu*q1*grad(v1)*dx
    a += -nu*q2*grad(v2)*dx
    a += (nu*InnerProduct(q1,n) + p*n[0] + nu*tau*(u1-u1hat))*v1*dS
    a += (nu*InnerProduct(q2,n) + p*n[1] + nu*tau*(u2-u2hat))*v2*dS

    #Pressure term
    a +=-p*v1.Diff(x)*dx 
    a +=-p*v2.Diff(y)*dx
    a += eps*p*q*dx
    
    #continuity equation
    a += -(u1*q.Diff(x) + u2*q.Diff(y))*dx
    a += (u1hat*n[0] + u2hat*n[1]) * q * dS
    # conservativity
    a += (
        nu*InnerProduct(q1,n)
        + p*n[0]
        + nu*tau*(u1-u1hat)
    )*v1hat*dS

    a += (
        nu*InnerProduct(q2,n)
        + p*n[1]
        + nu*tau*(u2-u2hat)
    )*v2hat*dS
    
    a.Assemble()

    return a
V = VectorL2(mesh, order=order)
W = L2(mesh, order=order)
M = FacetFESpace(
        mesh,
        order=order,
        dirichlet=".*"
    )
Q = L2(mesh, order=order)
fes = V * V * W * W * M * M * Q

Can you have a look at it and let me know if I am doing something wrong?