[color=#000000]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:[/color]
[color=#000000]# -- coding: utf-8 --[/color]
[color=#000000]from netgen.geom2d import SplineGeometry[/color]
[color=#000000]from netgen.meshing import Mesh[/color]
[color=#000000]from ngsolve import *[/color]
[color=#000000]background_domain = SplineGeometry()[/color]
[color=#000000]background_domain.AddRectangle( (-5/3, -5/3), (5/3, 5/3))[/color]
[color=#000000]mesh = Mesh(background_domain.GenerateMesh(maxh=0.5, quad_dominated=False))[/color]
[color=#000000]nmesh = specialcf.normal(mesh.dim)[/color]
[color=#000000]t = Parameter(0)[/color]
[color=#000000]exact_phi = CoefficientFunction((x-0.2*t)2 + y2 + z**2 - 1)[/color]
[color=#000000]flow = CoefficientFunction((0.2,0))[/color]
[color=#000000]V = L2(mesh, order=1, all_dofs_together=True, dgjumps = True) # DG space[/color]
[color=#000000]phi,psi = V.TnT()[/color]
[color=#000000]startfunc = GridFunction(V)[/color]
[color=#000000]startfunc.Set(exact_phi)[/color]
[color=#000000]res = startfunc.vec.CreateVector()[/color]
[color=#000000]t.Set(1/10)[/color]
[color=#000000]upw_flux = flownmesh * IfPos(flownmesh, phi, phi.Other(bnd=exact_phi))[/color]
[color=#000000]c_nonass = BilinearForm(V, nonassemble=True)[/color]
[color=#000000]c_nonass += (-flow * grad(psi) * phi).Compile() * dx[/color]
[color=#000000]c_nonass += (upw_flux * (psi-psi.Other())).Compile() * dx(skeleton=True)[/color]
[color=#000000]c_nonass += (upw_flux * psi).Compile() * ds(skeleton=True)[/color]
[color=#000000]c_ass = BilinearForm(V, nonassemble=False)[/color]
[color=#000000]c_ass += (-flow * grad(psi) * phi).Compile() * dx[/color]
[color=#000000]c_ass += (upw_flux * (psi-psi.Other())).Compile() * dx(skeleton=True)[/color]
[color=#000000]c_ass += (upw_flux * psi).Compile() * ds(skeleton=True)[/color]
[color=#000000]c_ass.Assemble()[/color]
[color=#000000]res.data = c_ass.mat * startfunc.vec - c_nonass.mat * startfunc.vec[/color]
[color=#000000]print(sum(res))[/color]
[color=#000000]Output: [/color]4.294511888968005
[color=#000000]The output should be zero since the bilinear forms are the same (except the nonassemble flag), shouldn’t it?
Can anybody explain this [/color]behavior [color=#000000]or is this a bug?
Thanks for your help!
Best wishes,
Paul[/color]