AeroSim Docs
3-D Taylor-Green Vortex (Case 00b)¶
This notebook validates the 3-D TGV case against the pseudo-spectral DNS of van Rees et al. (2011) at \(Re_\Gamma = 1/\nu = 1600\).
The initial condition on a periodic cube of side \(2\pi\) is (Eq. 10, \(\theta = 0\)):
where \(k = 2\pi/N\) in lattice units. Unlike the 2-D case, this flow is nonlinearly unstable and transitions to decaying turbulence. No exact solution exists for \(t > 0\), so validation is based on comparing the energy dissipation rate \(\varepsilon(t) = -dE_k/dt\) against the DNS.
The Reynolds number follows the van Rees convention \(Re = 1/\nu\) (characteristic length \(L^* = 1/k_0 = 1\), velocity \(V_0 = 1\)). In lattice units this corresponds to \(Re = V_0 N / (2\pi\nu)\).
A single \(N = 256\) grid is used with D3Q27/RR-BGK.
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/turbulence/01_taylor_green_vortex_3d/01_taylor_green_vortex_3d.nassu.yaml"
sim_cfgs = ConfigScheme.sim_cfgs_from_file_dct(filename)
sim_cfg = sim_cfgs["taylorGreenVortex3D", 0]
# Protocol constant
V0 = 0.04 # peak velocity amplitude (lattice units)
# Read quantities from config
N = sim_cfg.domain.domain_size.x
TAU = sim_cfg.models.LBM.tau
NU = sim_cfg.models.LBM.kinematic_viscosity
Ma = V0 * np.sqrt(3)
Re = V0 * N / (2.0 * np.pi * NU)
print(f"N={N}: tau={TAU:.6f}, nu={NU:.7f}, Re={Re:.0f}, Ma={Ma:.4f}")
N=256: tau=0.503056, nu=0.0010187, Re=1600, Ma=0.0693
DNS reference data (van Rees et al. 2011)¶
van Rees, W.M., Leonard, A., Pullin, D.I. (2011). A comparison of vortex and pseudo-spectral methods for the simulation of periodic vortical flows at high Reynolds numbers. Journal of Computational Physics, 230, 2794-2805.
The DNS was run at \(Re_\Gamma = 1/\nu = 1600\) on a \([0, 2\pi]^3\) periodic domain using a pseudo-spectral method (PS). The reference data is digitised from Figure 8a and 8b (768\(^3\) PS results, the most converged resolution). Figure 8a provides the overview (\(0 \le t \le 7\)) and Figure 8b the close-up of the peak region (\(7.5 \le t \le 10\)) at higher y-axis resolution.
The dissipation is \(\varepsilon(t) = -dE_k/dt\) in van Rees units (\(V_0 = 1\), \(k_0 = 1\)). To convert lattice quantities: \(t_{vR} = t_L \cdot 2\pi V_0/N\) and \(E_{vR} = E_L / V_0^2\).
[3]:
# van Rees et al. 2011, Re=1600, pseudo-spectral DNS (768^3)
# Digitised from Figure 8a (overview) and 8b (close-up of peak region).
# PS results (solid black curve).
import csv
import pathlib
_csv_path = pathlib.Path(
"validation/turbulence/01_taylor_green_vortex_3d/reference/vanrees2011_dns_eps.csv"
)
with open(_csv_path) as _f:
_reader = csv.DictReader(filter(lambda row: not row.startswith("#"), _f))
_rows = list(_reader)
dns_t = np.array([float(r["t"]) for r in _rows])
dns_eps = np.array([float(r["eps"]) for r in _rows])
print(f"Loaded {len(dns_t)} DNS reference points from {_csv_path.name}")
print(f"Peak dissipation: eps = {dns_eps.max():.5f} at t = {dns_t[np.argmax(dns_eps)]:.2f}")
Loaded 45 DNS reference points from vanrees2011_dns_eps.csv
Peak dissipation: eps = 0.01289 at t = 8.86
HDF5/XDMF output loader¶
Nassu writes snapshots to rolling HDF5 files indexed by an XDMF manifest. The loader below reads kinetic energy directly from those files without assembling full 3-D fields in memory.
[4]:
def load_kinetic_energy(src):
"""Volume-averaged kinetic energy E(t) over every exported step.
E(t) = mean( 0.5 * rho * (ux^2 + uy^2 + uz^2) ) over all cells.
Returns (times, E) as float arrays sorted by time.
"""
times = []
energies = []
for step in src.steps:
arrays, _ = src.read_arrays(step, ["rho", "ux", "uy", "uz"])
ke = 0.5 * arrays["rho"] * (arrays["ux"] ** 2 + arrays["uy"] ** 2 + arrays["uz"] ** 2)
times.append(src.steps_to_time[step])
energies.append(float(np.nanmean(ke)))
idx = np.argsort(times)
return np.array(times)[idx], np.array(energies)[idx]
Load simulation output¶
Kinetic energy is read directly from the HDF5 files and dissipation is estimated by central finite differences on the resulting time series.
[5]:
PROJECT_ROOT = common.find_project_root()
try:
src = common.FieldSource.from_cfg(sim_cfg, project_root=PROJECT_ROOT)
except FileNotFoundError as exc:
src = None
print(f"{exc} - run the simulation first.")
if src is not None:
times_L, E_L = load_kinetic_energy(src)
# Convert to van Rees units
times_vr = times_L * 2.0 * np.pi * V0 / N
E_vr = E_L / V0**2
print(f"N={N} (sim_id={sim_cfg.sim_id:03d}): {len(times_L)} snapshots loaded")
print(f"Lattice time range: {times_L[0]:.0f} - {times_L[-1]:.0f}")
print(f"Van Rees time range: {times_vr[0]:.2f} - {times_vr[-1]:.2f}")
N=256 (sim_id=000): 49 snapshots loaded
Lattice time range: 0 - 12000
Van Rees time range: 0.00 - 11.78
Dissipation rate¶
The dissipation rate is estimated from the energy time series by central finite differences:
[6]:
def compute_dissipation(E: np.ndarray, t: np.ndarray) -> np.ndarray:
"""Estimate -dE/dt by second-order central finite differences.
Args:
E: kinetic energy array, shape (n,).
t: time array, shape (n,).
Returns:
Dissipation rate array, shape (n,).
"""
eps = np.empty_like(E)
# Central differences for interior points
eps[1:-1] = -(E[2:] - E[:-2]) / (t[2:] - t[:-2])
# One-sided at boundaries
eps[0] = -(E[1] - E[0]) / (t[1] - t[0])
eps[-1] = -(E[-1] - E[-2]) / (t[-1] - t[-2])
return eps
if "times_vr" in dir() and len(times_vr) > 0:
# Dissipation in van Rees units: -dE_vr/dt_vr
epsilon_vr = compute_dissipation(E_vr, times_vr)
idx_peak = int(np.argmax(epsilon_vr))
print(f"N={N}: peak dissipation = {epsilon_vr[idx_peak]:.5f} at t = {times_vr[idx_peak]:.2f}")
N=256: peak dissipation = 0.01288 at t = 8.84
Plot dissipation rate - simulation vs DNS¶
The dissipation rate \(\varepsilon(t) = -dE_k/dt\) is plotted in van Rees units, matching Figure 8 of the paper. The peak occurs at \(t \approx 9\) with \(\varepsilon \approx 0.0129\).
[7]:
fig, ax = plt.subplots(figsize=(6, 5))
# DNS reference (van Rees 2011, 768^3 PS)
ax.plot(
dns_t,
dns_eps,
**common.markers.exp_line(linestyle="-.", linewidth=2),
label="PS DNS (van Rees 2011, $768^3$)",
)
if "epsilon_vr" in dir() and len(epsilon_vr) > 0:
ax.plot(
times_vr,
epsilon_vr,
**common.markers.sim_line(linestyle="--", linewidth=2),
label=f"AeroSim N={N} (RR-BGK)",
)
ax.set_xlabel(r"$t$")
ax.set_ylabel(r"$\varepsilon(t) = -dE_k/dt$")
ax.set_title("3-D TGV: energy dissipation rate")
ax.set_xlim(0, 11)
ax.set_ylim(0, 0.015)
ax.legend()
plt.tight_layout()
plt.show()
[8]:
# DNS reference peak (van Rees 2011)
DNS_PEAK_TIME = dns_t[np.argmax(dns_eps)]
DNS_PEAK_EPS = float(dns_eps.max())
if "epsilon_vr" in dir():
idx_peak = int(np.argmax(epsilon_vr))
sim_peak_time = float(times_vr[idx_peak])
sim_peak_eps = float(epsilon_vr[idx_peak])
print(f"Nassu (N={N}, RR-BGK): peak eps = {sim_peak_eps:.5f} at t = {sim_peak_time:.2f}")
print(f"DNS reference (768^3): peak eps = {DNS_PEAK_EPS:.5f} at t = {DNS_PEAK_TIME:.2f}")
print(f"Relative error in peak eps: {abs(sim_peak_eps - DNS_PEAK_EPS) / DNS_PEAK_EPS:.1%}")
print(f"Time shift: {abs(sim_peak_time - DNS_PEAK_TIME):.2f}")
Nassu (N=256, RR-BGK): peak eps = 0.01288 at t = 8.84
DNS reference (768^3): peak eps = 0.01289 at t = 8.86
Relative error in peak eps: 0.1%
Time shift: 0.03
Summary¶
A passing result shows:
The dissipation rate \(\varepsilon(t)\) reproduces the DNS peak location (\(t \approx 9\)) and magnitude (\(\varepsilon \approx 0.013\)) within acceptable tolerance for the \(N = 256\) grid.
\(E_k(t)\) decays smoothly, consistent with the turbulent cascade.
Significant under-prediction of the peak dissipation or a large time shift relative to the DNS indicates insufficient resolution or numerical issues.
Version¶
[9]:
sim_info = sim_cfg.output.read_info()
nassu_commit = sim_info["commit"]
nassu_version = sim_info["version"]
print("Version:", nassu_version)
print("Commit hash:", nassu_commit)
Version: 2.0.0a7
Commit hash: 5e47f2762575c2d285254d36bfd354b6af09fda1
Configuration¶
[10]:
from IPython.display import Code
Code(filename=filename)
[10]:
# Taylor-Green Vortex - 3D validation case
#
# Validates the 3D TGV at N=256 using D3Q27/RRBGK against the
# pseudo-spectral DNS of van Rees et al. (2011) at Re = 1600.
#
# Van Rees definition: Re = 1/nu (with V0=1, k=1 on [0,2*pi]^3).
# Lattice equivalent: Re = V0 * N / (2*pi * nu) = 1600
# V0=0.04, N=256: nu = V0*N/(2*pi*1600) = 0.001019, tau = 3*nu + 0.5 = 0.503056
# Ma = V0*sqrt(3) = 0.069
#
# Time mapping: t_vanRees = t_lattice * 2*pi*V0/N
# Peak dissipation from van Rees 2011 (768^3 PS): t ~ 9.
# In lattice steps: t_peak = 9 * N / (2*pi*V0) = 9168.
# n_steps = 12000 covers t ~ 11.8 in van Rees units.
#
# Reference: van Rees, W.M., Leonard, A., Pullin, D.I. (2011).
# J. Comput. Phys., 230, 2794-2805.
#
# Note: requires models.initialization.equations (PR #443 / feat/init-equation-sem-field).
variables:
V0: 0.04
simulations:
- name: taylorGreenVortex3D
save_path: ./validation/turbulence/01_taylor_green_vortex_3d/results
n_steps: 12000
report:
frequency: 1000
domain:
domain_size:
x: 256
y: 256
z: 256
block_size: 8
data:
exports:
default:
macrs: [rho, u]
interval:
frequency: 250
lvl: 0
target:
volumes:
default: {}
outputs:
instantaneous: true
plane_series:
macrs: [rho, u]
interval: {frequency: 1000, lvl: 0}
target:
planes:
# Mid-domain z slice of the vortex breakdown for field
# visualization; the volumetric dump above stays the
# quantitative source for the dissipation curves.
mid_z:
axis: z
axis_pos: 128
dist: 1
outputs:
instantaneous: true
models:
precision:
default: single
LBM:
tau: 0.503056
vel_set: D3Q27
coll_oper: RRBGK
engine:
name: CUDA
BC:
periodic_dims: [true, true, true]
initialization:
equations:
rho: "1"
ux: "0.04 * sin(2 * pi * x / 256) * cos(2 * pi * y / 256) * cos(2 * pi * z / 256)"
uy: "-0.04 * cos(2 * pi * x / 256) * sin(2 * pi * y / 256) * cos(2 * pi * z / 256)"
uz: "0"