Example 1: MMS on the Unit Square

Problem Setup¶

We verify the DG advection-diffusion implementation using the Method of Manufactured Solutions (MMS) on a unit square domain $\Omega = [0,1]^2$ with a prescribed uniform upward velocity field $\mathbf{w} = (0, w_\text{mag})$ .

The governing equation is:

\nabla \cdot (\mathbf{w} u) - \nabla \cdot (D \nabla u) = f \quad \text{in } \Omega

(1)

Boundary Conditions¶

All four sides carry a Dirichlet condition for the diffusion term ( $\Gamma_D = \partial\Omega$ ), imposed weakly via SIPG. The advection inflow/outflow split is determined by $\mathbf{w}\cdot\mathbf{n}$ :

Boundary	Location	$\mathbf{w}\cdot\mathbf{n}$	Advection type
Bottom	$y = 0$	$-w_\text{mag} < 0$	Inflow — $u = u_\text{exact}$ prescribed
Top	$y = 1$	$+w_\text{mag} > 0$	Outflow — interior trace used
Left / Right	$x = 0,\, 1$	0	Wall — no advective flux

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.patches as patches

matplotlib.rcParams['svg.fonttype'] = 'none'

fig, ax = plt.subplots(figsize=(3.5, 3.5))
ax.add_patch(patches.Rectangle((0, 0), 1, 1, facecolor='white', edgecolor='black', linewidth=1.5, zorder=1))
ax.text(0.5, 0.5, r'$\Omega$', ha='center', va='center', fontsize=18)
ax.text(0.5, -0.06, r'$\Gamma_D$ (inflow)',  ha='center', va='top',    fontsize=10)
ax.text(0.5,  1.06, r'$\Gamma_D$ (outflow)', ha='center', va='bottom', fontsize=10)
ax.text(-0.06, 0.5, r'$\Gamma_D$ (wall)',    ha='right',  va='center', fontsize=10)
ax.text( 1.06, 0.5, r'$\Gamma_D$ (wall)',    ha='left',   va='center', fontsize=10)
for xi in [0.3, 0.7]:
    for yi in [0.2, 0.5, 0.78]:
        ax.annotate('', xy=(xi, yi + 0.11), xytext=(xi, yi),
                    arrowprops=dict(arrowstyle='->', color='black', lw=1.2), zorder=2)
ax.text(0.82, 0.88, r'$\mathbf{w}$', fontsize=13, ha='center')
ax.set_xlim(-0.35, 1.35); ax.set_ylim(-0.18, 1.18)
ax.set_aspect('equal'); ax.axis('off')
plt.tight_layout()
fig.savefig('../images/domain_mwe1.svg', bbox_inches='tight')
plt.close(fig)

Manufactured Solution¶

The exact solution is chosen to contain both a smooth oscillatory component and an exponential boundary layer near the outflow boundary ( $y = 1$ ):

u_\text{exact}(x, y) = \cos(\pi x) \cdot \frac{1 - e^{(y-1)/D}}{1 - e^{-2/D}} + \frac{1}{2}\cos(\pi x)\sin(\pi y)

(2)

The first term produces a layer of thickness $\mathcal{O}(D)$ near $y=1$ ; for small $D$ (large Péclet number) this layer becomes sharp. The source term $f$ is derived analytically from $u_\text{exact}$ so that it satisfies the PDE exactly.

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import ufl

matplotlib.rcParams['svg.fonttype'] = 'none'

D = 0.1

def exact_solution(x, D):
    """UFL exact solution, used inside the weak form."""
    return (
        ufl.cos(ufl.pi * x[0]) * (1 - ufl.exp((x[1] - 1) / D)) / (1 - ufl.exp(-2 / D))
        + 0.5 * ufl.cos(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1])
    )

def u_exact_np(x, y, D):
    """NumPy exact solution, used for plotting."""
    return (
        np.cos(np.pi * x) * (1 - np.exp((y - 1) / D)) / (1 - np.exp(-2 / D))
        + 0.5 * np.cos(np.pi * x) * np.sin(np.pi * y)
    )

