VTKOutput() not works for stokes problem with Tayler-Hood

I try to solve an analytical solution of stokes problem and visualize it in Paraview by modifying the tutorial @ ngsolve/py_tutorials/mpi/mpi_navierstokes.py

But ngsolve seems stucks on export VTK file and prohibit me to load other python script with the following error:

ERROR: Thread already running
errinfo: Thread already running
    while executing
"NGS_LoadPy  $file"
    (menu invoke)

Please check attached script to reproduce the issue.

Thanks,
Bin

Attachment: Stokes-CG_2019-03-03.py

Dear Bin,

The problem is coming from the lines before. The problem is “V.Range(2)”. Your fespace only has 2 components, V.Range(2) does not make sense. It results in an indefinite loop.

Best, Christoph

P.S.:
Further note that for pressure order larger than 1 these lines are not doing what you might expect anyway. The min/max operation will give you the largest/smallest coefficient of your gridfunction. There are however no direct conclusions to draw from the coefficient and the min/max of the represented pressure function (not even at Lagrangian nodes (not a Lagrange basis!)).

Thanks for your reply. The netgen no longer hanging any more.

  1. Strange VTK output
    But it still can not output correct velocity into vtk file. I tested a simple 6 elements domain. The velocity value shown in netgen is (1.3447e-1,1.905e-01).

But in vtk file it wrote a lot of zero (such as 4.70064e-46). Paraview also report following error:

ERROR: In C:\bbd\ecd3383f\build\superbuild\paraview\src\VTK\IO\Legacy\vtkDataReader.cxx, line 1949
vtkUnstructuredGridReader (0000027E6DA553F0): Unsupported data type: 

Here is the python code I used for CG

from math import pi
from ngsolve import *
from netgen.geom2d import SplineGeometry
import numpy as np


u_exact = (sin(x)*sin(x)*cos(y)*sin(y),-cos(x)*sin(x)*sin(y)*sin(y))
p_exact = cos(x)*cos(y)
source = (-sin(x)*cos(y)-cos(x)*cos(x)*sin(2*y)+3*sin(x)*sin(x)*sin(2*y)
          ,-cos(x)*sin(y)+cos(y)*cos(y)*sin(2*x)-3*sin(y)*sin(y)*sin(2*x))

geo = SplineGeometry()
geo.AddRectangle((0,0),(pi,pi),bc="rectangle")
mesh = Mesh( geo.GenerateMesh(maxh = pi/20) )
order = 2
mesh.Curve(order)

print("ElementBoundary=",mesh.GetBoundaries())
print("NumEles=",mesh.ne)
print("NumEdges=",mesh.nedge)
print("NumVertex=",mesh.nv)

V = VectorH1(mesh,order=order, dirichlet="rectangle")
Q = H1(mesh,order=order-1)

X = FESpace([V,Q])

u,p = X.TrialFunction()
v,q = X.TestFunction()

print("V.ndof =", V.ndof,X.Range(0)) 
print("Q.ndof  =", Q.ndof,X.Range(1))
print("Total ndof=",X.ndof)

nu = 1.0 #fluid viscosity

stokes = nu*InnerProduct(grad(u), grad(v))-div(u)*q-div(v)*p - 1e-10*p*q
a = BilinearForm(X)
a += SymbolicBFI(stokes)
a.Assemble()

# Vector source term
f = LinearForm(X)   
source = CoefficientFunction(source)
f += SymbolicLFI(InnerProduct(source,v))
f.Assemble()

#Exact Dirichelt BC for velocity
gfu = GridFunction(X)

u_exact = CoefficientFunction(u_exact)
gfu.components[0].Set(u_exact, definedon=mesh.Boundaries("rectangle"))
print("Set Dirichelt velocity BC on ","rectangle")

inv_stokes = a.mat.Inverse(X.FreeDofs())

res = f.vec.CreateVector()
res.data = f.vec - a.mat*gfu.vec
gfu.vec.data += inv_stokes * res

Draw( gfu.components[0], mesh, "velocity" )
Draw( Norm(gfu.components[0]), mesh, "absvel(hdiv)")
Draw( gfu.components[1], mesh, "pressure" )

vel = CoefficientFunction(gfu.components[0])
pres = CoefficientFunction(gfu.components[1])

#Error analysis
print ("Pres L2-error:", sqrt (Integrate ( (pres-p_exact)*(pres-p_exact), mesh)))
print ("Velocity L2-error:", sqrt (Integrate ( (vel-u_exact )*(vel-u_exact ), mesh)))


#Output velocity and pressure to paraview
vtk = VTKOutput(ma=mesh,coefs=[vel,pres],names=["u","p"],filename="vtkout_CG",subdivision=0)
vtk.Do()
print("VTK file is generated!")


n = specialcf.normal(mesh.dim)
u = V.TrialFunction()

### Given a vectorial function u, compute the sum of outgoing fluxes per element
C = L2(mesh=mesh, order=0, dirichlet="rectangle")
gfc = GridFunction (space=C)
c = C.TestFunction()
a = BilinearForm(trialspace = V, testspace = C)
a += SymbolicBFI(c * u * n, element_boundary=True)
a.Apply(gfu.vec,gfc.vec)

print("Mass balance =",sum(gfc.vec))