Scalar Advection-Diffusion Strip (Case 12b)

This notebook validates the scalar boundary-condition reconstruction pair

  • {class}ScalarUniformInlet at the west inlet and {class}ScalarRegularizedDirichlet at the east outlet - against the exact steady-state 1D advection-diffusion profile. A uniform fluid \(u = (U, 0, 0)\) carries a passive scalar from a Dirichlet inlet (\(\phi = 1\)) to a Dirichlet outlet (\(\phi = 0\)); \(y\) and \(z\) are periodic so the problem is effectively 1D in \(x\).

The steady-state equation \(U\, d\phi/dx = D\, d^2\phi/dx^2\) with \(\phi(0) = 1\), \(\phi(L) = 0\) has the closed form

\[\phi(x) = \frac{e^{\mathrm{Pe}\,x/L} - e^{\mathrm{Pe}}}{1 - e^{\mathrm{Pe}}}, \qquad \mathrm{Pe} = \frac{U L}{D}.\]

Protocol: the Peclet number is held near \(\mathrm{Pe} = 8\) (mixed advection-diffusion regime with visible curvature) across four grids (\(N = 16, 32, 64, 128\) along \(x\)) by ramping \(D = N / 1600\) at fixed \(U = 0.005\). The two Dirichlet nodes sit exactly on the boundary lattice nodes \(x = 0\) and \(x = N - 1\), so the physical separation is \(L = N - 1\) lattice spacings and \(\mathrm{Pe} = U (N - 1) / D \to 8\) as \(N\) grows. The normalised \(L_2\) error against the analytical profile should decrease as \(O(\Delta x^2)\).

Setup

[1]:

import matplotlib.pyplot as plt
import numpy as np

import nassu.viz as common
from nassu.cfg.model import ConfigScheme

common.use_style()

Load simulation configuration

[2]:

filename = "validation/scalar_transport/01_passive_scalar_transport/01.1_scalar_advection_diffusion_strip.nassu.yaml"
sim_cfgs = ConfigScheme.sim_cfgs_from_file_dct(filename)

# Unrolled simulations: sim_id 0..3 correspond to N = 16, 32, 64, 128.
sim_cfgs_list = [sim_cfgs["advectionDiffusionStrip", i] for i in range(4)]

# Protocol constants read from the config (lattice units).
PHI_INLET = sim_cfgs_list[0].models.scalar_transports[
    "scalar"
].BC.BC_map[0].params["phi_inlet"]
PHI_OUTLET = sim_cfgs_list[0].models.scalar_transports[
    "scalar"
].BC.BC_map[1].params["phi_w"]
# Streamwise fluid velocity from the west UniformFlow inlet.
U_LBM = sim_cfgs_list[0].models.BC.BC_map[0].params["ux"]

# Per-grid quantities. The two Dirichlet nodes lie on x = 0 and x = N - 1,
# so the physical wall-to-wall separation is L = N - 1 lattice spacings and
# Pe = U (N - 1) / D approaches the nominal design value of 8 as N grows.
GRID_SIZES = [c.domain.domain_size.x for c in sim_cfgs_list]
N_STEPS = {c.domain.domain_size.x: c.n_steps for c in sim_cfgs_list}
D_LBM = {
    c.domain.domain_size.x: c.models.scalar_transports[
        "scalar"
    ].adv_diff_equation.D
    for c in sim_cfgs_list
}
PE = {N: U_LBM * (N - 1) / D_LBM[N] for N in GRID_SIZES}

print(f"phi_inlet={PHI_INLET}, phi_outlet={PHI_OUTLET}, U={U_LBM}")
for N in GRID_SIZES:
    print(
        f"N={N:3d}: D={D_LBM[N]:.4f}, L=N-1={N - 1:3d}, "
        f"Pe={PE[N]:.4f}, steps={N_STEPS[N]}"
    )