x = np.linspace(0, 1, 300)
y = np.linspace(0, 1, 300)
X, Y = np.meshgrid(x, y)
U = u_exact_np(X, Y, D)

fig, ax = plt.subplots(figsize=(4.5, 4))
cf = ax.pcolormesh(X, Y, U, cmap='RdBu_r', shading='auto',
                   edgecolors='face', rasterized=True)
cbar = plt.colorbar(cf, ax=ax, shrink=0.9)
cbar.set_label(r'$u_\mathrm{exact}$', fontsize=10)
ax.set_aspect('equal')
ax.set_xlabel(r'$x$', fontsize=11)
ax.set_ylabel(r'$y$', fontsize=11)
ax.set_title(r'Manufactured solution ($D = 0.1$)', fontsize=11)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.tight_layout()
fig.savefig('../images/solution_mwe1.svg', bbox_inches='tight', dpi=200)
plt.close(fig)

Problem formulation¶

We build the unit square mesh, define the DG function space and test function, and set the physical parameters: $D = 0.1$ , $w_\text{mag} = 1.0$ . The relevant measure of the balance between advection and diffusion is the mesh Péclet number $\text{Pe} = |\mathbf{w}|\,h\,/\,(2D)$ , where $h$ is the local mesh size. For a $64 \times 64$ mesh on the unit square, $h \approx 1/64$ , giving $\text{Pe} \approx 0.08$ , so this is a diffusion-dominated case. The penalty parameter is set to $\alpha = 10 p^2$ .

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from dolfinx import fem, mesh, geometry
from dolfinx.fem.petsc import NonlinearProblem

N = 64
degree = 1
D_value = 0.1
w_mag   = 1.0

msh = mesh.create_unit_square(
    MPI.COMM_WORLD, N, N,
    cell_type=mesh.CellType.triangle,
    diagonal=mesh.DiagonalType.crossed,
)

V = fem.functionspace(msh, ("DG", degree))
u = fem.Function(V)
v = ufl.TestFunction(V)

x = ufl.SpatialCoordinate(msh)
n = ufl.FacetNormal(msh)
h = ufl.CellDiameter(msh)

D       = fem.Constant(msh, PETSc.ScalarType(D_value))
w       = fem.Constant(msh, PETSc.ScalarType((0, w_mag)))
penalty = fem.Constant(msh, PETSc.ScalarType(10.0 * degree**2))

u_ex = exact_solution(x, D)
f    = -ufl.div(D * ufl.grad(u_ex)) + ufl.dot(w, ufl.grad(u_ex))

dx, ds, dS = ufl.dx, ufl.ds, ufl.dS

The indicator function $\lambda$ determines whether a given face is an outflow or inflow/wall face with respect to the advective velocity:

\lambda = \begin{cases} 1 & \text{if } \mathbf{w} \cdot \mathbf{n} > 0 \quad (\text{outflow}) \\ 0 & \text{otherwise} \quad (\text{inflow or wall}) \end{cases}

(3)

This allows the upwind advective flux to be written compactly. On an interior face $F$ with normal $\mathbf{n} = \mathbf{n}_0$ , the combination $2\langle \lambda \mathbf{w} u \rangle \cdot \mathbf{n}$ reduces to the standard upwind flux $( \mathbf{w} \cdot \mathbf{n})\, u_\text{up}$ , where $u_\text{up}$ is the value from the upwind element.

On the boundary $\partial\Omega$ , $\lambda$ is used to split the advective boundary contribution: the outflow term $\lambda(\mathbf{w} \cdot \mathbf{n})u$ enters the bilinear form using the interior trace, while the inflow term $(1 - \lambda)(\mathbf{w} \cdot \mathbf{n})u_\text{exact}$ prescribes the inflow data on the right-hand side.

lmbda   = ufl.conditional(ufl.gt(ufl.dot(w, n), 0), 1, 0)

Advection¶

The advection contribution groups the bulk volume term, the interior upwind flux, and both boundary cases together:

