Help with Periodic Taylor-Green Vortex Simulation — Pressure Remains Constant / Large Pressure Error

General
1

Hello everyone,

I am working on a 2D periodic Taylor-Green Vortex problem using NGSolve and trying to validate my finite element Navier–Stokes solver. My setup is:

  • Domain: [0,1]×[0,1][0,1]\times[0,1][0,1]×[0,1], with fully periodic boundary conditions.
  • Viscosity: ν=1/Re\nu = 1/\text{Re}ν=1/Re with Re=100\text{Re} = 100Re=100.
  • Finite Element Spaces:
    • Velocity: V = Periodic(VectorH1(mesh, order=k))
    • Pressure: Q = Periodic(H1(mesh, order=k-1))
    • Combined: X = FESpace([V, Q])

I am attempting a time-stepping approach where I assemble a bilinear form a that includes velocity and pressure terms:

python

复制代码

a = BilinearForm(
    (InnerProduct(u, v)
     + dt*(nu*InnerProduct(grad(u), grad(v)) 
            - div(v)*p 
            - div(u)*q))
    * dx
)

Then I create its inverse (inv = a.mat.Inverse(X.FreeDofs())) and perform partial updates for the velocity and pressure in each time step. However, I notice that my pressure solution effectively remains constant or has large errors compared to the analytical Taylor-Green pressure.

What I’ve Tried

  1. Removing the mean pressure:
    Since I know that in a fully periodic domain, the pressure is only determined up to an additive constant, I attempted to remove the global mean of the pressure after each time step:

python

复制代码

pgrid = gfu.components[1]
mean_val = Integrate(pgrid, mesh) / Integrate(CoefficientFunction(1), mesh)
pgrid.Set(pgrid - mean_val)

This does shift the pressure field so that its integral is zero, but the overall pressure error remains large (and sometimes increases for smaller time steps).
2. SetMeanValue(0.0) approach:
I tried to call Q.SetMeanValue(0.0) on the periodic pressure space but got:

plaintext

复制代码

AttributeError: 'ngsolve.comp.Periodic' object has no attribute 'SetMeanValue'

So it seems that this function is not exposed for a Periodic(...)-wrapped space in my NGSolve version.
3. Verification with L2 errors:

  • The velocity converges reasonably well with decreasing dt (showing about first-order in time, as expected for my approach).
  • The pressure, however, exhibits high or inconsistent errors. For example:

csharp

复制代码

dt = 0.10000, L2 Error (pressure) = 3.82e-01
dt = 0.05000, L2 Error (pressure) = 2.22e+00
dt = 0.02500, L2 Error (pressure) = 2.22e+00
dt = 0.01250, L2 Error (pressure) = 2.50e-01

The pattern is somewhat erratic, which suggests my time-stepping algorithm is not properly enforcing a pressure correction in each step.

My Questions

  1. How should I properly impose a zero-mean (or some form of constraint) on the pressure in a fully periodic problem in NGSolve?
  • Is there a recommended way to define a “mean-free” space for the pressure or add a constraint through a Lagrange multiplier in the saddle-point problem?
  1. Is my partial-step approach for velocity and pressure appropriate for an incompressible Navier–Stokes problem with periodic boundaries?
  • Right now, I update gfu1, gfu2, gfu3 and then combine them into a final solution using certain coefficients (alpha and beta), but it seems the pressure is never explicitly “coupled” at each step.
  1. Any suggestions on a stable time discretization (e.g., projection methods or fully coupled solves) that are known to handle periodic incompressible flow and yield accurate pressure fields?

Codes
from ngsolve import *
from netgen.geom2d import *
import numpy as np
import matplotlib.pyplot as plt

定义参数

k = 2
Re = 50
nu = 1 / Re # 动力黏性系数
dt_list = [0.1, 0.05, 0.025, 0.0125, 0.00625] # 选择不同的时间步长
tend = 1 # 计算时间

创建网格和有限元空间

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/64))

V = Periodic(VectorH1(mesh, order=k))
Q = Periodic(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(2pix)cos(2piy)decay, -cos(2pix)sin(2pi*y)*decay))

存储误差

errors =

遍历不同的时间步长

for dt in dt_list:
print(f"Computing for dt = {dt}")

# 初始速度
gfu = GridFunction(X)
velocity, pressure = gfu.components
velocity.Set(taylor_green_velocity(0))

# 定义 bilinear 和 linear forms
a = BilinearForm((nu * InnerProduct(grad(u),grad(v)) - div(v) * p - div(u) * q) * dx).Assemble()
mstar = BilinearForm(X)
mstar += InnerProduct(u, v) * dx  # 质量项
mstar += dt * nu * InnerProduct(grad(u), grad(v)) * dx  # 黏性扩散项
mstar += dt * (-div(v) * p * dx - div(u) * q * dx)  # 连续性约束（压力-速度耦合）
mstar.Assemble()

r = BilinearForm((InnerProduct(u,v) - dt * (Grad(u) * u) * v) * dx).Assemble()
inv = mstar.mat.Inverse(X.FreeDofs())
f = LinearForm(X)

t = 0
while t < tend - 0.5*dt:
    f.vec.data = r.mat * gfu.vec
    gfu.vec.data = inv * f.vec
    t += dt

# 计算 L2 误差
u_exact = taylor_green_velocity(tend)
l2_error = sqrt(Integrate(InnerProduct(velocity - u_exact, velocity - u_exact), mesh))
errors.append(l2_error)

计算误差阶数（时间收敛阶）

orders =
for i in range(len(dt_list) - 1):
order = np.log(errors[i] / errors[i + 1]) / np.log(dt_list[i] / dt_list[i + 1])
orders.append(order)

打印误差和收敛阶数

for i in range(len(dt_list)):
print(f"dt = {dt_list[i]:.5f}, Error = {errors[i]:.6e}")

for i in range(len(orders)):
print(f"Order between dt={dt_list[i]} and dt={dt_list[i+1]}: {orders[i]:.2f}")

绘制误差对 dt 的对数图

plt.figure()
plt.loglog(dt_list, errors, marker=“o”, linestyle=“–”, label=“L2 Error”)
plt.xlabel(“Time Step dt”)
plt.ylabel(“L2 Error”)
plt.legend()
plt.grid(True)
plt.show()