- Thank you received: 0
Forum Message
We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.
Notice
The forum is in read only mode.
VTKOutput() not works for stokes problem with Tayler-Hood
5 years 1 month ago #1482
by niushi200
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:
Please check attached script to reproduce the issue.
Thanks,
Bin
But ngsolve seems stucks on export VTK file and prohibit me to load other python script with the following error:
Code:
ERROR: Thread already running
errinfo: Thread already running
while executing
"NGS_LoadPy $file"
(menu invoke)
Attachment not found
Please check attached script to reproduce the issue.
Thanks,
Bin
Attachments:
5 years 1 month ago #1483
by schruste
Replied by schruste on topic VTKOutput() not works for stokes problem with Tayler-Hood
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!)).
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!)).
5 years 1 month ago - 5 years 1 month ago #1484
by niushi200
Replied by niushi200 on topic VTKOutput() not works for stokes problem with Tayler-Hood
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:
Here is the python code I used for CG
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:
Code:
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
Code:
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))
Last edit: 5 years 1 month ago by niushi200.
Time to create page: 0.111 seconds
