Download this testcase.

Advection–Diffusion or Fluid Transport in a 2D Periodic Channel

This example illustrates how to solve a 2D periodic channel problem using either the advection–diffusion equation (ADV_DIFF = 1) or the incompressible Navier–Stokes equations (ADV_DIFF = 0).

The domain is a rectangular channel with periodic lateral boundaries and Dirichlet conditions at the top and bottom.

Keywords

Advection-Diffusion, Periodic Boundaries, Mesh Adaptation

import os, shutil, gmsh, sys
import numpy as np
from migflow import fluid, advdiff

Output Directory

Create a clean output directory for simulation results.

outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)

Problem Configuration

ADV_DIFF = 1  # 1 → advection–diffusion, 0 → fluid simulation

Helper Class for Mesh Sizing

class MeshSize:
    """Callable mesh-size function for gmsh."""

    def __init__(self, lc):
        self._lc = lc

    def __call__(self, dim, tag, x, y, z, lc):
        return self._lc

Geometrical parameters

width, height = 4.0, 1.0
origin = [-width / 2, -height / 2]
lc = 0.1
mesh_size = MeshSize(lc)

Physical Parameters

g = np.array([0, 0])  # gravity (ignored for adv-diff)
rho = 1.204  # nitrogen density [kg/m³]
nu = 1.8e-5  # kinematic viscosity [m²/s]
mu = rho * nu  # dynamic viscosity [kg/(m·s)]
alpha = 3e-3  # thermal dilatation coefficient
k = 2e-2  # thermal conductivity [W/mK]
cp = 1010e-3  # heat capacity [kJ/kg·K]

Generate Initial Mesh

def gen_mesh(width, height, mesh_size, origin=np.array([0, 0])):
    """Generate a rectangular periodic mesh in gmsh.

    Parameters
    ----------
    width : float
        Domain width.
    height : float
        Domain height.
    mesh_size : MeshSize
        Callable object controlling element size.
    origin : np.ndarray, optional
        Origin of the rectangle.
    """
    origin = np.asarray(origin)
    gmsh.model.add("box")
    gmsh.model.occ.add_rectangle(origin[0], origin[1], 0, width, height)
    gmsh.model.occ.synchronize()

    def get_line(x0, x1, eps=1e-6):
        r = gmsh.model.get_entities_in_bounding_box(
            x0[0] - eps, x0[1] - eps, -eps, x1[0] + eps, x1[1] + eps, eps, 1
        )
        return [tag for _, tag in r]

    h, w = height, width
    l_tag = get_line(origin + [0, 0], origin + [0, h])
    r_tag = get_line(origin + [w, 0], origin + [w, h])
    # Define periodicity (right ↔ left)
    affine_matrix = np.eye(4, 4, dtype=np.float64)
    affine_matrix[0, 3] = width
    gmsh.model.mesh.set_periodic(1, r_tag, l_tag, affine_matrix.reshape(-1))
    gmsh.model.occ.synchronize()
    # Define physical groups
    gmsh.model.add_physical_group(
        1, get_line(origin + [0, 0], origin + [w, 0]), name="Bottom"
    )
    gmsh.model.add_physical_group(
        1, get_line(origin + [0, h], origin + [w, h]), name="Top"
    )
    gmsh.model.add_physical_group(2, [1], name="domain")
    gmsh.model.mesh.set_size_callback(mesh_size)


gmsh.initialize()
gen_mesh(width, height, mesh_size, origin)
gmsh.model.mesh.generate(2)

Problem Definition

def define_problem(adv_diff):
    """Create either an advection–diffusion or fluid problem."""
    if adv_diff:
        # --- Advection–Diffusion problem ---
        c = advdiff.AdvDiffProblem(2, 0, 0, 1.0, explicit=True, cg=False)
        c.load_msh(None)
        c.set_weak_boundary("Bottom")
        c.set_weak_boundary("Top")
        c.velocity()[:, 0] = 1.0  # uniform flow to the right
        x = c.coordinates_fields()
        c.solution()[x[:, 0] > -1e-8] = 1  # initialize concentration at inflow
        return c
    else:
        # --- Fluid problem (reference) ---
        f = fluid.FluidProblem(2, g, [mu, mu], [rho, rho])
        f.load_msh(None)
        f.set_mean_pressure(0)
        f.set_wall_boundary("Bottom")
        f.set_wall_boundary("Top")
        s0 = np.zeros_like(f.solution())
        s0[f.velocity_index()[:, 0]] = 1.0  # initial velocity in x-direction
        f.solution()[:] = s0
        c0 = np.zeros_like(f.coordinates()[:, 0])
        xc0 = f.coordinates()
        c0[xc0[:, 0] > -1e-8] = 1
        f.set_concentration_cg(c0)
        return f


p = define_problem(ADV_DIFF)

Mesh Adaptation

mesh_size = MeshSize(0.05)
gen_mesh(width, height, mesh_size, origin)
gmsh.model.mesh.clear()
gmsh.model.mesh.generate(2)
if ADV_DIFF:
    p.adapt_mesh()
    p.velocity()[:, 0] = 1.0
else:
    p.adapt_mesh()

Simulation Parameters

dt = 0.01  # time step [s]
tEnd = 1.0  # final time [s]
outf = 1  # output frequency

Time Integration Loop

t, i = 0, 0
while t < tEnd:
    print(f"{i:4d}, {t:.6g}/{tEnd:.6g}")
    if i % outf == 0:
        p.write_mig(outputdir, t)
    if ADV_DIFF:
        p.explicit_euler(dt, limiter=True)
    else:
        p.advance_concentration(dt)
    t += dt
    i += 1

Plot

python3 -m migflow.plot.migplot output --actors fluid --problem-type advection-diffusion --fluid-field solution --fluid-vmin 0 --fluid-vmax 0.25 --show-edges 1 --element-type triangle_p1dg