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.
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))