phi_inlet=1.0, phi_outlet=0.0, U=0.005
N= 16: D=0.0100, L=N-1= 15, Pe=7.5000, steps=12800
N= 32: D=0.0200, L=N-1= 31, Pe=7.7500, steps=25600
N= 64: D=0.0400, L=N-1= 63, Pe=7.8750, steps=51200
N=128: D=0.0800, L=N-1=127, Pe=7.9375, steps=102400

Analytical solution

The exact steady-state profile evaluated on the integer lattice positions \(x = 0, 1, \ldots, N - 1\), with the normalised coordinate \(\xi = x / (N - 1)\) running from the inlet node to the outlet node.

[3]:

def analytical_phi(N):
    """Exact steady advection-diffusion profile on the N lattice nodes.

    phi(xi) = (exp(Pe xi) - exp(Pe)) / (1 - exp(Pe)),  xi = x / (N - 1).
    """
    xi = np.arange(N) / (N - 1)
    Pe = PE[N]
    return (np.exp(Pe * xi) - np.exp(Pe)) / (1.0 - np.exp(Pe))


def l2_error(phi_lbm, phi_exact):
    """Normalised L2 error, E = sqrt(mean((phi_lbm - phi_exact)^2))."""
    return np.sqrt(np.mean((phi_lbm - phi_exact) ** 2))

Load simulation output

The steady scalar field is read from the default volume export at the final step through FieldSource (repo-root-relative paths from the loaded config). The field is uniform in \(y, z\) to within machine precision, so the 1D profile is taken as the spanwise mean \(\overline{\phi}(x)\).

[4]:

PROJECT_ROOT = common.find_project_root()

data = {}
for sim_cfg in sim_cfgs_list:
    N = sim_cfg.domain.domain_size.x
    try:
        src = common.FieldSource.from_cfg(sim_cfg, project_root=PROJECT_ROOT)
    except FileNotFoundError as exc:
        print(f"N={N}: {exc} - skipping.")
        continue

    step_end = src.steps[-1]
    arrays, _ = src.read_arrays(step_end, ["scalar_phi"])
    phi = arrays["scalar_phi"].astype(np.float64)
    # Spanwise-mean profile along x; spread across y, z quantifies 1D-ness.
    phi_x = phi.mean(axis=(1, 2))
    span_spread = phi.std(axis=(1, 2)).max()

    data[N] = {
        "phi_x": phi_x,
        "span_spread": span_spread,
        "t_used": src.steps_to_time[step_end],
    }
    print(
        f"N={N:3d} (sim_id={sim_cfg.sim_id:03d}): profile at "
        f"t={data[N]['t_used']:.0f}, max spanwise std={span_spread:.2e}"
    )

N= 16 (sim_id=000): profile at t=12800, max spanwise std=1.73e-07
N= 32 (sim_id=001): profile at t=25600, max spanwise std=9.81e-08
N= 64 (sim_id=002): profile at t=51200, max spanwise std=7.98e-07
N=128 (sim_id=003): profile at t=102400, max spanwise std=1.99e-07

Steady profile vs analytical

Overlay of the simulated spanwise-mean \(\overline{\phi}(x)\) on the exact 1D advection-diffusion solution for each grid. The near-outlet boundary layer sharpens with resolution as the discrete Peclet number approaches 8.

[5]:

fig, ax = common.fig_single()
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

for i, N in enumerate(GRID_SIZES):
    if N not in data:
        continue
    xi = np.arange(N) / (N - 1)
    c = colors[i % len(colors)]
    ax.plot(
        xi,
        data[N]["phi_x"],
        **common.markers.sim(shape="o", color=c),
        label=f"AeroSim N={N} (Pe={PE[N]:.2f})",
    )

