ngsolve HDG solver

hdg.py

import numpy as np

from ngsolve import *
from netgen.geom2d import SplineGeometry
from netgen.read_gmsh import ReadGmsh

from ngsolve.webgui import Draw

#Gmsh read support
#gmesh = ReadGmsh('gmsh_msh.msh')
#mesh = Mesh( gmesh ); 


geo = SplineGeometry()
geo.AddRectangle((0, 0), (2, 0.41), bcs = ("default", "outlet", "default", "inlet"))
geo.AddCircle((0.2, 0.2), r=0.05, leftdomain = 0, rightdomain = 1, bc = "default")
mesh = Mesh( geo.GenerateMesh(maxh = 0.08) ); 
#order = 3
#mesh.Curve(order)

#Draw(mesh)
print("ElementBoundary=",mesh.GetBoundaries())
print("NumEles=",mesh.ne)
print("NumEdges=",mesh.nedge)
print("NumVertex=",mesh.nv)

#PkPk-1 scheme

#scheme = 'PkP0'
scheme = 'PkPk-1'


order=2
V1 = HDiv ( mesh, order = order, dirichlet = "default")
V2 = FESpace ( "vectorfacet", mesh, order = order, dirichlet = "default" )
Q = L2(mesh, order = order-1, lowest_order_wb = True) 
V = FESpace ([V1,V2,Q])

#PkP0 scheme with the minimum DOF
if(scheme == 'PkP0'):
    V1 = HDiv ( mesh, order = order, dirichlet = "default", hodivfree = True )
    V2 = FESpace ( "vectorfacet", mesh, order = order, dirichlet = "default" )
    Q = L2(mesh, order = 0, lowest_order_wb = True) 
    V = FESpace ([V1,V2,Q])

#lowest_order_wb = following code
#hdiv dofs + TangentialFacet dofs
# 2d - version
#offset = (order +1) * (order + 2) + 3 * (order +1)
#for el in mesh.Elements():
#    dofs = V.GetDofNrs(el)
#    V.SetCouplingType(dofs[offset] , COUPLING_TYPE.INTERFACE_DOF )


#Change pressure as interface dof, we can not eliminate pressure out of the linear system
#only the local dofs from HDiv will be eliminated
#V.SetCouplingType(V.Range(2), COUPLING_TYPE.INTERFACE_DOF)

V.FinalizeUpdate()

u, uhat, p = V.TrialFunction()
v, vhat, q = V.TestFunction()

nu = 0.001
alpha = 4
condense=True

gradu = CoefficientFunction ( (grad(u),), dims=(mesh.dim,mesh.dim) )
gradv = CoefficientFunction ( (grad(v),), dims=(mesh.dim,mesh.dim) )

n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size

def tang(vec):
    return vec - (vec*n)*n

a = BilinearForm ( V, symmetric=True, condense=condense)
a += SymbolicBFI ( nu*InnerProduct ( gradu, gradv) )
a += SymbolicBFI ( nu*InnerProduct ( gradu.trans * n,  tang(vhat-v) ), element_boundary=True )
a += SymbolicBFI ( nu*InnerProduct ( gradv.trans * n,  tang(uhat-u) ), element_boundary=True )
a += SymbolicBFI ( nu*alpha*order*order/h * InnerProduct ( tang(vhat-v),  tang(uhat-u) ), element_boundary=True )
a += SymbolicBFI ( -div(u)*q -div(v)*p )
a.Assemble()

f = LinearForm ( V )

#Pressure Neumann BCs
p_in = CoefficientFunction( 1.0 )
f += - InnerProduct(p_in*n,v.Trace())*ds(definedon=mesh.Boundaries("inlet"))
p_out = CoefficientFunction( 0.2 )
f += - InnerProduct(p_out*n,v.Trace())*ds(definedon=mesh.Boundaries("outlet"))

f.Assemble()

#Check static condensation
ndof_total = V.ndof
ndof_local=0
for i in range(V.ndof):
    if(V.CouplingType(i)==COUPLING_TYPE.LOCAL_DOF):
        ndof_local+=1
print("V1.ndof =", V1.ndof,V.Range(0)) 
print("V2.ndof =", V2.ndof,V.Range(1)) 
print("Q.ndof  =", Q.ndof,V.Range(2))
print(f"Local DOFs: {ndof_local}, Global DOFs: {ndof_total-ndof_local}, "
      f"CompressRatio: {100-ndof_local/ndof_total*100:.2f}%")