\begin{align} F_\text{adv}(u;\,v) &= \underbrace{-\int_\Omega (\mathbf{w} u) \cdot \nabla v \, \mathrm{d}x}_{\text{bulk}} \\ &\quad + \underbrace{\int_{\mathcal{F}_I} (\mathbf{w} \cdot \mathbf{n})\, u_\text{up}\, \llbracket v \rrbracket \, \mathrm{d}s}_{\text{interior upwind flux}} \\ &\quad + \underbrace{\int_{\partial\Omega} \lambda\, (\mathbf{w} \cdot \mathbf{n})\, u\, v \, \mathrm{d}s}_{\text{outflow BC}} \\ &\quad - \underbrace{\int_{\partial\Omega} (1 - \lambda)\, (\mathbf{w} \cdot \mathbf{n})\, u_\text{exact}\, v \, \mathrm{d}s}_{\text{inflow BC}} \end{align}

(4)

F = 0

# Advection, bulk
F += -ufl.inner(w * u, ufl.grad(v)) * dx

# Advection, interior upwind flux
F += ufl.inner(2 * ufl.avg(lmbda * w * u), ufl.jump(v, n)) * dS

# Advection, outflow BC (interior trace)
F += ufl.inner(lmbda * ufl.dot(w, n) * u, v) * ds

# Advection, inflow BC (prescribed data)
F += -ufl.inner((1 - lmbda) * ufl.dot(w, n) * u_ex, v) * ds

Diffusion (SIPG)¶

The symmetric interior penalty formulation consists of the bulk gradient term, two face terms to restore consistency and symmetry, and an inter-element penalty:

\begin{align} F_\text{diff}(u;\,v) &= \underbrace{D \int_\Omega \nabla u \cdot \nabla v \, \mathrm{d}x}_{\text{bulk}} \\ &\quad - \underbrace{D \int_{\mathcal{F}_I} \langle \nabla u \rangle \cdot \mathbf{n}\, \llbracket v \rrbracket \, \mathrm{d}s}_{\text{consistency}} \\ &\quad - \underbrace{D \int_{\mathcal{F}_I} \llbracket u \rrbracket\, \langle \nabla v \rangle \cdot \mathbf{n} \, \mathrm{d}s}_{\text{symmetry}} \\ &\quad + \underbrace{\dfrac{\alpha D}{\langle h \rangle} \int_{\mathcal{F}_I} \llbracket u \rrbracket \llbracket v \rrbracket \, \mathrm{d}s}_{\text{penalty}} \end{align}

(5)

# Diffusion, bulk
F += D * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx

# Diffusion, SIPG interior faces
F += -D * ufl.inner(ufl.avg(ufl.grad(u)), ufl.jump(v, n)) * dS
F += -D * ufl.inner(ufl.jump(u, n), ufl.avg(ufl.grad(v))) * dS
F += D * (penalty / ufl.avg(h)) * ufl.inner(ufl.jump(u, n), ufl.jump(v, n)) * dS

Weak Dirichlet Condition and Source¶

The Dirichlet condition $u = u_\text{exact}$ on $\partial\Omega$ is imposed weakly by mirroring the SIPG treatment at boundary faces. The manufactured source closes the residual:

\begin{align} F_\text{Dir}(u;\,v) &= \underbrace{-D \int_{\partial\Omega} (\nabla u \cdot \mathbf{n})\, v \, \mathrm{d}s}_{\text{flux consistency}} \\ &\quad - \underbrace{D \int_{\partial\Omega} (\nabla v \cdot \mathbf{n})\, (u - u_\text{exact}) \, \mathrm{d}s}_{\text{symmetry}} \\ &\quad + \underbrace{\dfrac{\alpha D}{h} \int_{\partial\Omega} (u - u_\text{exact})\, v \, \mathrm{d}s}_{\text{penalty}} \\ &\quad - \underbrace{\int_\Omega f\, v \, \mathrm{d}x}_{\text{source}} \end{align}

(6)

# Weak Dirichlet (SIPG on boundary)
F += D * (
    -ufl.inner(ufl.grad(u), v * n) * ds
    - ufl.inner(ufl.grad(v), (u - u_ex) * n) * ds
    + (penalty / h) * ufl.inner(u - u_ex, v) * ds
)