# Analytical reference at the finest grid resolution.
N_fine = max(N for N in GRID_SIZES if N in data)
xi_fine = np.linspace(0.0, 1.0, 400)
Pe_fine = PE[N_fine]
phi_ana = (np.exp(Pe_fine * xi_fine) - np.exp(Pe_fine)) / (1.0 - np.exp(Pe_fine))
ana_label = (
    r"Analytical $\phi=(e^{\mathrm{Pe}\,\xi}-e^{\mathrm{Pe}})"
    r"/(1-e^{\mathrm{Pe}})$" + f" (Pe={Pe_fine:.2f})"
)
ax.plot(
    xi_fine,
    phi_ana,
    **common.markers.exp_line(linestyle='--'),
    label=ana_label,
)

ax.set_xlabel(r"$x / (N-1)$")
ax.set_ylabel(r"$\phi$")
ax.set_title("Advection-diffusion strip: steady profile vs analytical")
ax.legend()
plt.tight_layout()
plt.show()
../../../_images/validation_scalar_transport_01_passive_scalar_transport_01c_scalar_advection_diffusion_strip_10_0.png

\(L_2\) convergence

Normalised \(L_2\) error of the spanwise-mean profile against the exact solution, \(E_{L_2} = \sqrt{\frac{1}{N} \sum_i (\phi^{\mathrm{lbm}}_i - \phi^{\mathrm{exact}}_i)^2}\), evaluated at the final step. A second-order BC reconstruction gives slope \(-2\) on the log-log \(E_{L_2}\)-vs-\(N\) plot.

[6]:

N_values = []
errors = []
max_errors = []
for N in GRID_SIZES:
    if N not in data:
        continue
    phi_exact = analytical_phi(N)
    phi_sim = data[N]["phi_x"]
    err = l2_error(phi_sim, phi_exact)
    N_values.append(N)
    errors.append(err)
    max_errors.append(np.max(np.abs(phi_sim - phi_exact)))
    print(f"N={N:3d}: L2 error = {err:.4e}, max|err| = {max_errors[-1]:.4e}")

N_values = np.array(N_values)
errors = np.array(errors)
slope = np.polyfit(np.log(N_values), np.log(errors), 1)[0]
print(f"fitted convergence slope (E_L2 vs N) = {slope:.2f}")

N= 16: L2 error = 3.5355e-03, max|err| = 1.3119e-02
N= 32: L2 error = 5.9974e-04, max|err| = 1.8713e-03
N= 64: L2 error = 2.2234e-04, max|err| = 3.8862e-04
N=128: L2 error = 3.8747e-05, max|err| = 7.1427e-05
fitted convergence slope (E_L2 vs N) = -2.10

[7]:

fig, ax = common.fig_single()

ax.loglog(
    N_values,
    errors,
    **common.markers.sim(shape="o", linestyle="-"),
    label=f"AeroSim D3Q7 RR-BGK (slope {slope:.2f})",
)

# O(N^-2) reference line (second-order spatial accuracy).
N_ref = np.array([N_values.min() * 0.8, N_values.max() * 1.25])
scale = errors[0] * N_values[0] ** 2
ax.loglog(
    N_ref,
    scale * N_ref ** (-2.0),
    **common.markers.exp_line(linestyle="--"),
    label=r"$O(N^{-2})$",
)

ax.set_xlabel(r"$N$")
ax.set_ylabel(r"$E_{L_2}$")
ax.set_xticks(N_values)
ax.set_xticklabels([str(N) for N in N_values])
ax.set_title("Advection-diffusion strip: spatial convergence (Pe ~ 8)")
ax.legend()
plt.tight_layout()
plt.show()
../../../_images/validation_scalar_transport_01_passive_scalar_transport_01c_scalar_advection_diffusion_strip_13_0.png

Summary

A passing result shows:

  1. The simulated spanwise-mean profile overlays the exact advection-diffusion solution, with the inlet and outlet Dirichlet values captured exactly.

  2. The field is 1D to machine precision (negligible spanwise spread).

  3. The \(L_2\) error decreases with grid refinement at close to second order, confirming the {class}ScalarUniformInlet / {class}ScalarRegularizedDirichlet reconstruction pair is second-order accurate.

Version

[8]:

