Frank14f/CelerisLab
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
codex/block-root-bc
CelerisLab/research/surface_actuation_20261014/incremental.diff
T

1066 lines
49 KiB
Diff
Raw Permalink Blame History

--- a/src/CelerisLab/body/__init__.py
+++ b/src/CelerisLab/body/__init__.py
@@ -2,7 +2,7 @@
"""
Body / object management for immersed and rigid objects.
"""
-from .objects import SimObject, ObjectState, ObjectControl
+from .objects import SimObject, ObjectState, ObjectControl, SurfaceActuation
from importlib import import_module
from .registry import BodyRegistry
from .sync_plan import BodySyncPlan
@@ -11,7 +11,7 @@
from .coupling.soa_packer import pack_cut_links_to_soa, pack_sensor_to_soa
__all__ = [
- "SimObject", "ObjectState", "ObjectControl",
+ "SimObject", "ObjectState", "ObjectControl", "SurfaceActuation",
"ObjectManager",
"BodyRegistry",
"BodySyncPlan",
--- a/src/CelerisLab/body/manager.py
+++ b/src/CelerisLab/body/manager.py
@@ -108,6 +108,7 @@
# -- Registry delegation -------------------------------------------------
def add(self, obj: SimObject) -> int:
"""Register an object and return its id."""
+ obj.validate_surface_actuation(self.cfg)
body_id = self._registry.add(obj)
self._resize_buffers()
return body_id
@@ -160,6 +161,7 @@
The object is not added to the formal registry until ``commit_pending()``
is called. This implicitly activates the edit window.
"""
+ obj.validate_surface_actuation(self.cfg)
self._edit_active = True
self._pending_add.append(obj)
@@ -201,7 +203,8 @@
new = SimObject(obj_id=oid, geometry=obj.geometry,
center=obj.center, radius=obj.radius,
is_sensor=obj.is_sensor,
- is_force_region=obj.is_force_region)
+ is_force_region=obj.is_force_region,
+ surface_actuation=obj.surface_actuation)
new.state = obj.state
new.control = obj.control
return new
@@ -306,6 +309,7 @@
cl_host = plan.curved_host
field.curved.assign_host(*cl_host)
+ field.curved.assign_surface(plan.next_objects, self.cfg, plan.next_flags)
s_cells, s_ids = plan.sensor_host
field.sensors.assign_host(s_cells, s_ids)
fr_cells, fr_ids = plan.force_region_host
@@ -392,6 +396,7 @@
ex, ey = D2Q9_EX, D2Q9_EY
for obj in objects:
+ obj.validate_surface_actuation(self.cfg)
if (not obj.is_sensor) and (not obj.is_force_region) and hasattr(obj, 'get_curved_list'):
(
fluid_idx, dirs, q_vals, body_ids,
@@ -447,6 +452,13 @@
r_fr_cells = np.concatenate(fr_cells) if fr_cells else np.zeros(0, dtype=np.uint32)
r_fr_ids = np.concatenate(fr_ids) if fr_ids else np.zeros(0, dtype=np.int32)
+ if any(obj.surface_actuation is not None for obj in objects):
+ # Reject unresolved apertures/ownership before sync or restart mutates state.
+ from ..lbm.curved_links import CurvedLinkSoA
+ links = CurvedLinkSoA()
+ links.assign_host(r_idx, r_dir, r_q, r_bid, r_rx, r_ry, r_rz, r_fallback)
+ links.assign_surface(objects, self.cfg, domain_flags)
+
return (
r_idx, r_dir, r_q, r_bid,
r_rx, r_ry, r_rz, r_fallback,
@@ -466,21 +478,24 @@
"""
self._telemetry_field = field
- if rebuild_flags:
- field.flag = self.build_flags(field.build_channel_flags())
- field.upload_flags()
+ next_flags = self.build_flags(field.build_channel_flags()) if rebuild_flags else field.flag
(
cl_fluid_idx, cl_dir, cl_q, cl_body_id,
cl_rx, cl_ry, cl_rz, cl_fallback_class,
sensor_cells, sensor_obj_id,
fr_cells, fr_obj_id,
- ) = self.build_compact_lists(domain_flags=field.flag)
+ ) = self.build_compact_lists(domain_flags=next_flags)
+
+ if rebuild_flags:
+ field.flag = next_flags
+ field.upload_flags()
field.curved.assign_host(
cl_fluid_idx, cl_dir, cl_q, cl_body_id,
cl_rx, cl_ry, cl_rz, cl_fallback_class,
)
+ field.curved.assign_surface(self.objects, self.cfg, field.flag)
field.sensors.assign_host(sensor_cells, sensor_obj_id)
field.force_regions.assign_host(fr_cells, fr_obj_id)
@@ -595,6 +610,10 @@
def _upload_action_async(self, stream: cuda.Stream) -> None:
"""Async H2D of host action array to action_gpu."""
+ for obj in self.objects:
+ obj.validate_surface_actuation(self.cfg)
+ if obj.surface_actuation is not None and self.action[(obj.obj_id+1)*3*self.cfg.dim-1] != 0:
+ raise ValueError("Surface actuation cannot share the rotation action")
if self.action_gpu is None or self.action.size == 0:
return
if self._action_pinned is None or self._action_pinned.size != self.action.size:
@@ -712,6 +731,8 @@
The GPU action buffer is updated at the next ``run()`` call.
"""
self._validate_body_id(body_id)
+ obj = self.get(body_id)
+ obj.validate_surface_actuation(self.cfg, dict(vars(obj.state), omega=omega))
dim = self.cfg.dim if self.cfg else 2
slot = 3 * dim
base = body_id * slot
--- a/src/CelerisLab/body/objects.py
+++ b/src/CelerisLab/body/objects.py
@@ -15,6 +15,7 @@
from __future__ import annotations
from dataclasses import dataclass
+import math
from typing import Optional, Tuple
import numpy as np
@@ -53,6 +54,37 @@
mode: int = 0 # 0=force, 1=velocity, 2=displacement
+@dataclass(frozen=True)
+class SurfaceActuation:
+ """Fixed circular aperture: angles in radians, velocity in lattice units.
+
+ ``angle`` is measured counterclockwise from +x; ``half_width`` is the arc
+ half-angle. Positive normal_velocity blows from solid into fluid. The C1
+ profile (1-s*s)^2 uses the true hit angle, not a fluid-node aperture mask.
+ This immutable specification is set before initialization; no moving nozzle.
+ """
+ angle: float
+ half_width: float
+ normal_velocity: float
+
+ def __post_init__(self):
+ if not all(math.isfinite(v) for v in (self.angle, self.half_width, self.normal_velocity)):
+ raise ValueError("Surface actuation values must be finite")
+ if not 0 < self.half_width <= math.pi or abs(self.normal_velocity) >= 1:
+ raise ValueError("Require 0 < half_width <= pi and abs(normal_velocity) < 1")
+
+ def project(self, rx, ry):
+ """Project stored true hit-center vectors; returns velocity and aperture mask."""
+ radius = np.hypot(rx, ry)
+ if np.any(radius <= 0):
+ raise ValueError("Surface normal requires nonzero hit-center distance")
+ delta = np.arctan2(np.sin(np.arctan2(ry, rx)-self.angle),
+ np.cos(np.arctan2(ry, rx)-self.angle))
+ aperture = np.abs(delta) < self.half_width
+ speed = self.normal_velocity*np.maximum(0., 1-(delta/self.half_width)**2)**2
+ return np.stack((speed*rx/radius, speed*ry/radius), axis=-1).astype(np.float32), aperture
+
+
# ---------------------------------------------------------------------------
# SimObject container
# ---------------------------------------------------------------------------
@@ -77,6 +109,7 @@
radius: float = 0.0,
is_sensor: bool = False,
is_force_region: bool = False,
+ surface_actuation: Optional[SurfaceActuation] = None,
) -> None:
self.obj_id = obj_id
self.center = center if center is not None else (0.0, 0.0)
@@ -84,6 +117,10 @@
self._geometry = geometry
self._is_sensor = is_sensor
self._is_force_region = is_force_region
+ self._surface_actuation = (SurfaceActuation(**surface_actuation)
+ if isinstance(surface_actuation, dict) else surface_actuation)
+ if self._surface_actuation is not None and not isinstance(self._surface_actuation, SurfaceActuation):
+ raise TypeError("surface_actuation must be a SurfaceActuation or dict")
c = self.center
z0 = float(c[2]) if len(c) > 2 else 0.0
@@ -91,6 +128,22 @@
x=float(c[0]), y=float(c[1]), z=z0,
)
self.control = ObjectControl()
+
+ @property
+ def surface_actuation(self):
+ return self._surface_actuation
+
+ def validate_surface_actuation(self, cfg, state=None):
+ if self.surface_actuation is None:
+ return
+ from .geometry.circle import CircleGeometry
+ if (type(self.geometry) is not CircleGeometry or self.is_sensor or self.is_force_region
+ or cfg is None or not cfg.is_d2q9 or cfg.store_precision != "FP32"
+ or cfg.streaming != "double_buffer" or cfg.ddf_shifting):
+ raise ValueError("Surface actuation requires a circle and D2Q9/FP32/unshifted/double_buffer")
+ values = vars(self.state) if state is None else state
+ if any(values.get(k, 0.) != 0 for k in ("vx", "vy", "omega")):
+ raise ValueError("Surface actuation requires a stationary circle")
@property
def geometry(self) -> Geometry:
--- a/src/CelerisLab/common/checkpoint.py
+++ b/src/CelerisLab/common/checkpoint.py
@@ -24,7 +24,7 @@
from ..lbm.descriptors import FLUID, SOLID, OBSTACLE, SENSOR_FLAG, FRC_REGION
# File format version — bump on incompatible layout changes.
-CKPT_VERSION = 3
+CKPT_VERSION = 4
_SUPPORTED_STORE_PRECISIONS = ("FP32", "FP16S")
@@ -130,10 +130,20 @@
if not np.array_equal(flag_buf, expected_flags):
raise ValueError("Checkpoint flags do not match current object/domain geometry")
bodies.build_compact_lists_for(bodies.objects, domain_flags=expected_flags)
+ applied_action = None
+ if bodies.action_gpu is not None:
+ applied_action = np.empty_like(bodies.action)
+ cuda.memcpy_dtoh(applied_action, bodies.action_gpu)
+ for obj in bodies.objects:
+ if (obj.surface_actuation is not None
+ and applied_action[(obj.obj_id + 1) * 3 * lbm_cfg.dim - 1] != 0):
+ raise ValueError("Surface actuation requires zero applied omega")
with h5py.File(path, "w") as hf:
# Version and metadata
- hf.attrs["ckpt_version"] = CKPT_VERSION
+ # Older readers must not silently restart a driven surface as a wall.
+ hf.attrs["ckpt_version"] = (CKPT_VERSION if any(
+ obj.surface_actuation is not None for obj in bodies.objects) else 3)
hf.attrs["celerislab_version"] = __import__(
"CelerisLab").__version__
hf.attrs["timestamp"] = time.strftime("%Y-%m-%dT%H:%M:%S")
@@ -156,11 +166,13 @@
compression_opts=4)
hf.attrs["geometry_json"] = json.dumps(geometry)
+ hf.attrs["surface_actuation_json"] = json.dumps([
+ asdict(obj.surface_actuation) if obj.surface_actuation is not None else None
+ for obj in bodies.objects
+ ])
hf.attrs["runtime_params_json"] = json.dumps(field._runtime_params)
hf.attrs["obs_accum_steps"] = bodies._obs_accum_steps
- if bodies.action_gpu is not None:
- applied_action = np.empty_like(bodies.action)
- cuda.memcpy_dtoh(applied_action, bodies.action_gpu)
+ if applied_action is not None:
hf.create_dataset("applied_action", data=applied_action)
if bodies.obs_gpu is not None:
obs = np.empty(bodies.obs_total_floats, dtype=np.float32)
@@ -191,7 +203,8 @@
store_precision must match the checkpoint. Other parameters
(collision, LES, …) may differ for parameter-study restarts.
v3 restores the saved domain boundaries; v1/v2 retain their default-channel
- validation. Object geometry and roles must match in every version.
+ validation. v4 also requires matching static surface actuation. Object
+ geometry and roles must match in every version.
Args:
path: HDF5 checkpoint file path.
@@ -206,7 +219,7 @@
with h5py.File(path, "r") as hf:
# Version check
ver = hf.attrs.get("ckpt_version", 0)
- if ver not in (1, 2, CKPT_VERSION):
+ if ver not in (1, 2, 3, CKPT_VERSION):
raise ValueError(
f"Checkpoint version mismatch: file={ver}, "
f"expected={CKPT_VERSION}. "
@@ -253,6 +266,14 @@
if (len(obj_data) != bodies.count
or [od.get("obj_id") for od in obj_data] != list(range(bodies.count))):
raise ValueError("Checkpoint object IDs must match the dense current registry")
+ surface_specs = [asdict(obj.surface_actuation) if obj.surface_actuation is not None else None
+ for obj in bodies.objects]
+ if ver >= 4 and "surface_actuation_json" not in hf.attrs:
+ raise ValueError("Checkpoint is missing surface actuation specification")
+ saved_surface_specs = json.loads(hf.attrs.get(
+ "surface_actuation_json", json.dumps([None] * bodies.count)))
+ if saved_surface_specs != surface_specs:
+ raise ValueError("Checkpoint surface actuation mismatch")
domain_flags = None
if ver >= 3:
if "domain_flag" not in hf:
@@ -319,6 +340,11 @@
modes = [int(od.get("control", {}).get("mode", 0)) for od in obj_data]
if not all(np.isfinite(value) for record in states + controls for value in record.values()):
raise ValueError("Checkpoint object state/control must be finite")
+ for obj, state in zip(bodies.objects, states):
+ obj.validate_surface_actuation(lbm_cfg, state)
+ if (obj.surface_actuation is not None and applied_action is not None
+ and applied_action[(obj.obj_id + 1) * 3 * lbm_cfg.dim - 1] != 0):
+ raise ValueError("Checkpoint surface actuation requires zero applied omega")
for obj, state, control, mode in zip(bodies.objects, states, controls, modes):
for key, value in state.items():
--- a/src/CelerisLab/lbm/curved_links.py
+++ b/src/CelerisLab/lbm/curved_links.py
@@ -74,6 +74,48 @@
count: int = 0
capacity: int = 0
+
+ # Optional fixed surface velocity derived from true hit-center columns.
+ surface_velocity: Optional[np.ndarray] = None
+ aperture: Optional[np.ndarray] = None
+ surface_velocity_gpu: Optional[cuda.DeviceAllocation] = None
+ link_exchange_gpu: Optional[cuda.DeviceAllocation] = None
+
+ def assign_surface(self, objects, cfg, flags):
+ """Keep the original eight-column geometry/checkpoint contract unchanged."""
+ enabled = [obj for obj in objects if obj.surface_actuation is not None]
+ self.surface_velocity = self.aperture = None
+ if not enabled:
+ return
+ from .descriptors import FLUID, SOLID
+ if np.any((flags[self.fluid_idx] & FLUID) == 0) or np.any(flags[self.fluid_idx] & SOLID):
+ raise ValueError("Surface cut links must have active fluid receivers")
+ if len(np.unique(self.fluid_idx.astype(np.uint64)*cfg.nq+self.dir)) != len(self.dir):
+ raise ValueError("Surface boundary requires unique cut-link writers")
+ self.surface_velocity = np.zeros((len(self.dir), 2), np.float32)
+ self.aperture = np.zeros(len(self.dir), bool)
+ for obj in enabled:
+ obj.validate_surface_actuation(cfg)
+ selected = self.body_id == obj.obj_id
+ velocity, aperture = obj.surface_actuation.project(self.rx[selected], self.ry[selected])
+ if not aperture.any():
+ raise ValueError("Surface aperture hits no resolved cut links")
+ self.surface_velocity[selected] = velocity
+ self.aperture[selected] = aperture
+
+ def surface_exchange(self, body_id):
+ """Latest boundary step's (Fx,Fy,Tz), split by aperture; not nozzle thrust.
+
+ Caller synchronizes the CUDA stream/context. Rebuild/restart clears this
+ diagnostic; it is not the time-accumulated rigid-body observation.
+ """
+ if self.link_exchange_gpu is None:
+ raise RuntimeError("No surface actuation diagnostic is allocated")
+ values = np.empty((self.count, 3), np.float32)
+ cuda.memcpy_dtoh(values, self.link_exchange_gpu)
+ body = self.body_id == body_id
+ return {name: values[body & mask].sum(axis=0, dtype=np.float64)
+ for name, mask in (("wall", ~self.aperture), ("aperture", self.aperture))}
def assign_host(
self,
@@ -144,6 +186,7 @@
raise ValueError(
"CurvedLinkSoA column 'fallback_class' must use classes {0, 1}"
)
+ self.surface_velocity = self.aperture = None
def upload(self, stream: Optional[cuda.Stream] = None) -> None:
"""Upload host columns to GPU with capacity-based allocation reuse.
@@ -155,6 +198,9 @@
n = int(len(self.fluid_idx))
self.count = n
if n == 0:
+ _free_alloc(self.surface_velocity_gpu)
+ _free_alloc(self.link_exchange_gpu)
+ self.surface_velocity_gpu = self.link_exchange_gpu = None
return
if n > self.capacity or self.fluid_idx_gpu is None:
@@ -190,6 +236,14 @@
(self.rz_gpu, self.rz),
(self.fallback_class_gpu, self.fallback_class),
], stream)
+ _free_alloc(self.surface_velocity_gpu)
+ _free_alloc(self.link_exchange_gpu)
+ self.surface_velocity_gpu = self.link_exchange_gpu = None
+ if self.surface_velocity is not None:
+ self.surface_velocity_gpu = cuda.mem_alloc(self.surface_velocity.nbytes)
+ self.link_exchange_gpu = cuda.mem_alloc(n*3*4)
+ _upload_columns([(self.surface_velocity_gpu, self.surface_velocity),
+ (self.link_exchange_gpu, np.zeros((n, 3), np.float32))], stream)
def free_gpu_columns(self) -> None:
"""Release device columns only (host arrays unchanged)."""
@@ -201,6 +255,9 @@
_free_alloc(self.ry_gpu)
_free_alloc(self.rz_gpu)
_free_alloc(self.fallback_class_gpu)
+ _free_alloc(self.surface_velocity_gpu)
+ _free_alloc(self.link_exchange_gpu)
+ self.surface_velocity_gpu = self.link_exchange_gpu = None
self.fluid_idx_gpu = None
self.dir_gpu = None
self.q_gpu = None
--- a/src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh
+++ b/src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh
@@ -82,7 +82,7 @@
float Uw, float Vw,
float rx, float ry,
unsigned char fallback_class,
- float* __restrict__ obs, int body_id)
+ float* __restrict__ obs, int body_id, float* link_exchange = nullptr)
{
const unsigned int dir_opp = (unsigned int)opp_dir((int)dir);
@@ -120,7 +120,7 @@
// corrected incoming population into the adjacent fluid node.
store_ddf(fi, index_f(k_s, dir_opp), f_reflected);
- if (obs != nullptr) {
+ if (obs != nullptr || link_exchange != nullptr) {
#if USE_DDF_SHIFTING
const float w_i = (float)d_w[dir];
const float fx = (float)d_cx[dir] * ((f_toward + w_i) + (f_reflected + w_i));
@@ -129,10 +129,15 @@
const float fx = (float)d_cx[dir] * (f_toward + f_reflected);
const float fy = (float)d_cy[dir] * (f_toward + f_reflected);
#endif
- atomicAdd(&obs[obs_force_index(body_id, 0)], fx);
- atomicAdd(&obs[obs_force_index(body_id, 1)], fy);
const float tz = rx * fy - ry * fx;
- atomicAdd(&obs[obs_torque_index(body_id, 0)], tz);
+ if (link_exchange != nullptr) {
+ link_exchange[0] = fx; link_exchange[1] = fy; link_exchange[2] = tz;
+ }
+ if (obs != nullptr) {
+ atomicAdd(&obs[obs_force_index(body_id, 0)], fx);
+ atomicAdd(&obs[obs_force_index(body_id, 1)], fy);
+ atomicAdd(&obs[obs_torque_index(body_id, 0)], tz);
+ }
}
}
#endif // DIM == 2
--- a/src/CelerisLab/lbm/kernels/step/aux_kernels.cu
+++ b/src/CelerisLab/lbm/kernels/step/aux_kernels.cu
@@ -43,7 +43,8 @@
const int* cl_body_id,
const float* action,
float* obs,
- unsigned int n_curved)
+ unsigned int n_curved,
+ const float* surface_velocity = nullptr, float* link_exchange = nullptr)
{
unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
if (tid >= n_curved) return;
@@ -58,10 +59,16 @@
const float omega = action_omega(action, bid);
#if DIM == 2
- const float Uw = -omega * ry;
- const float Vw = omega * rx;
+ float Uw = -omega * ry;
+ float Vw = omega * rx;
+ if (surface_velocity != nullptr) {
+ // Independent prescribed fluid velocity on fixed geometry.
+ const float u = surface_velocity[2*tid], v = surface_velocity[2*tid+1];
+ if (u != 0.0f || v != 0.0f) { Uw += u; Vw += v; }
+ }
(void)cl_rz;
- apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, rx, ry, fallback_class, obs, bid);
+ apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, rx, ry, fallback_class, obs, bid,
+ link_exchange == nullptr ? nullptr : link_exchange + 3*tid);
#elif DIM == 3
const float rz = cl_rz[tid];
const float Uw = -omega * ry;
--- a/src/CelerisLab/lbm/stepper.py
+++ b/src/CelerisLab/lbm/stepper.py
@@ -153,6 +153,7 @@
cl.rx_gpu, cl.ry_gpu, cl.rz_gpu, cl.fallback_class_gpu, cl.body_id_gpu,
action_gpu, obs_gpu,
np.uint32(f.n_curved),
+ cl.surface_velocity_gpu or np.uintp(0), cl.link_exchange_gpu or np.uintp(0),
block=(tpb, 1, 1), grid=grid_c,
**launch_kw,
)
--- a/src/CelerisLab/simulation.py
+++ b/src/CelerisLab/simulation.py
@@ -98,7 +98,7 @@
# -- Object management ---------------------------------------------------
def add_body(self, type: str = "circle", center=None,
- radius: float = 0.0, **kwargs) -> int:
+ radius: float = 0.0, *, surface_actuation=None, **kwargs) -> int:
"""Add a simulation body. Returns body_id.
Before ``initialize()`` the object is added to the formal registry
@@ -110,6 +110,10 @@
(future: ``"polygon"``, ``"mesh"``).
center: (x, y) center coordinates (lattice units).
radius: body radius (lattice units).
+ surface_actuation: Fixed circle aperture dict/SurfaceActuation:
+ angle, half_width (radians), normal_velocity (positive outward).
+ D2Q9/FP32/unshifted/DB only; set before initialization or stage
+ a new body. No runtime motion or amplitude setter.
**kwargs: passed to the geometry constructor.
Returns:
@@ -118,6 +122,8 @@
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 surface_actuation is not None and type_lower != "circle":
+ raise ValueError("Surface actuation is only supported on circles")
if type_lower == "circle":
from .body.geometry.circle import CircleGeometry
@@ -125,7 +131,7 @@
geom = CircleGeometry(cx, cy, radius)
obj = SimObject(obj_id=-1, geometry=geom,
center=center, radius=radius,
- is_sensor=False)
+ is_sensor=False, surface_actuation=surface_actuation)
elif type_lower == "sensor":
from .body.geometry.circle import CircleGeometry
cx, cy = float(center[0]), float(center[1])
@@ -189,6 +195,8 @@
center = tuple(spec.get("center", ()))
radius = float(spec.get("radius", 0.0))
omega = spec.get("omega", None)
+ if spec.get("surface_actuation") is not None and obj_type != "cylinder":
+ raise ValueError("Surface actuation is only supported on cylinders in body config")
if obj_type == "":
raise ValueError(
@@ -200,9 +208,10 @@
)
if obj_type == "cylinder":
- body_id = self.add_body("circle", center=center, radius=radius)
+ body_id = self.add_body("circle", center=center, radius=radius,
+ surface_actuation=spec.get("surface_actuation"))
if omega is not None:
- self.bodies.get(body_id).state.omega = float(omega)
+ self.bodies.set_body_state(body_id, omega=float(omega))
elif obj_type == "sensor":
self.add_body("sensor", center=center, radius=radius)
else:
@@ -225,6 +234,10 @@
vy: translation velocity y (future). ``None`` = no change.
"""
obj = self.bodies.get(id)
+ proposed = dict(vars(obj.state))
+ proposed.update({k: float(v) for k, v in (("omega", omega), ("vx", vx), ("vy", vy))
+ if v is not None})
+ obj.validate_surface_actuation(self.lbm_cfg, proposed)
changed = False
if omega is not None:
obj.state.omega = float(omega)
@@ -290,6 +303,18 @@
if self.bodies.obs_pinned is None:
raise RuntimeError("No obs buffer. Call run() first.")
return self.bodies.read_sensor(id, normalize=normalize)
+
+ @context_bound
+ def read_surface_exchange(self, id: int):
+ """Latest step's wall/aperture (Fx,Fy,Tz) exchange, not total nozzle force.
+
+ Unlike read_body this is not accumulated/time-averaged. It resets after
+ geometry sync or checkpoint load; advance before interpreting it.
+ """
+ if self.bodies.get(id).surface_actuation is None:
+ raise ValueError("Body has no surface actuation")
+ self.stream.synchronize()
+ return self.field.curved.surface_exchange(id)
@context_bound
def read_body(self, id: int, *,
--- /dev/null
+++ b/tests/integration/test_surface_actuation_gpu.py
@@ -0,0 +1,249 @@
+"""Fixed surface actuation: discrete operator/ledger checks, not flux accuracy."""
+from dataclasses import asdict
+import json
+
+import h5py
+import numpy as np
+import pycuda.driver as cuda
+import pytest
+
+from CelerisLab.lbm.descriptors import FLUID, SOLID
+from ._cases import make_sim
+from .test_numerical_gpu import CX, CY, install_spatial_state, spatial_state
+
+
+WEIGHTS = np.array([4/9] + [1/9]*4 + [1/36]*4)
+OPPOSITE = np.array([0, 2, 1, 4, 3, 6, 5, 8, 7])
+GEOMETRY = ("fluid_idx", "dir", "q", "body_id", "rx", "ry", "rz", "fallback_class")
+
+
+def feq(rho, ux, uy):
+ """Independent FP64 polynomial; round only the uploaded populations."""
+ cu = CX[:, None, None]*ux + CY[:, None, None]*uy
+ return WEIGHTS[:, None, None]*rho*(1 + 3*cu + 4.5*cu**2 - 1.5*(ux**2 + uy**2))
+
+
+@pytest.mark.parametrize("speed", [-.01, .01])
+def test_actual_surface_kernel_equilibrium_branches(speed):
+ """Both q branches, q=.5 and donor fallback; non-unit rho and signed normal."""
+ sim = make_sim(nx=24, ny=16, velocity=0)
+ allocations = []
+ try:
+ sim.initialize()
+ normal = np.array([-.6, .8])
+ expected = feq(1.17, *(speed*normal))[:, 0, 0]
+ initial = np.broadcast_to(expected[:, None, None], (9, 16, 24)).astype(np.float32).copy()
+ # Distinct source/destination neighborhoods; no synthetic link race.
+ cells = np.array([4*24+12, 6*24+12, 8*24+12, 10*24+12], np.uint32)
+ directions = np.ones(4, np.uint8)
+ radii = np.tile(5*normal, (4, 1)).astype(np.float32)
+ columns = [cells, directions, np.array([.2, .5, .8, .2], np.float32),
+ radii[:, 0].copy(), radii[:, 1].copy(), np.zeros(4, np.float32),
+ np.array([0, 0, 0, 1], np.uint8), np.zeros(4, np.int32)]
+ with sim.ctx.activate():
+ for array in columns + [np.tile(speed*normal, (4, 1)).astype(np.float32),
+ np.zeros((4, 3), np.float32)]:
+ allocations.append(cuda.to_device(array))
+ cuda.memcpy_htod(sim.field.ddf_gpu, initial)
+ sim.stepper.curved_fn(sim.field.ddf_gpu, *allocations[:8],
+ np.uintp(0), np.uintp(0), np.uint32(4),
+ *allocations[8:], block=(32, 1, 1), grid=(1, 1))
+ actual = np.empty_like(initial)
+ exchange = np.empty((4, 3), np.float32)
+ cuda.memcpy_dtoh(actual, sim.field.ddf_gpu)
+ cuda.memcpy_dtoh(exchange, allocations[-1])
+ reflected = actual.reshape(9, -1)[2, cells+1]
+ np.testing.assert_allclose(reflected, expected[2], atol=7e-8, rtol=0)
+ # Equilibrium identity is independent of interpolation implementation.
+ fx = expected[1] + expected[2]
+ wanted = np.column_stack((np.full(4, fx), np.zeros(4), -radii[:, 1]*fx))
+ np.testing.assert_allclose(exchange, wanted, atol=4e-7, rtol=0)
+ finally:
+ with sim.ctx.activate():
+ for allocation in allocations:
+ allocation.free()
+ sim.close()
+
+
+def circle_sim(speed):
+ sim = make_sim(nx=32, ny=24, velocity=0)
+ spec = None if speed is None else dict(angle=.7, half_width=.85, normal_velocity=speed)
+ try:
+ bid = sim.add_body("circle", center=(16, 12), radius=4.3, surface_actuation=spec)
+ sim.initialize()
+ flags = sim.field.flag.copy()
+ flags.reshape(24, 32)[1:-1, [0, -1]] = FLUID # x periodic, y halfway walls
+ y, x = np.mgrid[:24, :32]
+ initial = feq(1.17+.001*np.cos(.2*x+.3*y),
+ .003*np.sin(.3*y), .002*np.cos(.4*x)).astype(np.float32)
+ install_spatial_state(sim, initial)
+ sim.field.flag[:] = flags # installer intentionally resets all flags
+ sim.field.upload_flags()
+ return sim, bid, initial
+ except BaseException:
+ sim.close()
+ raise
+
+
+@pytest.mark.parametrize("speed", [-.01, .01])
+def test_surface_population_ledger_and_public_diagnostics(speed):
+ sim, bid, initial = circle_sim(speed)
+ try:
+ cl = sim.field.curved
+ frozen = {name: getattr(cl, name).copy() for name in GEOMETRY}
+ flags = sim.field.flag.copy()
+ body_state = vars(sim.bodies.get(bid).state).copy()
+ physical = initial.astype(np.float64)
+ flat = physical.reshape(9, -1)
+ cell, direction = cl.fluid_idx.astype(int), cl.dir.astype(int)
+ ex, ey, q = CX[direction], CY[direction], cl.q.astype(float)
+ # Independently project the prescribed arc at true intersections.
+ angle = np.arctan2(cl.ry.astype(float), cl.rx.astype(float))
+ delta = (angle-.7+np.pi) % (2*np.pi)-np.pi
+ aperture = np.abs(delta) < .85
+ radial = np.column_stack((cl.rx, cl.ry)).astype(float)
+ normal = radial / np.linalg.norm(radial, axis=1)[:, None]
+ velocity = speed*np.maximum(0, 1-(delta/.85)**2)[:, None]**2*normal
+ np.testing.assert_array_equal(cl.aperture, aperture)
+ np.testing.assert_allclose(cl.surface_velocity, velocity, atol=1e-8, rtol=0)
+ assert aperture.any() and (~aperture).any()
+ assert (q < .5).any() and (q >= .5).any()
+ a = flat[direction, cell]
+ b = np.empty(len(cell))
+ passive = np.empty(len(cell))
+ for i, (k, d) in enumerate(zip(cell, direction)):
+ rho = flat[:, k].sum()
+ dot = CX[d]*velocity[i, 0] + CY[d]*velocity[i, 1]
+ if cl.fallback_class[i]:
+ passive[i] = a[i]
+ factor = 6
+ elif q[i] < .5:
+ passive[i] = 2*q[i]*a[i] + (1-2*q[i])*flat[d, k-CX[d]-32*CY[d]]
+ factor = 6
+ else:
+ passive[i] = a[i]/(2*q[i]) + (1-1/(2*q[i]))*flat[OPPOSITE[d], k]
+ factor = 3/q[i]
+ b[i] = passive[i]-factor*WEIGHTS[d]*rho*dot
+ force = np.column_stack((ex, ey))*(a+b)[:, None]
+ torque = cl.rx*force[:, 1]-cl.ry*force[:, 0]
+ wanted = np.column_stack((force, torque))
+ # Fluid-node and intersection lever arms differ by q*c: same link torque.
+ fluid_r = radial-q[:, None]*np.column_stack((ex, ey))
+ np.testing.assert_allclose(torque, fluid_r[:, 0]*force[:, 1]-fluid_r[:, 1]*force[:, 0],
+ atol=2e-15, rtol=0)
+ with sim.ctx.activate():
+ sim.stepper.curved_fn(sim.field.ddf_gpu,
+ *(getattr(cl, name+"_gpu") for name in GEOMETRY[:3]),
+ cl.rx_gpu, cl.ry_gpu, cl.rz_gpu, cl.fallback_class_gpu, cl.body_id_gpu,
+ np.uintp(0), np.uintp(0), np.uint32(cl.count),
+ cl.surface_velocity_gpu, cl.link_exchange_gpu,
+ block=(64, 1, 1), grid=((cl.count+63)//64, 1))
+ staged = np.empty_like(initial)
+ cuda.memcpy_dtoh(staged, sim.field.ddf_gpu)
+ incoming = staged.reshape(9, -1)[OPPOSITE[direction], cell+ex+32*ey].astype(float)
+ np.testing.assert_allclose(incoming, b, atol=8e-8, rtol=0)
+ # Incoming populations alone respond with the expected signed normal drive.
+ assert np.all((-np.sum(normal[aperture]*np.column_stack((ex, ey))[aperture], axis=1)) > 0)
+ core = aperture & (np.abs(delta) < .5*.85)
+ assert core.any()
+ assert np.all(speed*(incoming-passive)[core] > 0)
+
+ # Independent pull ledger, including all planar-wall links. Use actual
+ # rounded incoming values to isolate bookkeeping from CUDA arithmetic.
+ arrived = np.stack([np.roll(physical[d], (CY[d], CX[d]), (0, 1)) for d in range(9)])
+ arrived.reshape(9, -1)[OPPOSITE[direction], cell] = incoming
+ mass_exchange = float(np.sum(incoming-a))
+ momentum_exchange = -np.sum(np.column_stack((ex, ey))*(a+incoming)[:, None], axis=0)
+ for y, outgoing in ((1, (4, 6, 7)), (22, (3, 5, 8))):
+ for d in outgoing:
+ values = physical[d, y, :]
+ arrived[OPPOSITE[d], y, :] = values
+ momentum_exchange -= 2*values.sum()*np.array([CX[d], CY[d]])
+ fluid = ((flags & FLUID) != 0) & ((flags & SOLID) == 0)
+ change = (arrived-physical).reshape(9, -1)[:, fluid].sum(axis=1)
+ np.testing.assert_allclose([change.sum(), change@CX, change@CY],
+ [mass_exchange, *momentum_exchange], atol=3e-13, rtol=0)
+ # Start public run from the original state, not a previously staged field.
+ with sim.ctx.activate():
+ cuda.memcpy_htod(sim.field.ddf_gpu, initial)
+ sim.run(1)
+ actual = spatial_state(sim).astype(float).reshape(9, -1)[:, fluid]
+ expected = arrived.reshape(9, -1)[:, fluid]
+ moments = np.stack((np.ones(9), CX, CY))
+ np.testing.assert_allclose(moments@actual, moments@expected, atol=5e-7, rtol=0)
+ # FP32 collision conserves moments up to rounding, not bitwise globally.
+ change = (actual-physical.reshape(9, -1)[:, fluid]).sum(axis=1)
+ budget = 8*np.finfo(np.float32).eps*np.abs(expected).sum()
+ np.testing.assert_allclose([change.sum(), change@CX, change@CY],
+ [mass_exchange, *momentum_exchange], atol=budget, rtol=0)
+ diagnostic = sim.read_surface_exchange(bid)
+ for name, mask in (("wall", ~aperture), ("aperture", aperture)):
+ np.testing.assert_allclose(diagnostic[name], wanted[mask].sum(axis=0),
+ atol=4e-6, rtol=0)
+ total = diagnostic["wall"]+diagnostic["aperture"]
+ np.testing.assert_allclose(sim.read_force(bid, normalize=False), total[:2], atol=4e-6, rtol=0)
+ np.testing.assert_allclose(sim.read_torque(bid, normalize=False), total[2:], atol=4e-6, rtol=0)
+ for name, value in frozen.items():
+ np.testing.assert_array_equal(getattr(cl, name), value)
+ np.testing.assert_array_equal(sim.field.flag, flags)
+ assert vars(sim.bodies.get(bid).state) == body_state
+ finally:
+ sim.close()
+
+
+def test_zero_surface_drive_preserves_legacy_populations_and_geometry():
+ records = []
+ for speed in (None, 0.0):
+ sim, bid, _ = circle_sim(speed)
+ try:
+ sim.run(2)
+ records.append((spatial_state(sim), sim.field.flag.copy(),
+ {name: getattr(sim.field.curved, name).copy() for name in GEOMETRY}))
+ finally:
+ sim.close()
+ np.testing.assert_array_equal(records[0][0], records[1][0])
+ np.testing.assert_array_equal(records[0][1], records[1][1])
+ for name in GEOMETRY:
+ np.testing.assert_array_equal(records[0][2][name], records[1][2][name])
+
+
+def test_surface_checkpoint_continues_and_resets_latest_exchange(tmp_path):
+ sim, bid, _ = circle_sim(.01)
+ try:
+ sim.run(1)
+ saved = spatial_state(sim)
+ count = sim.stepper.step_count
+ spec = asdict(sim.bodies.get(bid).surface_actuation)
+ action = sim.bodies.action.copy()
+ path = sim.save_checkpoint(str(tmp_path / "surface.h5"))
+ with h5py.File(path, "r") as checkpoint:
+ assert checkpoint.attrs["ckpt_version"] == 4
+ geometry = json.loads(checkpoint.attrs["geometry_json"])
+ assert len(geometry[0]["projection_sha256"]) == 64
+ assert json.loads(checkpoint.attrs["surface_actuation_json"]) == [spec]
+ applied = checkpoint["applied_action"][:]
+ sim.run(2)
+ expected = spatial_state(sim)
+ expected_exchange = sim.read_surface_exchange(bid)
+ assert np.linalg.norm(expected_exchange["aperture"]) > 0
+ sim.load_checkpoint(path)
+ assert sim.stepper.step_count == count
+ np.testing.assert_array_equal(spatial_state(sim), saved)
+ np.testing.assert_array_equal(sim.bodies.action, action)
+ actual_action = np.empty_like(applied)
+ with sim.ctx.activate():
+ cuda.memcpy_dtoh(actual_action, sim.bodies.action_gpu)
+ np.testing.assert_array_equal(actual_action, applied)
+ assert asdict(sim.bodies.get(bid).surface_actuation) == spec
+ for values in sim.read_surface_exchange(bid).values():
+ np.testing.assert_array_equal(values, np.zeros(3))
+ restored = sim.save_checkpoint(str(tmp_path / "restored.h5"))
+ with h5py.File(restored, "r") as checkpoint:
+ assert json.loads(checkpoint.attrs["geometry_json"]) == geometry
+ assert json.loads(checkpoint.attrs["surface_actuation_json"]) == [spec]
+ sim.run(2)
+ np.testing.assert_array_equal(spatial_state(sim), expected)
+ for name, values in sim.read_surface_exchange(bid).items():
+ np.testing.assert_array_equal(values, expected_exchange[name])
+ finally:
+ sim.close()
--- /dev/null
+++ b/tests/unit/test_surface_actuation.py
@@ -0,0 +1,230 @@
+"""Static circle apertures: real geometry and host/device-memory protocol tests."""
+from dataclasses import FrozenInstanceError, asdict
+import json
+import math
+from types import SimpleNamespace
+
+import h5py
+import numpy as np
+import pytest
+
+from CelerisLab.body.geometry.circle import CircleGeometry
+from CelerisLab.body.manager import ObjectManager
+from CelerisLab.body.objects import SimObject, SurfaceActuation
+from CelerisLab.common.checkpoint import load_checkpoint, save_checkpoint
+from CelerisLab.config import LBMConfig
+from CelerisLab.lbm.curved_links import CurvedLinkSoA
+from CelerisLab.lbm.descriptors import FLUID, SOLID
+from CelerisLab.simulation import Simulation
+from .test_state_regressions import device, make_field, manager
+
+
+def circle(spec=None, *, center=(6., 5.), radius=2.999, **kwargs):
+ return SimObject(0, geometry=CircleGeometry(*center, radius), center=center,
+ radius=radius, surface_actuation=spec, **kwargs)
+
+
+def projected(obj):
+ cfg = LBMConfig(nx=12, ny=10, viscosity=.1)
+ cfg.validate()
+ flags = ObjectManager.build_flags_for(
+ [obj], np.full(cfg.nx*cfg.ny, FLUID, np.uint16), nx=cfg.nx, ny=cfg.ny, nz=1)
+ cl = CurvedLinkSoA()
+ cl.assign_host(*obj.get_curved_list(cfg.nx, cfg.ny, cfg.nq, domain_flags=flags))
+ cl.assign_surface([obj], cfg, flags)
+ return cl, cfg, flags
+
+
+@pytest.mark.parametrize('speed', [.002, -.002])
+def test_aperture_uses_true_hit_not_fluid_node_or_clamped_q(speed):
+ # Fluid (9,5), direction (-1,+1), circle radius 2.999. Actual t < .01,
+ # whereas the existing Bouzidi q is clipped to .01. The narrow opening
+ # contains the true hit but neither the fluid-node angle nor clipped hit.
+ t = (6-math.sqrt(36-8*(9-2.999**2)))/4
+ angle = math.atan2(t, 3-t)
+ spec = SurfaceActuation(angle, .0001, speed)
+ cl, _, _ = projected(circle(spec))
+ index, = np.flatnonzero((cl.fluid_idx == 9+5*12) & (cl.dir == 8))
+ assert cl.q[index] == np.float32(.01)
+ assert angle > spec.half_width
+ assert math.atan2(.01, 2.99)-angle > spec.half_width
+ assert cl.aperture[index]
+ expected = speed*np.array([3-t, t])/2.999
+ np.testing.assert_allclose(cl.surface_velocity[index], expected, rtol=2e-6, atol=1e-10)
+ np.testing.assert_allclose(np.hypot(cl.rx, cl.ry), 2.999, rtol=0, atol=3e-7)
+ np.testing.assert_array_equal(cl.surface_velocity[~cl.aperture], 0.)
+
+
+def test_surface_zero_preserves_eight_column_geometry_and_mask():
+ plain, zero = circle(), circle(SurfaceActuation(0., .7, 0.))
+ old, new = plain.get_curved_list(12, 10, 9), zero.get_curved_list(12, 10, 9)
+ assert len(old) == len(new) == 8
+ for before, after in zip(old, new):
+ np.testing.assert_array_equal(before, after)
+ np.testing.assert_array_equal(plain.get_flag_mask(12, 10), zero.get_flag_mask(12, 10))
+ cl, _, _ = projected(zero)
+ assert cl.aperture.any()
+ np.testing.assert_array_equal(cl.surface_velocity, 0.)
+ assert projected(plain)[0].surface_velocity is None
+
+
+def test_surface_spec_is_immutable_and_wraps_angles():
+ spec = SurfaceActuation(math.pi, .2, .002)
+ with pytest.raises(FrozenInstanceError):
+ spec.normal_velocity = .003
+ angles = np.array([math.pi-.1, -math.pi+.1, 0.])
+ velocity, mask = spec.project(np.cos(angles), np.sin(angles))
+ np.testing.assert_array_equal(mask, [True, True, False])
+ np.testing.assert_allclose(np.linalg.norm(velocity[:2], axis=1), .002*(1-.5**2)**2,
+ rtol=1e-6)
+ np.testing.assert_array_equal(velocity[2], 0.)
+
+
+@pytest.mark.parametrize('values', [(0., 0., .002), (0., math.pi+.1, .002),
+ (0., .2, 1.), (math.nan, .2, .002),
+ (0., math.inf, .002), (0., .2, math.nan)])
+def test_invalid_surface_spec(values):
+ with pytest.raises(ValueError):
+ SurfaceActuation(*values)
+
+
+@pytest.mark.parametrize('override', [dict(streaming='esopull'), dict(store_precision='FP16S'),
+ dict(ddf_shifting=True), dict(lattice_model='D3Q19')])
+def test_surface_rejects_unsupported_mode_before_registration(override):
+ cfg = LBMConfig(nx=12, ny=10, **override)
+ cfg.validate()
+ bodies = ObjectManager(cfg.nx, cfg.ny, cfg.nz, cfg.nq, cfg)
+ with pytest.raises(ValueError, match='Surface actuation requires'):
+ bodies.add(circle(SurfaceActuation(0., .7, .002)))
+ assert bodies.count == 0
+
+
+@pytest.mark.parametrize('role', ['is_sensor', 'is_force_region'])
+def test_surface_rejects_nonwall_roles(role):
+ with pytest.raises(ValueError):
+ circle(SurfaceActuation(0., .7, .002), **{role: True}).validate_surface_actuation(LBMConfig())
+
+
+@pytest.mark.parametrize('key', ['omega', 'vx', 'vy'])
+def test_surface_set_body_rejects_motion_before_mutation(key):
+ cfg = LBMConfig(nx=12, ny=10)
+ cfg.validate()
+ bodies = ObjectManager(cfg.nx, cfg.ny, cfg.nz, cfg.nq, cfg)
+ item = circle(SurfaceActuation(0., .7, .002))
+ bodies.add(item)
+ sim = Simulation.__new__(Simulation)
+ sim.lbm_cfg, sim.bodies = cfg, bodies
+ state, action = asdict(item.state), bodies.action.copy()
+ with pytest.raises(ValueError, match='stationary'):
+ sim.set_body(0, **{key: .01})
+ assert asdict(item.state) == state
+ np.testing.assert_array_equal(bodies.action, action)
+
+
+def test_surface_projection_rejects_unresolved_aperture_and_bad_writers():
+ with pytest.raises(ValueError, match='no resolved'):
+ projected(circle(SurfaceActuation(.123456, 1e-9, .002)))
+ item = circle(SurfaceActuation(0., .7, .002))
+ cl, cfg, flags = projected(item)
+ flags[cl.fluid_idx[0]] = SOLID
+ with pytest.raises(ValueError, match='active fluid'):
+ cl.assign_surface([item], cfg, flags)
+ flags[cl.fluid_idx[0]] = FLUID
+ cl.fluid_idx[1], cl.dir[1] = cl.fluid_idx[0], cl.dir[0]
+ with pytest.raises(ValueError, match='unique'):
+ cl.assign_surface([item], cfg, flags)
+
+
+def test_surface_public_dict_entry_and_sync_copy():
+ cfg = LBMConfig(nx=12, ny=10)
+ cfg.validate()
+ sim = Simulation.__new__(Simulation)
+ sim.lbm_cfg, sim._initialized = cfg, False
+ sim.bodies = ObjectManager(cfg.nx, cfg.ny, cfg.nz, cfg.nq, cfg)
+ spec = dict(angle=0., half_width=.7, normal_velocity=.002)
+ oid = sim.add_body('circle', center=(6., 5.), radius=2., surface_actuation=spec)
+ copied, = sim.bodies.build_next_objects()
+ assert asdict(copied.surface_actuation) == spec
+ assert copied.surface_actuation is sim.bodies.get(oid).surface_actuation
+
+
+def case(device, spec):
+ field = make_field(device)
+ bodies = manager(field, [circle(spec, radius=2.)])
+ bodies.sync_to_gpu(field)
+ return field, bodies
+
+
+def test_surface_gpu_buffer_lifecycle_with_memory_double(device):
+ field, bodies = case(device, SurfaceActuation(0., .7, .002))
+ cl = field.curved
+ velocity, exchange = cl.surface_velocity_gpu, cl.link_exchange_gpu
+ np.testing.assert_array_equal(np.frombuffer(velocity.data, np.float32).reshape(-1, 2),
+ cl.surface_velocity)
+ np.testing.assert_array_equal(np.frombuffer(exchange.data, np.float32), 0.)
+ cl.assign_surface([], field.cfg, field.flag)
+ cl.upload()
+ assert velocity.freed and exchange.freed
+ assert cl.surface_velocity_gpu is None and cl.link_exchange_gpu is None
+ cl.assign_surface(bodies.objects, field.cfg, field.flag)
+ cl.upload()
+ velocity, exchange = cl.surface_velocity_gpu, cl.link_exchange_gpu
+ cl.free()
+ assert velocity.freed and exchange.freed and cl.count == 0
+
+
+@pytest.mark.parametrize('damage', ['missing_spec', 'legacy_missing_spec', 'different_spec',
+ 'omega', 'vx', 'vy', 'applied_omega'])
+def test_surface_checkpoint_rejects_mismatch_before_any_write(device, tmp_path, damage):
+ field, bodies = case(device, SurfaceActuation(0., .7, .002))
+ path = save_checkpoint(field, SimpleNamespace(step_count=3), field.cfg, bodies,
+ str(tmp_path/'surface.h5'))
+ with h5py.File(path, 'a') as hf:
+ assert hf.attrs['ckpt_version'] == 4
+ if damage in ('missing_spec', 'legacy_missing_spec'):
+ del hf.attrs['surface_actuation_json']
+ if damage == 'legacy_missing_spec':
+ hf.attrs['ckpt_version'] = 3
+ elif damage == 'different_spec':
+ records = json.loads(hf.attrs['surface_actuation_json'])
+ records[0]['normal_velocity'] *= -1
+ hf.attrs['surface_actuation_json'] = json.dumps(records)
+ elif damage == 'applied_omega':
+ action = hf['applied_action'][:]
+ action[-1] = .01
+ hf['applied_action'][:] = action
+ else:
+ records = json.loads(hf.attrs['objects_json'])
+ records[0][damage] = .01
+ hf.attrs['objects_json'] = json.dumps(records)
+ writes, state = device.writes, asdict(bodies.get(0).state)
+ raw = [bytes(ptr.data) for ptr in (field.ddf_gpu, field.temp_gpu, field.flag_gpu)]
+ stepper = SimpleNamespace(_step_count=99)
+ with pytest.raises(ValueError):
+ load_checkpoint(path, field, stepper, field.cfg, bodies)
+ assert device.writes == writes and stepper._step_count == 99
+ assert asdict(bodies.get(0).state) == state
+ assert raw == [bytes(ptr.data) for ptr in (field.ddf_gpu, field.temp_gpu, field.flag_gpu)]
+
+
+@pytest.mark.parametrize('enabled', [False, True])
+def test_surface_checkpoint_rebuild_and_legacy_absence(device, tmp_path, enabled):
+ spec = SurfaceActuation(0., .7, .002) if enabled else None
+ field, bodies = case(device, spec)
+ path = save_checkpoint(field, SimpleNamespace(step_count=3), field.cfg, bodies,
+ str(tmp_path/'surface.h5'))
+ with h5py.File(path, 'a') as hf:
+ assert hf.attrs['ckpt_version'] == (4 if enabled else 3)
+ if not enabled:
+ del hf.attrs['surface_actuation_json']
+ restored, restored_bodies = case(device, spec)
+ stepper = SimpleNamespace(_step_count=0)
+ load_checkpoint(path, restored, stepper, restored.cfg, restored_bodies)
+ assert stepper._step_count == 3
+ for source, target in ((field.ddf_gpu, restored.ddf_gpu), (field.temp_gpu, restored.temp_gpu),
+ (field.flag_gpu, restored.flag_gpu)):
+ assert bytes(source.data) == bytes(target.data)
+ if enabled:
+ np.testing.assert_array_equal(field.curved.surface_velocity, restored.curved.surface_velocity)
+ else:
+ assert restored.curved.surface_velocity is None
Reference in New Issue View Git Blame Copy Permalink