# Source
F += -ufl.inner(f, v) * dx

We can then solve the problem and compare against the exct solution with the L2 error

problem = NonlinearProblem(
    F, u, J=ufl.derivative(F, u),
    petsc_options_prefix="advecdiff",
    petsc_options={
        "snes_type": "newtonls",
        "snes_linesearch_type": "none",
        "snes_rtol": 1e-10,
        "snes_atol": 1e-10,
        "snes_max_it": 20,
        "ksp_type": "preonly",
        "pc_type": "lu",
    },
)
u = problem.solve()
u.x.scatter_forward()

error_form = fem.form((u - u_ex) ** 2 * dx)
l2_error = np.sqrt(msh.comm.allreduce(fem.assemble_scalar(error_form), op=MPI.SUM))

print(f"L2 error: {l2_error:.3e}")

L2 error: 1.547e-03

def eval_at_points(uh, pts_2d):
    """Evaluate a DOLFINx Function at an (N, 2) array of (x, y) points."""
    msh = uh.function_space.mesh
    pts_3d = np.column_stack([pts_2d, np.zeros(len(pts_2d))])
    bb = geometry.bb_tree(msh, msh.topology.dim)
    candidates = geometry.compute_collisions_points(bb, pts_3d)
    colliding  = geometry.compute_colliding_cells(msh, candidates, pts_3d)
    cells = np.array(
        [colliding.links(i)[0] if len(colliding.links(i)) > 0 else -1
         for i in range(len(pts_2d))], dtype=np.int32
    )
    vals = np.full(len(pts_2d), np.nan)
    valid = cells >= 0
    if valid.any():
        vals[valid] = uh.eval(pts_3d[valid], cells[valid]).flatten()
    return vals


def compute_arc_length(xs, ys):
    points = np.vstack((xs, ys)).T
    dist   = np.linalg.norm(points[1:] - points[:-1], axis=1)
    return np.insert(np.cumsum(dist), 0, 0.0)


# Colourmap grid
Np = 200
xv, yv = np.linspace(0, 1, Np), np.linspace(0, 1, Np)
Xg, Yg = np.meshgrid(xv, yv)
U_ex_grid = u_exact_np(Xg, Yg, D_value)
U_h_grid  = eval_at_points(u, np.column_stack([Xg.ravel(), Yg.ravel()])).reshape(Np, Np)
vmin, vmax = U_ex_grid.min(), U_ex_grid.max()

# Sample profiles
profiles = [
    {"start": (0.0, 0.0), "end": (1.0, 1.0)},
    {"start": (0.2, 0.8), "end": (0.7, 0.2)},
    {"start": (0.2, 0.6), "end": (0.8, 0.8)},
]

fig, axs = plt.subplots(1, 3, figsize=(15, 5))

for ax, data, title in zip(axs[:2],
                            [U_ex_grid, U_h_grid],
                            ["Exact solution", "Computed solution"]):
    CS = ax.pcolormesh(Xg, Yg, data, cmap="RdBu_r", shading="auto",
                       vmin=vmin, vmax=vmax, rasterized=True)
    CS.set_edgecolor("face")
    plt.colorbar(CS, ax=ax, shrink=0.8)
    ax.set_title(title, fontsize=11)
    ax.set_xlabel("$x$", fontsize=10)
    ax.set_aspect("equal")

axs[0].set_ylabel("$y$", fontsize=10)
axs[0].sharey(axs[1])
plt.setp(axs[1].get_yticklabels(), visible=False)

