High L2 Error in Taylor-Green Vortex (Navier-Stokes) Using FEM with NGSolve
Hi all,
I’m currently implementing a finite element solver for the 2D incompressible Navier-Stokes equations using NGSolve and testing it against the classical Taylor-Green vortex benchmark. While the setup runs without any errors, I’m observing unexpectedly large L2 errors in both velocity and pressure, even after only one time step using explicit Euler.
Problem Summary:
- PDE: 2D incompressible Navier-Stokes
- Test case: Taylor-Green vortex (analytical solution known)
- Discretization:
- Finite element spaces:
VectorH1(order=k)andH1(order=k-1)withk=2 - Time-stepping: Explicit Euler
- Periodic boundary conditions on all sides
- Mesh curved with
mesh.Curve(3)
- Finite element spaces:
Unexpected Results:
Even with a small time step (dt = 0.001, final time T = 1.0), the output is:
java
复制代码
===== L2 Error Results (dt = 0.001) =====
L2 Velocity Error = 1.61e-01
L2 Pressure Error = 3.79e-01
These errors seem too large for such a simple benchmark problem.
Key Code Snippets:
python
复制代码
# Time-stepping loop
t = 0
while t < T - 0.5 * dt:
conv.Assemble()
res = a.mat * gfu.vec + conv.vec
gfu.vec.data -= dt * inv * res
t += dt
python
复制代码
# L2 error computation
velocity_error = sqrt(Integrate(InnerProduct(gfu.components[0] - velocity_exact,
gfu.components[0] - velocity_exact), mesh))
pressure_error = sqrt(Integrate((gfu.components[1] - pressure_exact)**2, mesh))
What I Have Tried:
Initial conditions match the analytical solution at t=0- Pressure normalization to eliminate constant offset
- Periodic boundary conditions using
copyandbc="periodic" - Mesh curvature enabled (
mesh.Curve(3)) - Reduced
dtto1e-4— error still significant - Tried both assembled and non-assembled versions of convective term
Main Questions:
- Is my implementation of the convective term correct?
I’m usinggrad(u)*u, but not sure if this is sufficient or needs skew-symmetric form/stabilization? - Is explicit Euler a poor choice for this Reynolds number (Re=100)?
Should I switch to an implicit method even for small dt? - Are the chosen FE spaces appropriate for velocity-pressure pairing?
UsingVectorH1andH1with orderskandk-1, respectively. - Do I need any pressure stabilization terms in this setup?
Full Code:
python
复制代码
from ngsolve import *
from netgen.geom2d import *
from ngsolve.webgui import Draw
import numpy as np
# 定义参数
k = 2
Re = 100
nu = 1 / Re
T = 1.0
dt = 0.001 # 仅一个时间步长
# 创建网格和有限元空间（周期性）
geo = SplineGeometry()
points = [(0,0), (1,0), (1,1), (0,1)]
pnums = [geo.AppendPoint(*p) for p in points]
bottom = geo.Append(["line", pnums[0], pnums[1]], leftdomain=1, rightdomain=0, bc="periodic")
right = geo.Append(["line", pnums[1], pnums[2]], leftdomain=1, rightdomain=0, bc="periodic")
top = geo.Append(["line", pnums[3], pnums[2]], leftdomain=0, rightdomain=1, copy=bottom, bc="periodic")
left = geo.Append(["line", pnums[0], pnums[3]], leftdomain=0, rightdomain=1, copy=right, bc="periodic")
mesh = Mesh(geo.GenerateMesh(maxh=1/32))
mesh.Curve(3)
V = VectorH1(mesh, order=k)
Q = H1(mesh, order=k-1)
X = FESpace([V, Q])
# 函数定义
u, p = X.TrialFunction()
v, q = X.TestFunction()
def taylor_green_velocity(t):
decay = exp(-8 * pi**2 * nu * t)
return CoefficientFunction((sin(2*pi*x)*cos(2*pi*y)*decay, -cos(2*pi*x)*sin(2*pi*y)*decay))
def taylor_green_pressure(t):
decay = exp(-16 * pi**2 * nu * t)
return 0.25 * (cos(4*pi*x) + cos(4*pi*y)) * decay
def NormalizePressure(gfu, mesh):
p = gfu.components[1]
mean_p = Integrate(p, mesh) / Integrate(1, mesh)
# 构造与 p.vec.data 同类型、同大小的常数向量
ones = p.vec.CreateVector()
ones[:] = 1.0
ones *= mean_p
# 执行减法：p -= mean_p * 1
p.vec.data -= ones
# 初始条件
gfu = GridFunction(X)
velocity, pressure = gfu.components # X由V，P组成，gfu.components即为[V，P]
gfu.components[0].Set(taylor_green_velocity(0))
gfu.components[1].Set(taylor_green_pressure(0))
# 系统矩阵定义（Euler）
a = BilinearForm((nu*InnerProduct(grad(u),grad(v)) - div(v)*p - div(u)*q)*dx).Assemble()
mstar = BilinearForm(X, symmetric=True)
mstar += InnerProduct(u, v) * dx
mstar += dt * nu * InnerProduct(grad(u), grad(v)) * dx
mstar += dt * (-div(v) * p * dx - q * div(u) * dx)
mstar.Assemble()
# conv = BilinearForm(X, nonassemble=True)
# conv += InnerProduct(grad(u)*u, v) * dx
conv = LinearForm(X)
conv += InnerProduct(grad(velocity)*velocity,v)*dx
inv = mstar.mat.Inverse(X.FreeDofs())
# 时间推进
t = 0
while t < T - 0.5 * dt:
conv.Assemble()
# res = (a.mat + conv.mat) * gfu.vec
res = a.mat * gfu.vec + conv.vec
gfu.vec.data -= dt * inv * res
t += dt
if int(t*10000) % 1000 == 0:
print(f"t = {t:.4f}")
NormalizePressure(gfu, mesh)
# 误差计算
velocity_exact = taylor_green_velocity(T)
pressure_exact = taylor_green_pressure(T)
velocity_error = sqrt(Integrate(InnerProduct(gfu.components[0] - velocity_exact,
gfu.components[0] - velocity_exact), mesh))
pressure_error = sqrt(Integrate((gfu.components[1] - pressure_exact)**2, mesh))
print("\n===== L2 误差结果 (dt = 0.0001) =====")
print(f"L2 Velocity Error = {velocity_error:.6e}")
print(f"L2 Pressure Error = {pressure_error:.6e}")
# 提取压力分量
p_num = gfu.components[1]
# 计算压力在整个域上的积分均值
total_volume = Integrate(1, mesh)
mean_pressure = Integrate(p_num, mesh) / total_volume
Any insight into why this setup gives such large errors would be greatly appreciated. I’m especially looking for advice on time integration schemes, correct form of the nonlinear term, and anything related to FE space compatibility for incompressible flow.
Thanks a lot in advance!
— A curious NGSolve user trying to validate Taylor-Green Vortex 