Frank14f/CelerisLab
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Preserve bounded validation evidence and historical failures

Browse Source
This commit is contained in:
Frank14f committed 2026-10-15 10:13:51 +08:00
1 parent cdbb9c42f6
commit 263f3b0c7a
239 files changed
+74191
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

research/j1-2026-10-12/incremental.diff
+629
View File
@@ -0,0 +1,629 @@
--- a/src/CelerisLab/common/checkpoint.py
+++ b/src/CelerisLab/common/checkpoint.py
@@ -214,6 +214,10 @@
# Config compatibility check
saved_cfg = json.loads(hf.attrs["config_json"])
+ if "on_site_velocity" in (saved_cfg.get("y_wall_bc"), lbm_cfg.y_wall_bc):
+ for key in ("y_wall_bc", "velocity", "wall_normal_velocity"):
+ if saved_cfg.get(key) != getattr(lbm_cfg, key):
+ raise ValueError(f"Checkpoint on-site boundary mismatch on '{key}'")
for key in ("lattice_model", "nx", "ny", "nz", "store_precision",
"streaming", "ddf_shifting"):
saved_val = saved_cfg.get(key)
@@ -279,6 +283,9 @@
if key in raw_params}
if not all(np.isfinite(value) for value in params.values()):
raise ValueError("Checkpoint runtime parameters must be finite")
+ if lbm_cfg.y_wall_bc == "on_site_velocity" and (
+ bodies.count or any(params.get(k, 0.0) != 0 for k in ("fx", "fy", "fz"))):
+ raise ValueError("on_site_velocity does not support body forces or objects")
obs_accum_steps = int(hf.attrs.get("obs_accum_steps", 0))
if obs_accum_steps < 0:
raise ValueError("Checkpoint obs_accum_steps must be nonnegative")
--- a/src/CelerisLab/config.py
+++ b/src/CelerisLab/config.py
@@ -34,7 +34,7 @@
"regularized": 3,
}
OUTLET_MAP = {"neq_extrap": 0, "zero_gradient": 1, "blended": 2}
-Y_WALL_BC_MAP = {"bounce_back": 0, "free_slip": 1}
+Y_WALL_BC_MAP = {"bounce_back": 0, "free_slip": 1, "on_site_velocity": 2}
DTYPE_MAP = {"FP32": "float", "FP64": "double"}
@@ -86,6 +86,7 @@
inlet_regularized_neq_damp: float = 0.5
outlet_srt_neq_damp: float = 0.5
y_wall_bc: str = "bounce_back"
+ wall_normal_velocity: float = 0.0
omega_min: float = 0.01
omega_max: float = 1.99
@@ -127,7 +128,7 @@
for name in ("viscosity", "velocity", "rho", "les_cs", "trt_magic_param",
"outlet_blend_alpha", "inlet_trt_neq_damp",
"inlet_regularized_neq_damp", "outlet_srt_neq_damp",
- "omega_min", "omega_max"):
+ "omega_min", "omega_max", "wall_normal_velocity"):
if not math.isfinite(getattr(self, name)):
raise ValueError(f"{name} must be finite.")
if self.viscosity <= 0 or self.rho <= 0:
@@ -150,6 +151,13 @@
)
if self.dim == 2 and self.nz != 1:
raise ValueError("nz must be 1 for 2D")
+ if self.y_wall_bc == "on_site_velocity":
+ if (not self.is_d2q9 or self.collision != "SRT"
+ or self.streaming != "double_buffer" or self.store_precision != "FP32"
+ or self.ddf_shifting or self.les_enabled):
+ raise ValueError("on_site_velocity requires D2Q9/SRT/FP32/unshifted/double_buffer without LES")
+ if abs(self.wall_normal_velocity) >= 1:
+ raise ValueError("on_site_velocity requires abs(wall_normal_velocity) < 1")
if not (0.0 <= self.inlet_trt_neq_damp <= 1.0):
raise ValueError("method.inlet.trt_neq_damp must lie in [0, 1].")
if not (0.0 <= self.inlet_regularized_neq_damp <= 1.0):
@@ -187,7 +195,9 @@
"""Parameters that cannot change while field allocations/state are retained."""
self.validate()
return (self.lattice_model, self.nx, self.ny, self.nz, self.data_type,
- self.store_precision, bool(self.ddf_shifting), self.streaming)
+ self.store_precision, bool(self.ddf_shifting), self.streaming,
+ (self.velocity, self.wall_normal_velocity)
+ if self.y_wall_bc == "on_site_velocity" else None)
def to_macros(self) -> Dict[str, Any]:
"""Return flat dict of macro_name → value for header generation."""
@@ -216,6 +226,7 @@
"OUTLET_BLEND_ALPHA": _float_literal(self.outlet_blend_alpha),
"OUTLET_BACKFLOW_CLAMP": int(self.outlet_backflow_clamp),
"Y_WALL_BC": Y_WALL_BC_MAP[self.y_wall_bc],
+ "WALL_NORMAL_VELOCITY": _float_literal(self.wall_normal_velocity),
"OMEGA_COLLISION_MIN": _float_literal(self.omega_min),
"OMEGA_COLLISION_MAX": _float_literal(self.omega_max),
"TRT_MAGIC_PARAM": _float_literal(self.trt_magic_param),
@@ -284,6 +295,7 @@
inlet_regularized_neq_damp=m["inlet"].get("regularized_neq_damp", 0.5),
outlet_srt_neq_damp=m["outlet"].get("srt_neq_damp", 0.5),
y_wall_bc=m.get("y_wall_bc", "bounce_back"),
+ wall_normal_velocity=m.get("wall_normal_velocity", 0.0),
omega_min=m["omega_guard"]["min"],
omega_max=m["omega_guard"]["max"],
threads_per_block=c["threads_per_block"],
--- a/src/CelerisLab/configs/CONFIG.md
+++ b/src/CelerisLab/configs/CONFIG.md
@@ -101,3 +101,18 @@
D2Q9/D3Q19 的 flag 扩展位 `0x0800` 定义为 `FRC_REGION`，用于标记 force_region 类型物体的格子。
- ForceRegionKernel 不依赖 flag 位执行，flag 只用于后处理可视化。
- `merge_flag_masks` 会正确合并 force_region 的 `FLUID | FRC_REGION` 标记。
+# J1 on-site velocity planes (limited reference case)
+
+Set `method.y_wall_bc` to `"on_site_velocity"` and
+`method.wall_normal_velocity` to V (default 0). This selects an x-periodic
+domain with physical planes through fluid rows y=1 and y=NY−2, H=NY−3.
+The bottom velocity is (0,V); the top is (`physics.velocity`,V).
+These are constant **Cartesian** velocities: positive V injects through the
+bottom and leaves through the top. The selected rows do not also receive
+halfway reflection. Initialization uses a linear tangential profile and V.
+
+Supported only for D2Q9/SRT/FP32/unshifted/double_buffer without LES, objects
+or body forces. Boundary targets and the mode are fixed at construction;
+checkpoint restart requires matching targets. This is an entire-plane
+Couette-throughflow reference, not a finite slot or curved jet API. The default
+`bounce_back` mode and existing configurations retain their previous path.
--- a/src/CelerisLab/cuda/compiler_v2.py
+++ b/src/CelerisLab/cuda/compiler_v2.py
@@ -65,6 +65,8 @@
macros = cfg.to_macros()
if isinstance(n_objects, bool) or not isinstance(n_objects, int) or n_objects < 0:
raise ValueError("n_objects must be a nonnegative integer")
+ if cfg.y_wall_bc == "on_site_velocity" and n_objects:
+ raise ValueError("on_site_velocity currently requires an object-free domain")
build_dir = tempfile.mkdtemp(prefix="celeris-build-", dir=build_root)
try:
shutil.copytree(_KERNEL_DIR, build_dir, dirs_exist_ok=True,
@@ -173,6 +175,7 @@
#define OUTLET_BLEND_ALPHA {m['OUTLET_BLEND_ALPHA']}
#define OUTLET_BACKFLOW_CLAMP {m['OUTLET_BACKFLOW_CLAMP']}
#define Y_WALL_BC {m['Y_WALL_BC']}
+#define WALL_NORMAL_VELOCITY {m['WALL_NORMAL_VELOCITY']}
#define OMEGA_COLLISION_MIN {m['OMEGA_COLLISION_MIN']}
#define OMEGA_COLLISION_MAX {m['OMEGA_COLLISION_MAX']}
--- a/src/CelerisLab/lbm/field.py
+++ b/src/CelerisLab/lbm/field.py
@@ -123,6 +123,8 @@
if self.cfg.is_d2q9:
flag[0:self.nx] = SOLID | BC_WALL
flag[(self.ny - 1) * self.nx:self.ny * self.nx] = SOLID | BC_WALL
+ if self.cfg.y_wall_bc == "on_site_velocity":
+ return flag # x periodic; physical y planes are the adjacent fluid rows.
flag[::self.nx] = SOLID | BC_INLET
flag[self.nx - 1::self.nx] = SOLID | BC_OUTLET
flag[0] = SOLID | BC_WALL
@@ -231,6 +233,10 @@
Accepted keys: omega, fx, fy, fz, n_objects.
"""
+ if self.cfg.y_wall_bc == "on_site_velocity" and (
+ any(float(kwargs.get(k, 0.0)) != 0 for k in ("fx", "fy", "fz"))
+ or int(kwargs.get("n_objects", 0)) != 0):
+ raise ValueError("on_site_velocity does not support body forces or objects")
self._runtime_params.update({
"omega": float(kwargs.get("omega", self._runtime_params["omega"])),
"fx": float(kwargs.get("fx", self._runtime_params["fx"])),
--- /dev/null
+++ b/src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh
@@ -0,0 +1,22 @@
+// On-site D2Q9 Zou-He velocity closure, paired lattice ordering.
+// Input is poststream f at the arrival time. Only unknown populations change;
+// the caller owns collision/store. No force, shifting or off-lattice walls.
+#ifndef CELERIS_BOUNDARY_PLANAR_VELOCITY_CUH
+#define CELERIS_BOUNDARY_PLANAR_VELOCITY_CUH
+
+__device__ __forceinline__ void reconstruct_bottom_velocity(float* f, float u, float v)
+{
+ const float rho = (f[0] + f[1] + f[2] + 2.0f*(f[4] + f[6] + f[7])) / (1.0f-v);
+ f[3] = f[4] + (2.0f/3.0f)*rho*v;
+ f[5] = f[6] + 0.5f*(f[2]-f[1]) + 0.5f*rho*u + (1.0f/6.0f)*rho*v;
+ f[8] = f[7] + 0.5f*(f[1]-f[2]) - 0.5f*rho*u + (1.0f/6.0f)*rho*v;
+}
+
+__device__ __forceinline__ void reconstruct_top_velocity(float* f, float u, float v)
+{
+ const float rho = (f[0] + f[1] + f[2] + 2.0f*(f[3] + f[5] + f[8])) / (1.0f+v);
+ f[4] = f[3] - (2.0f/3.0f)*rho*v;
+ f[6] = f[5] + 0.5f*(f[1]-f[2]) - 0.5f*rho*u - (1.0f/6.0f)*rho*v;
+ f[7] = f[8] + 0.5f*(f[2]-f[1]) + 0.5f*rho*u - (1.0f/6.0f)*rho*v;
+}
+#endif
--- a/src/CelerisLab/lbm/kernels/kernel_v2.cu
+++ b/src/CelerisLab/lbm/kernels/kernel_v2.cu
@@ -56,6 +56,9 @@
#include "boundary/bounce_back.cuh"
#include "boundary/curved_boundary.cuh"
#include "boundary/inlet_outlet.cuh"
+#if Y_WALL_BC == 2
+#include "boundary/planar_velocity.cuh"
+#endif
#include "boundary/ibm_kernels.cuh"
// ============================================================================
--- a/src/CelerisLab/lbm/kernels/step/init_flow.cu
+++ b/src/CelerisLab/lbm/kernels/step/init_flow.cu
@@ -49,10 +49,12 @@
{
if (y == 0u || y == (unsigned int)(NY - 1))
return (uint16_t)(FLAG_SOLID | FLAG_BC_WALL);
+#if Y_WALL_BC != 2
if (x == 0u)
return (uint16_t)(FLAG_SOLID | FLAG_BC_INLET);
if (x == (unsigned int)(NX - 1))
return (uint16_t)(FLAG_SOLID | FLAG_BC_OUTLET);
+#endif
return FLAG_FLUID;
}
@@ -85,7 +87,15 @@
if (is_solid(fl)) {
write_rest_equilibrium(fi, k);
} else {
+#if Y_WALL_BC == 2
+ // Initial g(0)=feq: both boundary velocities are already satisfied.
+ // Start from linear Couette, not the sought exponential steady profile.
+ float f[NQ];
+ compute_feq(RHO, U0*(float)(y-1u)/(float)(NY-3), WALL_NORMAL_VELOCITY, f);
+ stream_pull_store(k, f, fi);
+#else
write_equilibrium(fi, k, inlet_target_u((float)y));
+#endif
}
#elif DIM == 3
unsigned int x, y, z; unsigned long k;
@@ -138,4 +148,4 @@
#endif
}
-#endif // CELERIS_STEP_INIT_FLOW_CU
\ No newline at end of file
+#endif // CELERIS_STEP_INIT_FLOW_CU
--- a/src/CelerisLab/lbm/kernels/step/one_step_double.cu
+++ b/src/CelerisLab/lbm/kernels/step/one_step_double.cu
@@ -88,7 +88,12 @@
bounce_back_swap(f);
if (is_fluid(fl) && (y == 1u || y == (unsigned int)(NY - 2))) {
-#if Y_WALL_BC == 1
+#if Y_WALL_BC == 2
+ // Reconstruct f(t+dt) on the physical plane, then collide below.
+ // This replaces, rather than compounds, the halfway wall correction.
+ if (y == 1u) reconstruct_bottom_velocity(f, 0.0f, WALL_NORMAL_VELOCITY);
+ else reconstruct_top_velocity(f, U0, WALL_NORMAL_VELOCITY);
+#elif Y_WALL_BC == 1
apply_wall_freeslip_d2q9_y_pull(y, f, fi_in, k);
#else
apply_wall_bb_d2q9_y_pull(y, f, fi_in, k);
--- a/src/CelerisLab/simulation.py
+++ b/src/CelerisLab/simulation.py
@@ -115,6 +115,8 @@
Returns:
int body_id (formal id before initialize; -1 if staged).
"""
+ if self.lbm_cfg.y_wall_bc == "on_site_velocity":
+ raise ValueError("on_site_velocity currently requires an object-free domain")
type_lower = str(type).strip().lower()
if type_lower == "circle":
--- /dev/null
+++ b/tests/unit/cpp/planar_velocity.cpp
@@ -0,0 +1,30 @@
+#include <cstdio>
+#include <initializer_list>
+#define __device__
+#define __forceinline__ inline
+#include "boundary/planar_velocity.cuh"
+
+int main() {
+ const int cx[9]={0,1,-1,0,0,1,-1,1,-1};
+ const int cy[9]={0,0,0,1,-1,1,-1,-1,1};
+ const float w[9]={4.f/9,1.f/9,1.f/9,1.f/9,1.f/9,1.f/36,1.f/36,1.f/36,1.f/36};
+ for (int top : {0,1}) for (float rho : {1.f,1.04f})
+ for (float u : {-.02f,0.f,.02f}) for (float v : {-.01f,-.001f,0.f,.001f,.01f})
+ for (int neq : {0,1}) {
+ float f[9];
+ printf("%d %.9g %.9g %.9g %d\n",top,rho,u,v,neq);
+ for (int i=0;i<9;++i) {
+ float cu=cx[i]*u+cy[i]*v;
+ f[i]=w[i]*rho*(1+3*cu+4.5f*cu*cu-1.5f*(u*u+v*v));
+ f[i]+=neq*.0001f*(i-4);
+ // Unknown inputs must be completely overwritten, never consumed.
+ if (cy[i]==(top?-1:1)) f[i]=-7.f-i;
+ printf("%.9g ",f[i]);
+ }
+ puts("");
+ if (top) reconstruct_top_velocity(f,u,v);
+ else reconstruct_bottom_velocity(f,u,v);
+ for (float value : f) printf("%.9g ",value);
+ puts("");
+ }
+}
--- /dev/null
+++ b/tests/unit/test_planar_velocity.py
@@ -0,0 +1,95 @@
+"""J1 actual-header moments and bounded host capability contracts."""
+import json
+from pathlib import Path
+from types import SimpleNamespace
+
+import h5py
+import numpy as np
+import pytest
+
+import CelerisLab
+from CelerisLab.config import LBMConfig, load_lbm_config
+from CelerisLab.cuda import compiler_v2 as compiler
+from CelerisLab.common.checkpoint import save_checkpoint, load_checkpoint
+from .test_numerical_operators import compile_run
+from .test_state_regressions import device, make_field, manager
+
+
+def test_planar_helper_moments_and_known_populations(tmp_path):
+ root = Path(__file__).resolve().parents[2]
+ assert Path(CelerisLab.__file__).resolve() == root/'src/CelerisLab/__init__.py'
+ assert Path(compiler.kernel_path('kernel_v2.cu')).resolve() == root/'src/CelerisLab/lbm/kernels/kernel_v2.cu'
+ lines = compile_run(tmp_path, 'planar_velocity.cpp', {})
+ cx = np.array([0,1,-1,0,0,1,-1,1,-1])
+ cy = np.array([0,0,0,1,-1,1,-1,-1,1])
+ w = np.array([4/9]+[1/9]*4+[1/36]*4)
+ for start in range(0, len(lines), 3):
+ top, rho, u, v, neq = np.fromstring(lines[start], sep=' ')
+ before, after = [np.fromstring(s, sep=' ') for s in lines[start+1:start+3]]
+ known = cy != (-1 if top else 1)
+ np.testing.assert_array_equal(after[known], before[known])
+ density = after.sum()
+ np.testing.assert_allclose([after@cx, after@cy], density*np.array([u,v]), atol=3e-7, rtol=0)
+ if not neq:
+ cu = cx*u+cy*v
+ expected = w*rho*(1+3*cu+4.5*cu**2-1.5*(u*u+v*v))
+ np.testing.assert_allclose(after, expected, atol=3e-7, rtol=0)
+ assert abs(density-rho) <= 3e-7
+
+
+@pytest.mark.parametrize('override', [
+ dict(collision='TRT'), dict(collision='MRT'), dict(streaming='esopull'),
+ dict(store_precision='FP16S'), dict(ddf_shifting=True), dict(les_enabled=True),
+ dict(lattice_model='D3Q19'), dict(wall_normal_velocity=1.),
+ dict(wall_normal_velocity=float('nan')),
+])
+def test_planar_unsupported_configuration(override):
+ with pytest.raises(ValueError):
+ LBMConfig(y_wall_bc='on_site_velocity', **override).validate_runtime()
+
+
+def test_planar_config_roundtrip_and_fixed_targets(tmp_path):
+ path = Path(__file__).resolve().parents[2]/'src/CelerisLab/configs/config_lbm.json'
+ data = json.loads(path.read_text())
+ data['method'].update(y_wall_bc='on_site_velocity', wall_normal_velocity=-.01,
+ collision='SRT', streaming='double_buffer',
+ store_precision='FP32', ddf_shifting=False)
+ data['method']['les']['enabled'] = False
+ path = tmp_path/'config.json'
+ path.write_text(json.dumps(data))
+ cfg = load_lbm_config(str(path))
+ assert float(cfg.to_macros()['WALL_NORMAL_VELOCITY'][:-1]) == -.01
+ with pytest.raises(ValueError, match='object-free'):
+ compiler.generate_config(cfg, n_objects=1, build_root=tmp_path)
+ before = cfg.allocation_signature()
+ cfg.wall_normal_velocity = .01
+ assert cfg.allocation_signature() != before
+
+
+def test_planar_force_rejected_before_state_change(device):
+ field = make_field(device, y_wall_bc='on_site_velocity')
+ before, writes = field._runtime_params.copy(), device.writes
+ with pytest.raises(ValueError, match='body forces'):
+ field.update_params(fx=.001)
+ assert field._runtime_params == before and device.writes == writes
+
+
+@pytest.mark.parametrize('damage', ['target', 'mode', 'force'])
+def test_planar_checkpoint_contract_before_writes(device, tmp_path, damage):
+ field = make_field(device, y_wall_bc='on_site_velocity')
+ bodies = manager(field)
+ bodies.sync_to_gpu(field)
+ path = save_checkpoint(field, SimpleNamespace(step_count=3), field.cfg, bodies,
+ str(tmp_path/'plane.h5'))
+ with h5py.File(path, 'a') as hf:
+ if damage == 'force':
+ hf.attrs['runtime_params_json'] = '{"fx":0.001}'
+ else:
+ cfg = json.loads(hf.attrs['config_json'])
+ cfg['wall_normal_velocity' if damage == 'target' else 'y_wall_bc'] = (
+ .01 if damage == 'target' else 'bounce_back')
+ hf.attrs['config_json'] = json.dumps(cfg)
+ writes = device.writes
+ with pytest.raises(ValueError):
+ load_checkpoint(path, field, SimpleNamespace(_step_count=0), field.cfg, bodies)
+ assert device.writes == writes
--- /dev/null
+++ b/tests/validation/run_j1.py
@@ -0,0 +1,239 @@
+"""Bounded J1 Couette-throughflow validation; gates fixed before first run.
+
+See research/j1-2026-10-12/REVIEW.md. No production audit kernel or extra launch
+is added: a temporary test build records actual incoming populations and the
+roundoff of the real collision, then is compared with uninstrumented stepping.
+"""
+import argparse
+import hashlib
+import json
+import os
+from pathlib import Path
+import sys
+from tempfile import TemporaryDirectory
+from unittest.mock import patch
+
+import numpy as np
+import pycuda.driver as cuda
+
+import CelerisLab
+from CelerisLab.cuda import compiler_v2 as compiler
+from CelerisLab.simulation import Simulation
+
+ROOT = Path(__file__).resolve().parents[2]
+CX = np.array([0,1,-1,0,0,1,-1,1,-1])
+CY = np.array([0,0,0,1,-1,1,-1,-1,1])
+GATES = dict(u_relative_l2=1e-3, v_max_error=2e-6, flux_relative_error=1e-3,
+ mass_offset=1e-3, mass_window_span=6e-6, mass_chunk_change=2e-6,
+ u_change=2e-5, v_change=2e-7, budget_relative_residual=2e-7)
+
+
+def make_sim(v, mode='on_site_velocity'):
+ config = json.loads((ROOT/'src/CelerisLab/configs/config_lbm.json').read_text())
+ config['grid'].update(nx=16, ny=35, nz=1, lattice_model='D2Q9')
+ config['physics'].update(viscosity=.1, velocity=.02, rho=1.)
+ config['method'].update(collision='SRT', streaming='double_buffer',
+ store_precision='FP32', ddf_shifting=False,
+ y_wall_bc=mode, wall_normal_velocity=v)
+ config['method']['les']['enabled'] = False
+ with TemporaryDirectory(prefix='celeris-j1-config-') as folder:
+ lbm, body = Path(folder)/'lbm.json', Path(folder)/'body.json'
+ lbm.write_text(json.dumps(config))
+ body.write_text('{"objects":[]}')
+ return Simulation(str(lbm), str(body), device_id=0)
+
+
+def populations(sim):
+ return sim.get_ddf().reshape(9, sim.field.ny, sim.field.nx).astype(np.float64)
+
+
+def observables(f):
+ interior = f[:, 1:-1, :]
+ rho = interior.sum(axis=0)
+ ux = np.einsum('q,qyx->yx', CX, interior)/rho
+ uy = np.einsum('q,qyx->yx', CY, interior)/rho
+ return rho, ux, uy
+
+
+AUDIT = r'''
+__device__ double j1_incoming[2*NX];
+__device__ double j1_collision[NX*NY];
+__device__ void j1_collide_audit(float* f, float rho, float ux, float uy) {
+ const unsigned x=threadIdx.x+blockIdx.x*blockDim.x, y=blockIdx.y;
+ double before=0.0;
+ for(int i=0;i<NQ;++i) before+=(double)f[i];
+ if(y==1) j1_incoming[x]=(double)f[3]+f[5]+f[8];
+ if(y==NY-2) j1_incoming[NX+x]=(double)f[4]+f[6]+f[7];
+ collide_dispatch(f,rho,ux,uy);
+ double after=0.0;
+ for(int i=0;i<NQ;++i) after+=(double)f[i];
+ j1_collision[y*NX+x]=after-before;
+}
+'''
+
+
+def population_budget(sim):
+ """16 real steps, separate incoming/outgoing/collision records, no targets."""
+ saved = sim.get_ddf().copy()
+ count = sim.stepper.step_count
+ sim.run(16)
+ plain = sim.get_ddf().copy()
+ sim.field.ddf[:] = saved
+ sim.field.upload_ddf()
+ sim.stepper._step_count = count
+ build = compiler.generate_config(sim.lbm_cfg)
+ original_step = sim.stepper.step_fn
+ try:
+ source = Path(build)/'step/one_step_double.cu'
+ code = source.read_text().replace('collide_dispatch(f, rho_n, ux, uy);',
+ 'j1_collide_audit(f, rho_n, ux, uy);')
+ marker = '__global__ __launch_bounds__(NT, 2)'
+ assert code.count(marker) == 1
+ source.write_text(code.replace(marker, AUDIT+'\n'+marker))
+ audit_hash = hashlib.sha256(source.read_bytes()).hexdigest()
+ ptx = compiler.compile_kernel(build_dir=build)
+ rows = []
+ with sim.ctx.activate():
+ module = compiler.load_module(ptx)
+ params, size = module.get_global('d_params')
+ cuda.memcpy_dtod(params, sim.field.module.get_global('d_params')[0], size)
+ incoming_ptr, _ = module.get_global('j1_incoming')
+ collision_ptr, size = module.get_global('j1_collision')
+ cuda.memset_d8(collision_ptr, 0, size)
+ incoming = np.empty(2*sim.field.nx, dtype=np.float64)
+ collision = np.empty(sim.field.n, dtype=np.float64)
+ sim.stepper.step_fn = module.get_function('OneStep')
+ for _ in range(16):
+ before = populations(sim)
+ mass = before[:, 1:-1].sum()
+ outgoing_bottom = before[[4,6,7], 1, :].sum()
+ outgoing_top = before[[3,5,8], -2, :].sum()
+ sim.run(1)
+ after = populations(sim)
+ cuda.memcpy_dtoh(incoming, incoming_ptr)
+ cuda.memcpy_dtoh(collision, collision_ptr)
+ delta = after[:, 1:-1].sum()-mass
+ bottom = incoming[:sim.field.nx].sum()-outgoing_bottom
+ top = incoming[sim.field.nx:].sum()-outgoing_top
+ physical_residual = delta-bottom-top
+ residual = physical_residual-collision.sum()
+ rows.append(dict(step=sim.stepper.step_count, mass=mass, delta_mass=delta,
+ bottom_net_in=bottom, top_net_in=top,
+ collision_roundoff=float(collision.sum()),
+ physical_residual=physical_residual,
+ budget_residual=residual))
+ identical = np.array_equal(sim.get_ddf(), plain)
+ return dict(steps=rows, instrumented_matches_plain=identical,
+ audit_source_sha256=audit_hash,
+ max_relative_residual=max(abs(r['budget_residual'])/r['mass'] for r in rows),
+ max_physical_relative_residual=max(abs(r['physical_residual'])/r['mass'] for r in rows))
+ finally:
+ sim.stepper.step_fn = original_step
+ compiler.cleanup_build_dir(build)
+
+
+def run_case(v):
+ height, nu, speed, nx = 32, .1, .02, 16
+ chunk = int(np.ceil(.5*height**2/(np.pi**2*nu)))
+ with make_sim(v) as sim:
+ sim.initialize()
+ # Reject unsupported entrypoints before observing the reference flow.
+ try:
+ sim.add_body('circle', center=(8,16), radius=2)
+ except ValueError:
+ pass
+ else:
+ raise AssertionError('J1 accepted an object')
+ rho, ux, uy = observables(populations(sim))
+ initial_mass = rho.sum()
+ previous_u, previous_v, previous_mass = ux.mean(axis=1), uy.mean(axis=1), initial_mass
+ history = [dict(step=0, mass=float(initial_mass))]
+ steady = False
+ for sample in range(1, 41):
+ sim.run(chunk)
+ f = populations(sim)
+ if not np.isfinite(f).all():
+ raise AssertionError('Nonfinite J1 populations')
+ rho, ux, uy = observables(f)
+ u, normal, mass = ux.mean(axis=1), uy.mean(axis=1), rho.sum()
+ item = dict(step=sim.stepper.step_count, mass=float(mass),
+ u_change=float(np.sqrt(np.mean((u-previous_u)**2))/speed),
+ v_change=float(np.sqrt(np.mean((normal-previous_v)**2))),
+ mass_change=float(abs(mass-previous_mass)/initial_mass))
+ history.append(item)
+ previous_u, previous_v, previous_mass = u.copy(), normal.copy(), mass
+ tail = history[-4:]
+ steady = sample >= 8 and all(
+ r['u_change'] <= GATES['u_change'] and r['v_change'] <= GATES['v_change']
+ and r['mass_change'] <= GATES['mass_chunk_change'] for r in tail)
+ steady = steady and np.ptp([r['mass'] for r in history[-5:]])/initial_mass <= GATES['mass_window_span']
+ if steady:
+ break
+ y = np.arange(height+1)
+ exact = speed*y/height if v == 0 else speed*np.expm1(v*y/nu)/np.expm1(v*height/nu)
+ relative_l2 = float(np.linalg.norm(u[1:-1]-exact[1:-1])/np.linalg.norm(exact[1:-1]))
+ v_error = float(np.max(np.abs(uy-v)))
+ sections = [8,16,24] # Distances above physical bottom; none is a BC row.
+ flux = [float(np.sum(rho[s]*uy[s])) for s in sections]
+ target_flux = float(rho.mean()*v*nx)
+ flux_error = max(abs(q-target_flux)/abs(target_flux) for q in flux) if v else max(abs(q)/nx for q in flux)
+ mass_span = float(np.ptp([r['mass'] for r in history[-5:]])/initial_mass)
+ offset = float(abs(mass-initial_mass)/initial_mass)
+ result = dict(v=v, height=height, nx=nx, nu=nu, top_u=speed,
+ pe=v*height/nu, re=speed*height/nu,
+ ma=float(np.hypot(speed,v)*np.sqrt(3)),
+ steps=sim.stepper.step_count, max_steps=40*chunk, steady=bool(steady),
+ relative_l2=relative_l2, v_max_error=v_error,
+ mass_offset=offset, mass_window_span=mass_span,
+ rho_range=[float(rho.min()),float(rho.max())],
+ section_y=sections, measured_flux=flux, reference_flux=target_flux,
+ flux_error=flux_error, profile_u=u.tolist(), profile_v=normal.tolist(),
+ analytic_u=exact.tolist(), history=history)
+ budget = population_budget(sim)
+ result['population_budget'] = budget
+ result['passed'] = bool(steady and relative_l2 <= GATES['u_relative_l2']
+ and v_error <= GATES['v_max_error']
+ and flux_error <= (GATES['flux_relative_error'] if v else GATES['v_max_error'])
+ and offset <= GATES['mass_offset'] and mass_span <= GATES['mass_window_span']
+ and budget['instrumented_matches_plain']
+ and budget['max_relative_residual'] <= GATES['budget_relative_residual'])
+ return result
+
+
+def off_path_check(baseline):
+ outcomes = []
+ for kernel_root in (baseline, Path(compiler.kernel_path(''))):
+ with patch.object(compiler, '_KERNEL_DIR', str(kernel_root)):
+ with make_sim(0., mode='bounce_back') as sim:
+ sim.initialize()
+ sim.run(7)
+ outcomes.append((sim.get_ddf().copy(), sim.get_flags().copy()))
+ return bool(all(np.array_equal(a,b) for a,b in zip(*outcomes)))
+
+
+def main():
+ parser = argparse.ArgumentParser(description=__doc__)
+ parser.add_argument('--output', type=Path, required=True)
+ parser.add_argument('--baseline-kernels', type=Path,
+ help='Optional pre-change CUDA source for the off-mode comparison')
+ args = parser.parse_args()
+ assert Path.cwd() == ROOT == Path('/home/frank14f/CelerisLab')
+ assert Path(CelerisLab.__file__).resolve() == ROOT/'src/CelerisLab/__init__.py'
+ assert Path(compiler.kernel_path('kernel_v2.cu')).resolve() == ROOT/'src/CelerisLab/lbm/kernels/kernel_v2.cu'
+ assert os.environ.get('CUDA_VISIBLE_DEVICES') in ('0','1')
+ report = dict(python=sys.executable, package=CelerisLab.__file__,
+ cuda_source=compiler.kernel_path('kernel_v2.cu'),
+ physical_gpu=os.environ['CUDA_VISIBLE_DEVICES'], gates=GATES, cases=[])
+ report['off_path_bitwise_equal'] = (off_path_check(args.baseline_kernels)
+ if args.baseline_kernels else None)
+ for v in (0., .001, -.001, .01, -.01):
+ case = run_case(v)
+ report['cases'].append(case)
+ args.output.write_text(json.dumps(report, indent=2)+'\n')
+ print(json.dumps({k:case[k] for k in ('v','steps','steady','relative_l2','v_max_error',
+ 'mass_offset','mass_window_span','flux_error','passed')}), flush=True)
+ return 0 if report['off_path_bitwise_equal'] is not False and all(c['passed'] for c in report['cases']) else 1
+
+
+if __name__ == '__main__':
+ raise SystemExit(main())