Download this testcase.
Bidimensional Poiseuille Flow — TFEM Convergence Test¶
This example evaluates the convergence of a Two-Field Finite Element Method (TFEM) applied to a 2D Poiseuille flow. The domain is meshed using Gmsh and solved with MigFlow’s fluid module. Two mesh configurations are compared: - a regular structured mesh - a perturbed mesh containing isolated triangular “caps”.
Keywords¶
TFEM, FEM, Gmsh, Convergence, Poiseuille flow, Mesh perturbation
Description¶
This benchmark demonstrates: - How to generate a 2D rectangular channel mesh with named boundaries in Gmsh. - How to introduce local mesh perturbations (“caps”) to assess solver robustness. - How to impose an analytical Poiseuille velocity profile at the inlet. - How to solve a steady-state viscous flow using MigFlow’s TFEM formulation. - How to measure the numerical error with respect to the analytical solution.
Physical Setup¶
The flow corresponds to classical 2D Poiseuille flow driven by a pressure difference. The analytical velocity profile is:
u_x(y) = (1 / (20 μ)) * y * (1 − y)
Boundary conditions: - Left boundary: imposed parabolic inflow velocity profile. - Right boundary: zero relative pressure (open outlet). - Top/Bottom: no-slip solid walls (velocity = 0).
Numerical parameters: - ρ = 1000 kg/m³ (fluid density) - ν = 1e−3 m²/s (kinematic viscosity) - μ = ρν (dynamic viscosity) - dt = 10 (pseudo–time step) - tEnd = 1000 (final pseudo–time)
Convergence Study¶
The test loops over several mesh resolutions N = {10, 20, 40, 80, 160}. For each resolution, two simulations are run: - with caps (perturbed mesh) - without caps (regular mesh)
At the end of each simulation, the L² error norm between numerical and analytical velocity is computed:
error = || u_num − u_exact ||₂
The code then plots the error in log–log scale to assess convergence order and to compare robustness with/without mesh perturbations.
from migflow import fluid, time_integration
import numpy as np
import os
import time
import shutil
import gmsh
# Geometry
# --------
# The computational domain is a 2D channel of height H = 1 and length L = 1.
# Depending on configuration:
# - A single rectangle is generated for the regular mesh.
# - Two adjacent rectangles are generated and remeshed to insert triangular caps
# in the perturbed case.
#
# Geometrical parameters
H = 1 # domain height
L = 1 # domain width
Ox = -5
Oy = 0
# Mesh Generation
# ---------------
# Mesh generation relies on Gmsh with:
# - transfinite curves along top/bottom boundaries to control refinement;
# - automatic Delaunay meshing (Mesh.Algorithm = 6);
# - optional internal triangular elements (“caps”) that locally distort the mesh
# for convergence robustness testing.
# Meshes are exported to .msh files and loaded by MigFlow.
def mesh(
N=100, L0=0.5, D=0, withcaps=True, needles=False, one_side=False, special=False
):
L = 2 * L0
gmsh.model.mesh.set_size_callback(lambda *x: L / (N))
def make_square(x0, y0, x1, y1):
xs = [(x0, y0), (x0, y1), (x1, y1), (x1, y0)]
p = list(gmsh.model.geo.add_point(x[0], x[1], 0) for x in xs)
l = list(gmsh.model.geo.add_line(p0, p1) for (p0, p1) in zip(p, np.roll(p, 1)))
s = gmsh.model.geo.add_plane_surface([gmsh.model.geo.add_curve_loop(l)])
return s, l
if not withcaps:
s0, ls0 = make_square(0, 0, 4 * L, L)
gmsh.model.geo.synchronize()
gmsh.model.add_physical_group(1, [ls0[1]], -1, "Left")
gmsh.model.add_physical_group(1, [ls0[2]], -1, "Top")
gmsh.model.add_physical_group(1, [ls0[3]], -1, "Right")
gmsh.model.add_physical_group(1, [ls0[0]], -1, "Bottom")
gmsh.model.add_physical_group(2, [s0], -1, "domain")
gmsh.option.set_number("Mesh.Algorithm", 6)
gmsh.model.mesh.generate(2)
# gmsh.fltk.run()
return
s0, ls0 = make_square(0, 0, 4 * L0, L)
s1, ls1 = make_square(4 * L0 + D, 0, 4 * L, L)
gmsh.model.geo.synchronize()
l0 = gmsh.model.get_entities_in_bounding_box(L, 0, 0, L, L, 0, 1)
gmsh.model.mesh.set_transfinite_curve(ls0[3], N + 1)
if needles:
gmsh.model.mesh.set_transfinite_curve(ls1[1], N + 1)
elif one_side:
gmsh.model.mesh.set_transfinite_curve(ls1[1], 2 * N + 1)
elif special:
gmsh.model.mesh.set_transfinite_curve(ls1[1], 3 * N + 2)
else:
gmsh.model.mesh.set_transfinite_curve(ls1[1], N + 2)
gmsh.option.set_number("Mesh.Algorithm", 6)
gmsh.model.mesh.generate(1)
l1tag, l1, _ = gmsh.model.mesh.get_nodes(1, ls1[1], False, False)
l0tag, _, _ = gmsh.model.mesh.get_nodes(1, ls0[3], False, False)
l1 = l1.reshape(-1, 3)
if needles:
l1[:, 1] = np.linspace(L / N, L * (1 - 1 / N), N - 1)[::-1]
elif one_side:
l1[:, 1] = np.linspace(L / (2 * N), L - (L / (2 * N)), 2 * N - 1)[::-1]
elif special:
l1[1 : 3 * N : 3, 1] = np.linspace(L / (2 * N), L - (L / (2 * N)), N)[::-1]
l1[0 : 3 * N - 1 : 3, 1] = l1[1 : 3 * N : 3, 1] + D
l1[2 : 3 * N + 1 : 3, 1] = l1[1 : 3 * N : 3, 1] - D
else:
l1[:, 1] = np.linspace(L / (2 * N), L - (L / (2 * N)), N)[::-1] #
for t, x in zip(l1tag, l1):
gmsh.model.mesh.set_node(t, x, [])
# gmsh.model.mesh.synchronize()
gmsh.model.mesh.generate(2)
con = gmsh.model.add_discrete_entity(2)
l1, _, l1param = gmsh.model.mesh.get_nodes(1, ls1[1], True, True)
l0, _, l0param = gmsh.model.mesh.get_nodes(1, ls0[3], True, True)
l1 = l1[np.argsort(l1param)]
l0 = l0[np.argsort(-l0param)]
if needles:
tri = np.zeros((N, 3), np.int32)
tri[:, 1] = l1[:-1]
tri[:, 0] = l1[1:]
tri[:, 2] = l0[:-1]
gmsh.model.mesh.add_elements_by_type(con, 2, [], tri.flat)
tri2 = np.zeros((N, 3), np.int32)
tri2[:, 1] = l0[:-1]
tri2[:, 0] = l1[1:]
tri2[:, 2] = l0[1:]
gmsh.model.mesh.add_elements_by_type(con, 2, [], tri2.flat)
elif one_side:
tri = np.zeros((N, 3), np.int32)
tri[:, 0] = l1[0 : 2 * N - 1 : 2]
tri[:, 1] = l1[1 : 2 * N : 2]
tri[:, 2] = l0[:-1]
gmsh.model.mesh.add_elements_by_type(con, 2, [], tri.flat)
tri1 = np.zeros((N, 3), np.int32)
tri1[:, 0] = l1[1 : 2 * N : 2]
tri1[:, 1] = l1[2 : 2 * N + 1 : 2]
tri1[:, 2] = l0[1:]
gmsh.model.mesh.add_elements_by_type(con, 2, [], tri1.flat)
tri2 = np.zeros((N, 3), np.int32)
tri2[:, 0] = l0[:-1]
tri2[:, 1] = l1[1 : 2 * N : 2]
tri2[:, 2] = l0[1:]
gmsh.model.mesh.add_elements_by_type(con, 2, [], tri2.flat)
elif special:
tri = np.zeros((N + 1, 3), np.int32)
tri[:, 0] = l1[0 : 3 * N + 1 : 3]
tri[:, 1] = l1[1 : 3 * N + 2 : 3]
tri[:, 2] = l0[:]
gmsh.model.mesh.add_elements_by_type(con, 2, [], tri.flat)
tri1 = np.zeros((N, 3), np.int32)
tri1[:, 0] = l1[1 : 3 * N - 1 : 3]
tri1[:, 1] = l1[2 : 3 * N : 3]
tri1[:, 2] = l0[:-1]
gmsh.model.mesh.add_elements_by_type(con, 2, [], tri1.flat)
tri2 = np.zeros((N, 3), np.int32)
tri2[:, 0] = l0[:-1]
tri2[:, 1] = l1[2 : 3 * N : 3]
tri2[:, 2] = l0[1:]
gmsh.model.mesh.add_elements_by_type(con, 2, [], tri2.flat)
tri3 = np.zeros((N, 3), np.int32)
tri3[:, 0] = l1[2 : 3 * N : 3]
tri3[:, 1] = l1[3 : 3 * N + 1 : 3]
tri3[:, 2] = l0[1:]
gmsh.model.mesh.add_elements_by_type(con, 2, [], tri3.flat)
else:
tri = np.zeros((N + 1, 3), np.int32)
tri[:, 0] = l1[:-1]
tri[:, 2] = l1[1:]
tri[:, 1] = l0
gmsh.model.mesh.add_elements_by_type(con, 2, [], tri.flat)
tri2 = np.zeros((N, 3), np.int32)
tri2[:, 0] = l0[:-1]
tri2[:, 2] = l1[1:-1]
tri2[:, 1] = l0[1:]
gmsh.model.mesh.add_elements_by_type(con, 2, [], tri2.flat)
conl_top = gmsh.model.add_discrete_entity(1)
conl_bottom = gmsh.model.add_discrete_entity(1)
gmsh.model.mesh.add_elements_by_type(conl_top, 1, [], [l1[0], l0[0]])
gmsh.model.mesh.add_elements_by_type(conl_bottom, 1, [], [l0[-1], l1[-1]])
gmsh.model.add_physical_group(1, [ls0[1]], -1, "Left")
gmsh.model.add_physical_group(1, [ls0[2], ls1[2], conl_top], -1, "Top")
gmsh.model.add_physical_group(1, [ls1[3]], -1, "Right")
gmsh.model.add_physical_group(1, [ls0[0], ls1[0], conl_bottom], -1, "Bottom")
gmsh.model.add_physical_group(2, [s0, s1, con], -1, "domain")
# gmsh.fltk.run()
def mesh_isolated_caps(N=100, L=1, D=0, N_caps=-1):
gmsh.model.mesh.set_size_callback(lambda *x: L / (N))
def make_square(x0, y0, x1, y1):
xs = [(x0, y0), (x0, y1), (x1, y1), (x1, y0)]
p = list(gmsh.model.geo.add_point(x[0], x[1], 0) for x in xs)
l = list(gmsh.model.geo.add_line(p0, p1) for (p0, p1) in zip(p, np.roll(p, 1)))
s = gmsh.model.geo.add_plane_surface([gmsh.model.geo.add_curve_loop(l)])
return s, l
s0, ls0 = make_square(0, 0, 4 * L, L)
gmsh.model.geo.synchronize()
gmsh.option.set_number("Mesh.Algorithm", 6)
gmsh.model.mesh.generate(2)
if N_caps < 0:
N_caps = (N // 2) * 5
nodestag, nodes, nodesparam = gmsh.model.mesh.get_nodes()
nodesorder = np.argsort(nodestag)
nodes = nodes.reshape(-1, 3)[:, :]
[etype], [etags], [elements] = gmsh.model.mesh.get_elements(2)
elements = elements.reshape(etags.shape[0], -1)
n = elements.shape[1]
elements = np.searchsorted(nodestag, elements, sorter=nodesorder)
_, bnd_edges_top = gmsh.model.mesh.get_elements_by_type(1, ls0[2])
_, bnd_edges_bottom = gmsh.model.mesh.get_elements_by_type(1, ls0[0])
_, bnd_edges_left = gmsh.model.mesh.get_elements_by_type(1, ls0[1])
_, bnd_edges_right = gmsh.model.mesh.get_elements_by_type(1, ls0[3])
random_iel = np.random.choice(elements.shape[0], N_caps, replace=False)
random_el = elements[random_iel]
new_nodes = np.mean(nodes[random_el[:, :2]], axis=1)
dists = np.linalg.norm(nodes[random_el[:, 2]] - new_nodes, axis=1)[:, None]
Ds = np.minimum(D, dists / 2)
new_nodes += Ds * (nodes[random_el[:, 2]] - new_nodes) / dists
final_nodes = np.vstack([nodes, new_nodes])
new_ids = nodes.shape[0] + np.arange(N_caps)
new_elements = np.array([random_el[:, 0], random_el[:, 1], new_ids])
new_elements2 = np.array([random_el[:, 1], random_el[:, 2], new_ids])
new_elements3 = np.array([random_el[:, 2], random_el[:, 0], new_ids])
new_elements = np.vstack([new_elements.T, new_elements2.T, new_elements3.T])
elements = np.delete(elements, random_iel, axis=0)
final_elements = np.vstack([elements, new_elements]) + 1
gmsh.model.mesh.clear()
# Inject transformed mesh
gmsh.model.mesh.add_nodes(
2, 1, np.arange(final_nodes.shape[0]) + 1, final_nodes.flatten()
)
# modified
dom_entity = gmsh.model.add_discrete_entity(2)
gmsh.model.mesh.add_elements_by_type(
dom_entity, 2, np.arange(final_elements.shape[0]) + 1, final_elements.flatten()
)
bnd_entity_top = gmsh.model.add_discrete_entity(1)
gmsh.model.mesh.add_elements_by_type(bnd_entity_top, 1, [], bnd_edges_top.flatten())
gmsh.model.add_physical_group(1, [bnd_entity_top], 10, "Top")
bnd_entity_bottom = gmsh.model.add_discrete_entity(1)
bnd_entity_left = gmsh.model.add_discrete_entity(1)
bnd_entity_right = gmsh.model.add_discrete_entity(1)
gmsh.model.mesh.add_elements_by_type(
bnd_entity_bottom, 1, [], bnd_edges_bottom.flatten()
)
gmsh.model.add_physical_group(1, [bnd_entity_bottom], 11, "Bottom")
gmsh.model.mesh.add_elements_by_type(
bnd_entity_left, 1, [], bnd_edges_left.flatten()
)
gmsh.model.add_physical_group(1, [bnd_entity_left], 12, "Left")
gmsh.model.mesh.add_elements_by_type(
bnd_entity_right, 1, [], bnd_edges_right.flatten()
)
gmsh.model.add_physical_group(1, [bnd_entity_right], 13, "Right")
gmsh.model.add_physical_group(2, [dom_entity], -1, "domain")
Run the simulation for different mesh configurations¶
def run(N, withcaps):
# Output Directory
# ----------------
# Create a clean output directory for simulation results.
outputdir = "output"
shutil.rmtree(outputdir, True)
if not os.path.isdir(outputdir):
os.makedirs(outputdir)
# Mesh Generation
# ---------------
# Depending on the `withcaps` flag, either a regular mesh or a perturbed
# mesh with triangular caps is generated and saved to .msh files.
if withcaps:
mesh_isolated_caps(N=N)
filename = outputdir + f"/mesh_{N}_caps.msh"
else:
mesh(N=N, D=0, withcaps=False)
filename = outputdir + f"/mesh_{N}_nocaps.msh"
gmsh.write(filename)
# Physical parameters
# -------------------
g = 0 # gravity
rho = 1000 # fluid density
nu = 1e-3 # kinematic viscosity
mu = nu * rho # dynamic viscosity
# Numerical parameters
# --------------------
#
# FLUID PROBLEM
#
# Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
f = fluid.FluidProblem(2, g, mu, rho)
f.load_msh(filename)
f.set_open_boundary(
"Left", velocity=[lambda x: 1 / (20 * mu) * x[:, 1] * (1 - x[:, 1]), 0]
)
f.set_wall_boundary("Bottom", velocity=[0, 0])
f.set_wall_boundary("Top", velocity=[0, 0])
f.set_strong_boundary("Bottom", velocity=[0, 0])
f.set_strong_boundary("Top", velocity=[0, 0])
f.set_open_boundary("Right", pressure=0)
# Simulation Parameters
# ---------------------
t = 0 # initial time
i = 0 # iteration counter
outf = 1 # number of iterations between output files
dt = 10 # time step
tEnd = 500 # final time
# Simulation Loop
# ---------------
while t < tEnd:
print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
time_integration.iterate(f, None, dt)
t += dt
# Output files writting
if i % outf == 0:
f.write_mig(outputdir, t)
i += 1
# Error Computation
# -----------------
# After reaching the final time, the error at the nodes is computed between the numerical
# and analytical velocity profiles to assess convergence.
s = f.velocity()
x = f.coordinates_fields()[f.velocity_index()[:, 0]]
vel = (s[:, 0] - 1 / (20 * nu * rho) * x[:, 1] * (1 - x[:, 1])) ** 2
vS = np.sum((1 / (20 * nu * rho) * x[:, 1] * (1 - x[:, 1])) ** 2)
print(f"Error criteria : {(vel.sum())**.5} \t < {(vS**0.5)/50}")
print(f"Succeed ? {(vel.sum())**.5 < (vS**0.5)/50}")
gmsh.model.remove()
return vel.sum() ** 0.5
Main Convergence Loop¶
The main script loops over different mesh resolutions and both mesh configurations (with/without caps), collecting error norms for convergence analysis.
if __name__ == "__main__":
Ns = np.array([10, 20, 40, 80, 160])
data_caps = []
data_nocaps = []
plot = False
for N in Ns:
for withcaps in [False, True]:
print(f"N = {N}, withcaps = {withcaps}")
sol = run(N, withcaps)
if withcaps:
data_caps.append(sol)
else:
data_nocaps.append(sol)
if plot:
import matplotlib.pyplot as plt
plt.figure()
# log log plot
plt.loglog(1.0 / Ns, data_caps, label="with caps", marker="o")
plt.loglog(1.0 / Ns, data_nocaps, label="without caps", marker="o")
plt.legend()
plt.show()
Plot¶
python3 -m migflow.plot.migplot output --actors fluid --fluid-field velocity --fluid-vmin 0 --fluid-vmax 0.0125