Scalar Wall Dirichlet Conduction (Case 12c)

This notebook validates the {class}ScalarRegularizedDirichlet BC by driving pure 1D conduction between two opposing Dirichlet walls with the fluid at rest. The west wall holds \(\phi_w = 0\) and the east wall \(\phi_w = 1\); \(y\) and \(z\) are periodic, so the steady state reduces to the textbook linear ramp.

The steady pure-diffusion equation \(D\, d^2\phi/dx^2 = 0\) with \(\phi(0) = 0\), \(\phi(L) = 1\) is the linear profile

\[\phi(x) = \frac{x}{L}.\]

Because a linear field has zero second derivative, a second-order scheme reproduces the steady solution exactly (no spatial-truncation error), so the headline check is that the simulated profile is linear to \(R^2 \approx 1\) and matches \(\phi = x/L\). The residual \(L_2\) error at the snapshot is set by the incomplete conduction transient (time constant \(\sim L^2 / (\pi^2 D)\)), not by spatial discretisation.

Protocol: constant diffusivity \(D = 0.05\) (\(\tau_\phi = 0.7\)) across three grids (\(N = 16, 32, 64\) along \(x\)), with \(n_{\text{steps}} \propto N^2\) so each grid reaches the same fraction of steady state. The two Dirichlet nodes sit on the boundary lattice nodes \(x = 0\) and \(x = N - 1\), so \(L = N - 1\).

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.2_scalar_wall_dirichlet_conduction.nassu.yaml"
sim_cfgs = ConfigScheme.sim_cfgs_from_file_dct(filename)

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

# Dirichlet wall values read from the scalar BC map (west, east).
scalar_bc = sim_cfgs_list[0].models.scalar_transports["scalar"].BC.BC_map
PHI_W = scalar_bc[0].params["phi_w"]   # west wall
PHI_E = scalar_bc[1].params["phi_w"]   # east wall

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
}

print(f"phi_west={PHI_W}, phi_east={PHI_E}")
for N in GRID_SIZES:
    D = D_LBM[N]
    tau_transient = (N - 1) ** 2 / (np.pi**2 * D)
    print(
        f"N={N:3d}: D={D:.3f}, steps={N_STEPS[N]}, "
        f"transient tau~L^2/(pi^2 D)={tau_transient:.0f}"
    )

phi_west=0.0, phi_east=1.0
N= 16: D=0.050, steps=4096, transient tau~L^2/(pi^2 D)=456
N= 32: D=0.050, steps=16384, transient tau~L^2/(pi^2 D)=1947
N= 64: D=0.050, steps=65536, transient tau~L^2/(pi^2 D)=8043

Analytical solution

The exact steady profile on the integer lattice nodes is the linear ramp \(\phi(x_i) = i / (N - 1)\) from the west wall (\(\phi = 0\)) to the east wall (\(\phi = 1\)).

[3]:

def analytical_phi(N):
    """Exact linear conduction profile, phi(x_i) = i / (N - 1)."""
    return np.arange(N) / (N - 1)


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. The field is uniform in \(y, z\), 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)
    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=4096, max spanwise std=5.55e-08
N= 32 (sim_id=001): profile at t=16384, max spanwise std=1.13e-07
N= 64 (sim_id=002): profile at t=65536, max spanwise std=3.94e-08

Steady profile vs analytical linear ramp

Overlay of the simulated spanwise-mean \(\overline{\phi}(x)\) on the analytical linear conduction profile \(\phi = x / L\), annotated with the coefficient of determination \(R^2\) of the simulated profile against the exact ramp.

[5]:

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

r2 = {}
for i, N in enumerate(GRID_SIZES):
    if N not in data:
        continue
    xi = np.arange(N) / (N - 1)
    phi_sim = data[N]["phi_x"]
    ss_res = np.sum((phi_sim - xi) ** 2)
    ss_tot = np.sum((phi_sim - phi_sim.mean()) ** 2)
    r2[N] = 1.0 - ss_res / ss_tot
    c = colors[i % len(colors)]
    ax.plot(
        xi,
        phi_sim,
        **common.markers.sim(shape="o", color=c),
        label=f"AeroSim N={N} ($R^2$={r2[N]:.6f})",
    )

xi_line = np.linspace(0.0, 1.0, 200)
ax.plot(
    xi_line,
    xi_line,
    **common.markers.exp_line(linestyle="--"),
    label=r"Analytical $\phi = x / L$",
)

ax.set_xlabel(r"$x / (N-1)$")
ax.set_ylabel(r"$\phi$")
ax.set_title("Wall Dirichlet conduction: steady profile vs analytical")
ax.legend()
plt.tight_layout()
plt.show()
../../../_images/validation_scalar_transport_01_passive_scalar_transport_01d_scalar_wall_dirichlet_conduction_10_0.png

\(L_2\) error vs grid