# Profiles
for profile in profiles:
    sx, sy = profile["start"]
    ex, ey = profile["end"]

    (line,) = axs[1].plot([sx, ex], [sy, ey])
    c = line.get_color()

    # Exact: fewer points, markers
    xs_ex = np.linspace(sx, ex, 30)
    ys_ex = np.linspace(sy, ey, 30)
    arc_ex = compute_arc_length(xs_ex, ys_ex)
    u_ex_vals = u_exact_np(xs_ex, ys_ex, D_value)

    # Computed: more points, solid line
    xs_h = np.linspace(sx, ex, 100)
    ys_h = np.linspace(sy, ey, 100)
    arc_h = compute_arc_length(xs_h, ys_h)
    u_h_vals = eval_at_points(u, np.column_stack([xs_h, ys_h]))

    axs[2].plot(arc_ex, u_ex_vals, marker="o", linestyle="None",
                color=c, alpha=0.3)
    axs[2].plot(arc_h,  u_h_vals,  color=c)

axs[2].set_xlabel("Arc length", fontsize=10)
axs[2].set_ylabel("Solution",   fontsize=10)
axs[2].grid(alpha=0.3)
axs[2].spines[["right", "top"]].set_visible(False)

legend_marker = mpl.lines.Line2D([], [], color="black", marker="o",
                                  linestyle="None", alpha=0.5, label="Exact")
legend_line   = mpl.lines.Line2D([], [], color="black", label="Computed")
axs[2].legend(handles=[legend_marker, legend_line], fontsize=9)

plt.tight_layout()
fig.savefig('../images/mwe1_comparison.svg', bbox_inches='tight', dpi=200)
plt.close(fig)

Convergence Test¶

The mesh is refined uniformly ( $N = 8, 16, \ldots, 128$ ) and the $L^2$ error $\|u_h - u_\text{exact}\|_{L^2(\Omega)}$ is tracked. For a degree- $p$ DG space with upwind advection flux, the expected $L^2$ convergence rate on general unstructured meshes is $\mathcal{O}(h^{p+1/2})$ . This is one half-order below the optimal elliptic rate $\mathcal{O}(h^{p+1})$ obtained by SIPG for pure diffusion, and reflects the hyperbolic character of the upwind flux. For $p = 1$ the expected rate is therefore $\mathcal{O}(h^{1.5})$ , which is used as the reference slope in the plot below.

def solve_advection_diffusion(N, degree=1, D_value=0.1, w_mag=1.0):
    msh = mesh.create_unit_square(
        MPI.COMM_WORLD, N, N,
        cell_type=mesh.CellType.triangle,
        diagonal=mesh.DiagonalType.crossed,
    )

    V = fem.functionspace(msh, ("DG", degree))
    u = fem.Function(V)
    v = ufl.TestFunction(V)

    x = ufl.SpatialCoordinate(msh)
    n = ufl.FacetNormal(msh)
    h = ufl.CellDiameter(msh)

    D  = fem.Constant(msh, PETSc.ScalarType(D_value))
    w  = fem.Constant(msh, PETSc.ScalarType((0, w_mag)))

    u_ex    = exact_solution(x, D)
    f       = -ufl.div(D * ufl.grad(u_ex)) + ufl.dot(w, ufl.grad(u_ex))
    penalty = fem.Constant(msh, PETSc.ScalarType(10.0 * degree**2))
    lmbda   = ufl.conditional(ufl.gt(ufl.dot(w, n), 0), 1, 0)

    dx, ds, dS = ufl.dx, ufl.ds, ufl.dS

    F = 0

    # Advection, upwind flux
    F += -ufl.inner(w * u, ufl.grad(v)) * dx
    F += ufl.inner(2 * ufl.avg(lmbda * w * u), ufl.jump(v, n)) * dS
    F += ufl.inner(lmbda * ufl.dot(w, n) * u, v) * ds

    # Diffusion, SIPG
    F += D * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
    F += -D * ufl.inner(ufl.avg(ufl.grad(u)), ufl.jump(v, n)) * dS
    F += -D * ufl.inner(ufl.jump(u, n), ufl.avg(ufl.grad(v))) * dS
    F += D * (penalty / ufl.avg(h)) * ufl.inner(ufl.jump(u, n), ufl.jump(v, n)) * dS
    
    # Weak Dirichlet (SIPG on boundary)
    F += D * (
        -ufl.inner(ufl.grad(u), v * n) * ds
        - ufl.inner(ufl.grad(v), (u - u_ex) * n) * ds
        + (penalty / h) * ufl.inner(u - u_ex, v) * ds
    )

    # Inflow BC
    F += -ufl.inner((1 - lmbda) * ufl.dot(w, n) * u_ex, v) * ds
    
    # Source
    F += -ufl.inner(f, v) * dx

    problem = NonlinearProblem(
        F, u, J=ufl.derivative(F, u),
        petsc_options_prefix="advecdiff",
        petsc_options={
            "snes_type": "newtonls",
            "snes_linesearch_type": "none",
            "snes_rtol": 1e-10,
            "snes_atol": 1e-10,
            "snes_max_it": 20,
            "ksp_type": "preonly",
            "pc_type": "lu",
        },
    )
    u = problem.solve()
    u.x.scatter_forward()

    error_form = fem.form((u - u_ex) ** 2 * dx)
    l2_error = np.sqrt(msh.comm.allreduce(fem.assemble_scalar(error_form), op=MPI.SUM))

    return l2_error, u


