128 lines
4.1 KiB
Python
128 lines
4.1 KiB
Python
"""DFG Re100 periodic reference comparison with declared pressure-NEQ outlet, not exact do-nothing.
|
|
|
|
CPU helpers import without CUDA. GPU execution uses unchanged d79e3f6 production
|
|
sources. Default pilot is deliberately too short to certify periodic convergence.
|
|
"""
|
|
|
|
import argparse
|
|
import csv
|
|
import hashlib
|
|
import json
|
|
import math
|
|
import os
|
|
from pathlib import Path
|
|
import time
|
|
from tempfile import TemporaryDirectory
|
|
import numpy as np
|
|
|
|
ROOT = Path(__file__).resolve().parents[2]
|
|
PIN = "d79e3f6f0a424f5bcc3a1a6d610313e49daf97a2"
|
|
|
|
|
|
def geometry(n):
|
|
if isinstance(n, bool) or n < 10 or n % 10:
|
|
raise ValueError("Diameter must be a positive multiple of 10, at least 10")
|
|
return dict(
|
|
diameter=n,
|
|
nx=22 * n + 1,
|
|
ny=41 * n // 10 + 2,
|
|
center=[2 * n, 2 * n + 0.5],
|
|
radius=n / 2,
|
|
dx=0.1 / n,
|
|
bottom=0.5,
|
|
top=41 * n / 10 + 0.5,
|
|
inlet_source_x=0.0,
|
|
outlet_source_x=22.0 * n,
|
|
)
|
|
|
|
|
|
def parameters(n, speed):
|
|
g = geometry(n)
|
|
if not 0 < speed <= 0.04:
|
|
raise ValueError("Require 0 < lattice mean speed <= .04")
|
|
return dict(
|
|
**g,
|
|
mean_speed=speed,
|
|
peak_speed=1.5 * speed,
|
|
viscosity=speed * n / 100,
|
|
dt=speed * g["dx"],
|
|
peak_mach=1.5 * speed * math.sqrt(3),
|
|
re=100.0,
|
|
)
|
|
|
|
|
|
def make_config(n, speed, collision, outlet, inlet="regularized"):
|
|
p = parameters(n, speed)
|
|
cfg = json.loads((ROOT / "src/CelerisLab/configs/config_lbm.json").read_text())
|
|
cfg["grid"].update(nx=p["nx"], ny=p["ny"], nz=1, lattice_model="D2Q9")
|
|
cfg["physics"].update(velocity=speed, viscosity=p["viscosity"], rho=1.0)
|
|
cfg["method"].update(
|
|
collision=collision,
|
|
streaming="double_buffer",
|
|
store_precision="FP32",
|
|
ddf_shifting=True,
|
|
y_wall_bc="bounce_back",
|
|
)
|
|
cfg["method"]["les"]["enabled"] = False
|
|
cfg["method"]["inlet"].update(profile="parabolic", scheme=inlet)
|
|
cfg["method"]["outlet"].update(mode=outlet, backflow_clamp=False, blend_alpha=0.7)
|
|
body = {
|
|
"objects": [
|
|
dict(type="cylinder", center=p["center"], radius=p["radius"], omega=0.0)
|
|
]
|
|
}
|
|
return cfg, body, p
|
|
|
|
|
|
def bilinear_fluid(a, fluid, x, y):
|
|
ix, iy = math.floor(x), math.floor(y)
|
|
tx, ty = x - ix, y - iy
|
|
result = 0.0
|
|
for xx, yy, w in (
|
|
(ix, iy, (1 - tx) * (1 - ty)),
|
|
(ix + 1, iy, tx * (1 - ty)),
|
|
(ix, iy + 1, (1 - tx) * ty),
|
|
(ix + 1, iy + 1, tx * ty),
|
|
):
|
|
if w == 0:
|
|
continue
|
|
if not (0 <= yy < a.shape[0] and 0 <= xx < a.shape[1]) or not fluid[yy, xx]:
|
|
raise ValueError("Pressure interpolation touched a non-fluid/outside node")
|
|
result += w * float(a[yy, xx])
|
|
return result
|
|
|
|
|
|
def pressure_probes(rho, fluid, p):
|
|
# Linear extrapolation from distances h and 2h outside the cylinder to
|
|
# each exact wall probe, after bilinear interpolation using only fluid.
|
|
cx, cy = p["center"]
|
|
r = p["radius"]
|
|
values = []
|
|
for sign in (-1, 1):
|
|
wall = cx + sign * r
|
|
near = bilinear_fluid(rho, fluid, wall + sign, cy)
|
|
far = bilinear_fluid(rho, fluid, wall + 2 * sign, cy)
|
|
values.append(2 * near - far)
|
|
scale = (1.0 / p["mean_speed"]) ** 2 / 3
|
|
return dict(
|
|
rho_front=values[0], rho_back=values[1], dp=(values[0] - values[1]) * scale
|
|
)
|
|
|
|
|
|
def outlet_diagnostic(rho, ux, uy, p):
|
|
# Extrapolate fluid values to the declared x=NX-1 source plane. Diagnostic
|
|
# excludes two rows beside each wall; not an exact discrete BC identity.
|
|
rows = slice(3, -3)
|
|
i = -2
|
|
dudx = (5 * ux[rows, i] - 8 * ux[rows, i - 1] + 3 * ux[rows, i - 2]) / 2
|
|
dvdx = (5 * uy[rows, i] - 8 * uy[rows, i - 1] + 3 * uy[rows, i - 2]) / 2
|
|
rho_plane = 2 * rho[rows, i] - rho[rows, i - 1]
|
|
pressure = (rho_plane - 1) / 3
|
|
scale = p["mean_speed"] ** 2
|
|
residual = np.stack((p["viscosity"] * dudx - pressure, p["viscosity"] * dvdx))
|
|
return dict(
|
|
rms_over_rho_U2=float(np.sqrt(np.mean(residual**2)) / scale),
|
|
max_over_rho_U2=float(np.max(abs(residual)) / scale),
|
|
definition="fluid-side one-sided estimate of nu*d_x(u)-p*n at declared outlet",
|
|
)
|