4. Convection equation

from ngsolve import *
from ngsolve.webgui import Draw
from time import sleep, time

mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))
# mesh = Mesh(unit_cube.GenerateMesh(maxh=0.05))

wind = CF( (y-0.5, 0.5-x))

order = 5
fesT = L2(mesh, order=order)
fesF = FacetFESpace(mesh, order=order)
fes = fesT*fesF
eT,eF = fes.embeddings
u,uhat = fes.TrialFunction()
v = fesT.TestFunction()
gfu = GridFunction(fesT)
fesT.ndof, fesF.ndof

(19320, 8520)

traceop = fesT.TraceOperator(fesF, average=True)
projinner = Projector(fesF.GetDofs(mesh.Boundaries(".*")), False)
makeboth = (eF@projinner@traceop+eT)

n = specialcf.normal(mesh.dim)
dS = dx(element_vb=BND)
wn = wind*n

convequ = -wind*u*grad(v)*dx + IfPos(wn, wn*u, wn*(2*uhat-u))*v*dS
bf = BilinearForm(convequ, nonlinear_matrix_free_bdb=True, nonassemble=False).Assemble()

massinv = fesT.Mass(1).Inverse() 

convop = massinv@bf.mat@makeboth

gfu.Set (exp(-100*((x-0.75)**2+(y-0.5)**2)))

scene = Draw(gfu, order=order, euler_angles=(-50,0,-10), deformation=True, scale=0.3)

hv = gfu.vec.CreateVector()
tau = 0.002
t = 0
tend = 10
i = 0
ts = time()
with TaskManager(): # pajetrace=10**8):
    while t < tend:
        t += tau
        hv.data = gfu.vec - tau/2 * convop * gfu.vec
        gfu.vec.data -= tau * convop * hv
        i += 1
        if i%20 == 0:
            scene.Redraw()
print ("time =", time()-ts)

time = 4.7213051319122314

# print (bf1.mat.GetOperatorInfo())

# for t in Timers():
#    if "SoA" in t["name"]:
#        print (t)