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