sim_info = sim_cfgs_list[0].output.read_info()

print("Version:", sim_info["version"])
print("Commit hash:", sim_info["commit"])

Version: 2.0.0a7
Commit hash: 5e47f2762575c2d285254d36bfd354b6af09fda1

Configuration

[9]:

from IPython.display import Code

Code(filename=filename)

[9]:

# Scalar advection-diffusion strip - ScalarUniformInlet validation
#
# Variant of case 12 (passive scalar transport) that exercises the
# scalar `ScalarUniformInlet` BC at the west face combined with a
# `ScalarRegularizedDirichlet` outlet at the east face. The fluid is
# a uniform u = (U, 0, 0) carried by `UniformFlow` / Neumann-outlet
# fluid BCs; y and z are periodic so the problem is effectively 1D
# in x.
#
# Analytical steady-state solution of d/dx(U phi) = D d^2 phi/dx^2
# with phi(0) = 1 and phi(L) = 0:
#
#     phi(x) = (exp(Pe x / L) - exp(Pe)) / (1 - exp(Pe))
#     Pe = U L / D
#
# Four grid resolutions (N = 16, 32, 64, 128) hold the Peclet number
# constant at Pe = 8 (mixed advection-diffusion regime with visible
# curvature in the steady profile):
#   - U_lbm = 0.005   (constant; Mach safely below 0.1)
#   - D_lbm = N / 1600 -> 0.01, 0.02, 0.04, 0.08
#     (tau_phi - 1/2 = 4 D for D3Q7 cs2_phi = 1/4 -> tau_phi ramps
#      0.54 to 0.82 across grids)
#   - n_steps ~ 4 L / U so the advective transient clears

simulations:
  - name: advectionDiffusionStrip
    save_path: ./validation/scalar_transport/01_passive_scalar_transport/results/advection_diffusion_strip

    n_steps: !unroll [12800, 25600, 51200, 102400]

    report:
      frequency: 1000

    domain:
      domain_size:
        x: !unroll [16, 32, 64, 128]
        y: !unroll [8, 8, 8, 8]
        z: !unroll [8, 8, 8, 8]
      block_size: 8

    data:
      exports:
        default:
          macrs: [rho, u, scalar_phi, scalar_q_neq]
          interval:
            frequency: !unroll [640, 1280, 2560, 5120]
            lvl: 0
          target:
            volumes:
              default: {}
          outputs:
            instantaneous: true
        plane_series:
          macrs: [scalar_phi]
          interval: {frequency: !unroll [640, 1280, 2560, 5120], lvl: 0}
          target:
            planes:
              # Streamwise plane through the advected scalar front.
              mid_span:
                axis: y
                axis_pos: 4
                dist: 1
          outputs:
            instantaneous: true

    models:
      precision:
        default: single

      LBM:
        tau: 0.6
        vel_set: D3Q27
        coll_oper: RRBGK

      engine:
        name: CUDA

      # West and east drive inlet / outlet; y and z are periodic so
      # the problem stays 1D in x.
      BC:
        periodic_dims: [false, true, true]
        BC_map:
          - pos: W
            BC: UniformFlow
            params:
              rho: 1.0
              ux: 0.005
              uy: 0.0
              uz: 0.0
          - pos: E
            BC: RegularizedNeumannOutlet
            wall_normal: E
            params:
              rho: 1.0
      scalar_transports:
        scalar:
          velocity_set: D3Q7
          collision_operator: RRBGK
          adv_diff_equation:
            D: !unroll [0.01, 0.02, 0.04, 0.08]
            S: "0"
          # Initial field is zero; the inlet drives the scalar into the
          # domain over time until the steady-state profile is set.
          initial_field: "0"
          BC:
            BC_map:
              - pos: W
                BC: ScalarUniformInlet
                params:
                  phi_inlet: 1.0
              - pos: E
                BC: ScalarRegularizedDirichlet
                wall_normal: E
                params:
                  phi_w: 0.0