Problem with nonassemble = True in DG scheme

More
1 year 10 months ago - 1 year 10 months ago #4602 by Schwering
Hello everybody,

I am currently working on a DG scheme to solve a transport equation for a level set function. Thereby, I found the following issue: Defining a bilinear form with the flag nonassemble = True and calculating a matrix-vector product leads to a different result than defining the bilinear form without the flag, assembling it, and performing the matrix-vector product.

Adding the code as a file is not working, so I paste the code here:


# -*- coding: utf-8 -*-
from netgen.geom2d import SplineGeometry
from netgen.meshing import Mesh
from ngsolve import *
 
background_domain = SplineGeometry()
background_domain.AddRectangle( (-5/3, -5/3), (5/3, 5/3))
mesh = Mesh(background_domain.GenerateMesh(maxh=0.5, quad_dominated=False))
nmesh = specialcf.normal(mesh.dim)
t = Parameter(0)
exact_phi = CoefficientFunction((x-0.2*t)**2 + y**2 + z**2 - 1)
flow = CoefficientFunction((0.2,0)) 

V = L2(mesh, order=1, all_dofs_together=True, dgjumps = True)  # DG space
phi,psi = V.TnT() 
startfunc = GridFunction(V)
startfunc.Set(exact_phi) 
res = startfunc.vec.CreateVector() 

t.Set(1/10) 
upw_flux = flow*nmesh * IfPos(flow*nmesh, phi, phi.Other(bnd=exact_phi)) 

c_nonass = BilinearForm(V, nonassemble=True)
c_nonass += (-flow * grad(psi) * phi).Compile() * dx
c_nonass += (upw_flux * (psi-psi.Other())).Compile() * dx(skeleton=True)
c_nonass += (upw_flux * psi).Compile() * ds(skeleton=True) 

c_ass = BilinearForm(V, nonassemble=False)
c_ass += (-flow * grad(psi) * phi).Compile() * dx
c_ass += (upw_flux * (psi-psi.Other())).Compile() * dx(skeleton=True)
c_ass += (upw_flux * psi).Compile() * ds(skeleton=True) 
c_ass.Assemble() 

res.data = c_ass.mat * startfunc.vec - c_nonass.mat * startfunc.vec 
print(sum(res))

Output: 4.294511888968005

The output should be zero since the bilinear forms are the same (except the nonassemble flag), shouldn't it?

Can anybody explain this
behavior or is this a bug?

Thanks for your help!

Best wishes,
Paul

 
Last edit: 1 year 10 months ago by Schwering.
Time to create page: 0.095 seconds