Retain shifted density perturbation in opt-in blocked collision

tests/unit/cpp/shifted_block_operators.cpp

@@ -0,0 +1,36 @@
+#include <cmath>
+#include <cstdio>
+#include <initializer_list>
+#define __device__
+#define __constant__
+#define __forceinline__ inline
+#define LATTICE_D2Q9 0
+#define LATTICE_MODEL 0
+#define NQ 9
+#define DIM 2
+#define USE_DDF_SHIFTING 1
+#define WALL_NORMAL_VELOCITY 0.f
+#define WALL_SLOT_WIDTH 0
+#define NX 16
+#include "core/descriptors.cuh"
+#include "operators/macro.cuh"
+#include "operators/collision_srt.cuh"
+#include "boundary/planar_velocity.cuh"
+void emit(const float* f) {for(int q=0;q<9;q++) printf("%.9g ",f[q]);puts("");}
+int main() {
+    for(float delta : {0.f, 0x1p-30f, -0x1p-30f, .0002f})
+    for(float u : {0.f,.005f}) {
+        float eq[9]; compute_feq(1.f+delta,u,-.2f*u,eq,delta);
+        printf("E %.9g %.9g %.9g\n",delta,u,-.2f*u);emit(eq);
+    }
+    for(int mode : {0,1,2}) for(float rho : {1.f,1.04f}) {
+        float f[9];const float u=.007f,v=.003f;
+        compute_feq(rho,u,v,f,rho-1.f);
+        for(int q=0;q<9;q++)f[q]+=.00001f*(q-4);
+        printf("B %d %.9g %.9g %.9g\n",mode,rho,u,v);emit(f);
+        if(mode==0)reconstruct_bottom_velocity(f,u,v);
+        else if(mode==1)reconstruct_top_velocity(f,u,v);
+        else reconstruct_top_reservoir(f,rho);
+        emit(f);
+    }
+}

tests/unit/test_shifted_block_operators.py

@@ -0,0 +1,32 @@
+"""Actual-header CPU algebra, not CUDA arithmetic/performance evidence."""
+import numpy as np
+from .test_numerical_operators import compile_run
+
+
+def test_retained_density_and_shifted_planar_headers(tmp_path):
+    rows=compile_run(tmp_path,'shifted_block_operators.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)
+    it=iter(rows)
+    for header in it:
+        kind,*numbers=header.split()
+        a=np.fromstring(next(it),sep=' ',dtype=np.float32).astype(float)
+        if kind=='E':
+            delta,u,v=np.array(numbers,dtype=np.float32).astype(float)
+            rho=float(np.float32(1+delta));cu=3*(cx*u+cy*v)
+            expected=w*(delta+rho*(cu+.5*cu*cu-1.5*(u*u+v*v)))
+            np.testing.assert_allclose(a,expected,rtol=3e-6,atol=1e-10)
+            if u==0:
+                np.testing.assert_allclose(a.sum(),delta,rtol=5e-8,atol=0)
+                if delta: assert np.any(a!=0)
+        else:
+            mode=int(numbers[0]);rho,u,v=np.array(numbers[1:],dtype=np.float32).astype(float)
+            b=np.fromstring(next(it),sep=' ',dtype=np.float32).astype(float)
+            known=cy!=(1 if mode==0 else -1)
+            np.testing.assert_array_equal(a[known],b[known])
+            physical_rho=1+b.sum()
+            if mode==2:
+                assert abs(physical_rho-rho)<2e-8
+                assert abs(b@cx)<2e-9
+            else:
+                np.testing.assert_allclose([b@cx,b@cy],physical_rho*np.array([u,v]),atol=3e-9,rtol=0)