u = GridFunction(V)

#we don't need to explict setup no slip BC here
#noslip = CoefficientFunction( (0,0) )
#gfu.components[0].Set(noslip, definedon=mesh.Boundaries("default"))

r = f.vec.CreateVector()
#Dirichlet BC setup
r.data = f.vec - a.mat * u.vec 
#Condense RHS
if condense:
    r.data += a.harmonic_extension_trans * r
#Solve a condensed system
u.vec.data += a.mat.Inverse(V.FreeDofs(condense)) * r
#Recover the full solution
if condense:
    u.vec.data += a.harmonic_extension * u.vec
    u.vec.data += a.inner_solve * r

vel = u.components[0]
pres = u.components[2]
    
#gives a list of the integral on each part of the bnd of u*n
flux_bd = Integrate(vel * n, mesh, BND,  region_wise = True) 
Markers = list(set(mesh.GetBoundaries()))
flux = {m:0.0 for m in Markers}
for i,bd in enumerate(mesh.GetBoundaries()):
    flux[bd] += flux_bd[i]
print(flux)
    
Draw( Norm(vel), mesh, "absvel(hdiv)")
Draw( Norm(pres), mesh, "pressure")

postProcess.py

#Cell-wise interpolation (continue from hdg.py)
vel = u.components[0]

ele=ElementId(VOL, 118)

#Get element pts
Pts=[]
for Vi in mesh[ele].vertices:
    Pts+=[np.array(mesh[Vi].point)]
Pts=np.array(Pts)

#Get the centroid of cell 118
cntr= np.mean(Pts,axis=0) #physical coordinate
cntr_ngs = mesh(*cntr)    
cntr_local = cntr_ngs.pnt #local coordinate

#----------NGS Velocity interpolation ------
vel = u.components[0]
print('Ngs Interp=', vel(cntr_ngs) ) 

#----------Manully Velocity interpolation ------

#Get shape function of BDMk
vel = u.components[0]
FE=vel.space.GetFE(ele)
shapeFunc=FE.CalcShape(cntr_local[0],cntr_local[1]).NumPy()

#Get solution coeff for BDMk
uIdx=np.array(vel.space.GetDofNrs(ele))
print(uIdx)
coeff = u.vec.FV().NumPy()[uIdx]

#Get piola trans (not sure how to get this in ngsolve)
Jac=np.array( [[ Pts[1,0]-Pts[0,0],Pts[2,0]-Pts[0,0] ],    #[(x2-x1, x3-x1),
               [ Pts[1,1]-Pts[0,1],Pts[2,1]-Pts[0,1]] ] )  # (y2-y1, y3-y1)]
trans=Jac/np.linalg.det(Jac)
print(Jac)

#Interpolate velocity by shapefunc*coeff*trans_mat
shapeFunc_trans = np.dot(trans,shapeFunc.T)
vel_interp = np.dot(shapeFunc_trans, coeff)
print('Manual Interp=',vel_interp)

#Check if two interpolation has the same value
assert np.all( np.isclose(vel_interp, vel(cntr_ngs)) )

TODO list:

• Read mesh from gmsh with phsyscial names

     gmesh = ReadGmsh('gmsh_msh.msh')
     mesh = Mesh( gmesh ); 
     
     print("ElementBoundary=",mesh.GetBoundaries())
     print("NumEles=",mesh.ne)
     print("NumEdges=",mesh.nedge)
     print("NumVertex=",mesh.nv)
     #Draw(mesh)

• Static condensation (hdg.py)
  The codes runs smoothly, but the results are not correct.
• Compute flux on boundary
  In fenics we can do flux_bd = Integrate(vel * n, mesh, BND, region_wise = True)
• Spatially varing fluid viscosity
  Not sure how to set a piecewise material coefficient. Namely, we need to set material properties for each cell.
  In fenics, we can use DG0 and set value for each cell below.

      nu = Function(FunctionSpace(mesh, "DG", 0)) 
      nu.vector()[0] = 0.001 #cell 0 has viscosity 0.001
      nu.vector()[1] = 0.002 #cell 1 has viscosity 0.002
      ...

• Manully cell-wise velocity interpolation (postProcess.py)
  I can not get the same velocity from my manually velocity interpoaltion
  Ngsolve seems to use different trans matrix other than piola transform
• MPI implementation with direct solver of MUMPs
  Installing ngsolve on cluster now.. Will try this later on
• MPI implementation with a preconditioned iterative solver