Normalised \(L_2\) error of the spanwise-mean profile against the exact linear ramp, \(E_{L_2} = \sqrt{\frac{1}{N} \sum_i (\phi^{\mathrm{lbm}}_i - i/(N-1))^2}\), at the final step. Unlike a curved solution, a linear steady field carries no spatial-truncation error for the second-order scheme, so the residual here is dominated by the incomplete conduction transient rather than \(\Delta x\); the \(O(N^{-1})\) and \(O(N^{-2})\) guides are shown for reference. The residual stays below \(3 \times 10^{-4}\) on every grid while \(R^2 \approx 1\).

[6]:

N_values = []
errors = []
mid_values = []
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)
    mid_values.append(phi_sim[N // 2])
    print(
        f"N={N:3d}: L2 error = {err:.4e}, R^2 = {r2[N]:.8f}, "
        f"phi(mid) = {phi_sim[N // 2]:.5f} (exact {phi_exact[N // 2]:.5f})"
    )

N_values = np.array(N_values)
errors = np.array(errors)

N= 16: L2 error = 5.5849e-05, R^2 = 0.99999997, phi(mid) = 0.53325 (exact 0.53333)
N= 32: L2 error = 1.0713e-04, R^2 = 0.99999987, phi(mid) = 0.51598 (exact 0.51613)
N= 64: L2 error = 2.3597e-04, R^2 = 0.99999935, phi(mid) = 0.50762 (exact 0.50794)

[7]:

fig, ax = common.fig_single()

ax.loglog(
    N_values,
    errors,
    **common.markers.sim(shape="o", linestyle="-"),
    label="AeroSim D3Q7 RR-BGK (transient residual)",
)

N_ref = np.array([N_values.min() * 0.8, N_values.max() * 1.25])
scale1 = errors[0] * N_values[0]
ax.loglog(
    N_ref,
    scale1 * N_ref ** (-1.0),
    **common.markers.exp_line(linestyle=":"),
    label=r"$O(N^{-1})$",
)
scale2 = errors[0] * N_values[0] ** 2
ax.loglog(
    N_ref,
    scale2 * 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("Wall Dirichlet conduction: L2 residual vs grid")
ax.legend()
plt.tight_layout()
plt.show()
../../../_images/validation_scalar_transport_01_passive_scalar_transport_01d_scalar_wall_dirichlet_conduction_13_0.png

Summary

A passing result shows:

  1. The simulated profile is a straight line from \(\phi = 0\) at the west wall to \(\phi = 1\) at the east wall, matching \(\phi = x / L\) with \(R^2 \approx 1\).

  2. The mid-domain value is \(\phi(L/2) \approx 0.5\) to within the BC truncation.

  3. The \(L_2\) residual stays at the \(10^{-4}\) level, set by the incomplete conduction transient (a linear steady field has zero spatial-truncation error), confirming the {class}ScalarRegularizedDirichlet walls impose the Dirichlet values exactly.

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 wall Dirichlet conduction - ScalarRegularizedDirichlet validation
#
# Variant of case 12 (passive scalar transport) that exercises the
# scalar `ScalarRegularizedDirichlet` BC on two opposing walls with
# the fluid at rest. y and z are periodic; west / east walls are
# `RegularizedHWBB` (no-slip on the fluid) and Dirichlet on the scalar.
#
# Analytical steady-state solution of d^2 phi/dx^2 = 0 with
# phi(0) = 0 and phi(L) = 1:
#
#     phi(x) = x / L
#
# Three grid resolutions (N = 16, 32, 64) on a fixed scalar
# diffusivity D_lbm = 0.05 (tau_phi = 0.7, well-stable). Time to
# reach 99.9 % of steady state is t ~ 7 L^2 / (pi^2 D); n_steps is
# sized to ~ 8 L^2 / D so the profile is firmly settled before the
# validation snapshot.

simulations:
  - name: wallDirichletConduction
    save_path: ./validation/scalar_transport/01_passive_scalar_transport/results/wall_dirichlet_conduction

    n_steps: !unroll [4096, 16384, 65536]

    report:
      frequency: 1000

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

    data:
      exports:
        default:
          macrs: [rho, u, scalar_phi, scalar_q_neq]
          interval:
            frequency: !unroll [256, 1024, 4096]
            lvl: 0
          target:
            volumes:
              default: {}
          outputs:
            instantaneous: true
        plane_series:
          macrs: [scalar_phi]
          interval: {frequency: !unroll [256, 1024, 4096], lvl: 0}
          target:
            planes:
              # Mid-span plane of the wall-to-wall conduction ramp.
              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

      # Fluid at rest. West and east walls are no-slip (RegularizedHWBB);
      # y / z are periodic.
      BC:
        periodic_dims: [false, true, true]
        BC_map:
          - pos: W
            BC: RegularizedHWBB
            wall_normal: W
          - pos: E
            BC: RegularizedHWBB
            wall_normal: E

      scalar_transports:
        scalar:
          velocity_set: D3Q7
          collision_operator: RRBGK
          adv_diff_equation:
            D: 0.05
            S: "0"
          initial_field: "0"
          BC:
            BC_map:
              - pos: W
                BC: ScalarRegularizedDirichlet
                wall_normal: W
                params:
                  phi_w: 0.0
              - pos: E
                BC: ScalarRegularizedDirichlet
                wall_normal: E
                params:
                  phi_w: 1.0