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 numpy as np
import matplotlib.pyplot as plt
from nassu.cfg.model import ConfigScheme
from tests.validation.notebooks import common
common.use_style()
Load simulation configuration¶
[2]:
filename = "tests/validation/cases/00b_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 pathlib
import csv
_csv_path = pathlib.Path("tests/validation/comparison/Taylor_Green/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]:
import h5py
from nassu.cfg.schemes.simul import MacrExport
def load_h5_kinetic_energy(macr_export: MacrExport, dim=3):
"""Compute volume-averaged kinetic energy E(t) from XDMF/HDF5 output.
Reads all rolling HDF5 files under the macrs/instantaneous directory.
Returns (times, E) as float arrays sorted by time.
"""
macrs_dir = macr_export.xdmf_filename.parent
export_name = macr_export.base_path.name
h5_files = sorted(macrs_dir.glob(f"{export_name}.*.h5"))
times = []
energies = []
for h5_path in h5_files:
with h5py.File(h5_path, "r") as f:
for ts_key in sorted(f.keys()):
t = float(ts_key[1:]) # strip leading "t"
total_E = 0.0
total_cells = 0
for block_key in f[ts_key].keys():
rho = f[ts_key][block_key]["rho"][()]
ux = f[ts_key][block_key]["ux"][()]
uy = f[ts_key][block_key]["uy"][()]
uz = f[ts_key][block_key]["uz"][()]
ke = 0.5 * rho * (ux**2 + uy**2 + uz**2)
total_E += float(ke.sum())
total_cells += ke.size
times.append(t)
energies.append(total_E / total_cells)
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]:
macr_export = sim_cfg.output.instantaneous["default"]
if not macr_export.xdmf_filename.exists():
print(f"XDMF output not found ({macr_export.xdmf_filename}) - run the simulation first.")
else:
times_L, E_L = load_h5_kinetic_energy(macr_export, dim=3)
# 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}" f" 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: 1.6.45
Commit hash: ef533df4387f6c78c3790b6168e6170bf0b0b2ca
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: ./tests/validation/results/00b_taylor_green_vortex_3d
n_steps: 12000
report:
frequency: 1000
domain:
domain_size:
x: 256
y: 256
z: 256
block_size: 8
data:
divergence: { frequency: 500 }
instantaneous:
default:
# ~48 snapshots for smooth dissipation curves
interval: { frequency: 250 }
macrs: [rho, u]
statistics:
interval: { frequency: 0 }
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"