Hi Joachim
Thanks for your help! I was able to fix it eventually. I am putting the code here so other people can also use it.
from ngsolve import *
from netgen.geom2d import SplineGeometry
from ngsolve.webgui import Draw
from netgen.occ import *
from netgen.webgui import Draw as DrawGeo
shape = Rectangle(2,2).Face().Move((-1,-1,0))
shape.edges.Min(X).name="left"
shape.edges.Max(X).name="right"
shape.edges.Min(Y).name="bottom"
shape.edges.Max(Y).name="top"
time = 0.0
dt = 0.01
mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.25))
fes_main = H1(mesh, order=3, dirichlet="bottom|right|left|top")
gfu_main = GridFunction(fes_main)
new_mesh = Mesh(mesh.ngmesh.Copy())
fes = H1(new_mesh, order=3, dirichlet="bottom|right|left|top", autoupdate=True)
u,v = fes.TnT()
gfu_new = GridFunction(fes, autoupdate=True)
gfu_old = GridFunction(fes, autoupdate=True)
gfu_old.Set((1-y*y)*x)
def Solve(gfu_old, fes, u, v):
b = CoefficientFunction((2*y*(1-x*x),-2*x*(1-y*y)))
gaussp = exp(-6*((x+0.5)*(x+0.5)+y*y))-exp(-6*((x-0.5)*(x-0.5)+y*y))
a = BilinearForm(fes, symmetric=False)
a += 0.01*grad(u)*grad(v)*dx + b*grad(u)*v*dx
a.Assemble()
m = BilinearForm(fes, symmetric=False)
m += u*v*dx
m.Assemble()
mstar = m.mat.CreateMatrix()
mstar.AsVector().data = m.mat.AsVector() + dt * a.mat.AsVector()
invmstar = mstar.Inverse(freedofs=fes.FreeDofs())
f = LinearForm(fes)
f += gaussp*v*dx
f.Assemble()
res = dt * f.vec - dt * a.mat * gfu_old.vec
gfu_new.vec.data = gfu_old.vec + invmstar * res
return gfu_new
def CalcError():
flux = 0.01*grad(gfu_new)
# interpolate finite element flux into H(div) space:
gf_flux.Set(flux)
# Gradient-recovery error estimator
err = 1/0.01*(flux-gf_flux)*(flux-gf_flux)
elerr = Integrate (err, new_mesh, VOL, element_wise=True)
maxerr = max(elerr)
#print ("maxerr = ", maxerr)
for el in new_mesh.Elements():
new_mesh.SetRefinementFlag(el, elerr[el.nr] > 0.25*maxerr)
cnt = 1;time = cnt*dt; tend = 2; nsamples = 10; sample_int = int(floor(tend / dt / nsamples)+1)
gfut = GridFunction(gfu_main.space,multidim=0)
gfu_main.Set(gfu_old)
gfut.AddMultiDimComponent(gfu_main.vec)
while time < tend - 0.5 * dt:
space_flux = HDiv(new_mesh, order=2, autoupdate=True)
gf_flux = GridFunction(space_flux, "flux", autoupdate=True)
with TaskManager():
while fes.ndof < 100000:
Solve(gfu_old, fes, u, v)
CalcError()
new_mesh.Refine()
gfu_main.Set(Solve(gfu_old, fes, u, v))
print("\r",time,end="")
new_mesh = Mesh(mesh.ngmesh.Copy())
fes = H1(new_mesh, order=3, dirichlet="bottom|right|left|top", autoupdate=True)
gfu_old = GridFunction(fes, autoupdate=True)
gfu_new = GridFunction(fes, autoupdate=True)
gfu_old.Set(gfu_main)
if cnt % sample_int == 0:
gfut.AddMultiDimComponent(gfu_old.vec)
u,v = fes.TnT()
cnt += 1; time = cnt * dt
Draw(gfut, mesh, interpolate_multidim=True, animate=True)