errors = []
Ns = [8, 16, 32, 64, 128]

for N in Ns:
    err, _ = solve_advection_diffusion(degree=1, w_mag=0.1, N=N, D_value=0.1)
    errors.append(err)
    print(f"N = {N:3d},  L2 error = {err:.6e}")

h = np.array([1 / N for N in Ns], dtype=float)
errors = np.array(errors)

fig, ax = plt.subplots(figsize=(5, 4))
ax.loglog(h, errors, "o-")
ax.loglog(h, 2 * h**1.5, "--k")
ax.annotate("Order 1.5", (h[0], 2 * h[0]**1.5),
            textcoords="offset points", xytext=(8, 0), fontsize=9)
ax.set_ylabel("$L^2$ error", fontsize=10)
ax.set_xlabel("Element size ($h$)", fontsize=10)
ax.grid(True, which="major", ls="--", lw=0.5)
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
plt.tight_layout()
fig.savefig('../images/mwe1_conv_mesh.svg', bbox_inches='tight', dpi=200)
plt.close(fig)

N =   8,  L2 error = 8.397482e-03
N =  16,  L2 error = 2.425757e-03
N =  32,  L2 error = 6.980604e-04
N =  64,  L2 error = 2.065179e-04
N = 128,  L2 error = 6.386110e-05

The second test fixes the mesh at $N = 32$ and sweeps the polynomial degree $p = 1, 2, 3, 4$ , with $D = 0.1$ and $|\mathbf{w}| = 0.1$ . At this mesh size the cell Péclet number is $\text{Pe} = |\mathbf{w}|\,h\,/\,(2D) \approx 0.016$ , placing the problem firmly in the diffusion-dominated regime. In this regime the SIPG method recovers the optimal elliptic rate $\mathcal{O}(h^{p+1})$ per degree, so for a fixed mesh the $L^2$ error should decrease monotonically and rapidly with increasing $p$ . The manufactured solution is analytic on $[0,1]^2$ , which favours fast algebraic (and in principle exponential) convergence in $p$ .

errors = []
degrees = [1, 2, 3, 4]
N = 32

for deg in degrees:
    err, _ = solve_advection_diffusion(degree=deg, w_mag=0.1, N=N, D_value=0.1)
    errors.append(err)
    print(f"Deg = {deg:3d},  L2 error = {err:.6e}")

errors = np.array(errors)

fig, ax = plt.subplots(figsize=(5, 4))
ax.plot(degrees, errors, "o-")
ax.set_yscale("log")
ax.set_ylabel("$L^2$ error", fontsize=10)
ax.set_xlabel("Element degree", fontsize=10)
ax.set_xticks(degrees)
ax.grid(True, which="major", ls="--", lw=0.5)
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
plt.tight_layout()
fig.savefig('../images/mwe1_conv_deg.svg', bbox_inches='tight', dpi=200)
plt.close(fig)

Deg =   1,  L2 error = 6.980604e-04
Deg =   2,  L2 error = 6.145345e-05
Deg =   3,  L2 error = 2.194764e-05
Deg =   4,  L2 error = 1.075512e-05