From 209f3deb3ad62660f68f9a08a2c97103b10fedf3 Mon Sep 17 00:00:00 2001 From: Peter Sharpe Date: Sat, 19 Sep 2026 17:34:43 -0400 Subject: [PATCH 1/2] Keep near-wall normals stable under SDF sign roundoff --- CHANGELOG.md | 2 + .../unified_external_aero_recipe/src/sdf.py | 18 +++--- .../tests/test_sdf_transform.py | 61 ++++++++++++++++--- 3 files changed, 65 insertions(+), 16 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 02cb3aeedc..f62d85a9ec 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -62,6 +62,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- Unified external aero recipe: near-wall SDF normals ignore subresolution sign + noise, preventing tiny negative distances from reversing the outward face normal. - Fixes mesh dtype handling: preserves integer-coordinate precision, normalizes connectivity safely, and rejects integer `.to()` casts. Floating/complex casts preserve the source mesh. diff --git a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/sdf.py b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/sdf.py index 78d546b548..aec111dd99 100644 --- a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/sdf.py +++ b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/sdf.py @@ -60,7 +60,8 @@ class ComputeSDFFromBoundary(MeshTransform): on the surface. Points essentially *on* the surface (boundary-layer points at sub-micron wall distances) instead use the oriented normal of the hit face, since at those distances the closest-point direction - is float32 rounding noise. + and SDF sign can be float32 rounding noise. Within this uncertainty + band, the fallback uses the outward face normal. Parameters ---------- @@ -147,27 +148,28 @@ def apply_to_domain(self, domain: DomainMesh) -> DomainMesh: # rounding noise (or exactly zero) and its direction is # meaningless. Substitute the oriented normal of the hit face: # the exact limit of the closest-point direction at the wall. - # Sign-align with the SDF so the rare interior point keeps - # pointing into the body like its neighbors (the SDF treats - # on-surface as outside, so dist == 0 gets the outward normal). # The band is scale-aware: the closest point carries rounding # noise ~ eps * |coordinate|, so an absolute cutoff under-covers # geometry far from the origin. 128 eps (~1.5e-5 per unit - # coordinate) clears that noise floor with a wide margin, and - # widening the band is free because the substitute is exact. + # coordinate) clears that noise floor with a wide margin. # Computed unconditionally and selected with a mask rather than # branching on ``near_surface.any()`` -- that host readback would # stall the prefetch stream. dist = torch.norm(normals, dim=-1) coord_scale = query_points.abs().amax(dim=-1).clamp(min=1.0) - near_surface = dist < (128.0 * torch.finfo(torch.float32).eps * coord_scale) + surface_tolerance = 128.0 * torch.finfo(torch.float32).eps * coord_scale + near_surface = dist < surface_tolerance face_normals = surface.cell_normals.to(query_points.dtype)[hit_faces] # A degenerate (zero-area) hit face has no meaningful normal -- # ``cell_normals`` returns a zero vector for it. Keep the raw # closest-point direction there instead of substituting zeros. face_normal_ok = (face_normals * face_normals).sum(-1) > 0.5 + # A tiny negative SDF at the wall can be rounding noise. Only + # reverse the fallback when the signed distance resolves a point + # inside the body beyond the same uncertainty band. + genuinely_inside = sdf_values < -surface_tolerance oriented = torch.where( - (sdf_values >= 0).unsqueeze(-1), face_normals, -face_normals + genuinely_inside.unsqueeze(-1), -face_normals, face_normals ) normals = torch.where( (near_surface & face_normal_ok).unsqueeze(-1), oriented, normals diff --git a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_sdf_transform.py b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_sdf_transform.py index 3c9d522787..b6a460f765 100644 --- a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_sdf_transform.py +++ b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_sdf_transform.py @@ -27,6 +27,8 @@ from __future__ import annotations +import pytest +import sdf as sdf_module import torch from sdf import ComputeSDFFromBoundary @@ -78,7 +80,19 @@ def _domain_with_interior(interior_points: torch.Tensor) -> DomainMesh: ) -def test_sdf_normals_near_wall_use_face_normal(): +@pytest.mark.parametrize( + "device", + [ + "cpu", + pytest.param( + "cuda:0", + marks=pytest.mark.skipif( + not torch.cuda.is_available(), reason="CUDA not available" + ), + ), + ], +) +def test_sdf_normals_near_wall_use_face_normal(device): """Near-wall normals equal the oriented face normal, not tangent junk.""" torch.manual_seed(0) n_query = 500 @@ -88,7 +102,7 @@ def test_sdf_normals_near_wall_use_face_normal(): q_xy = torch.rand(n_query, 2) * torch.tensor([0.9, 0.8]) + torch.tensor([9.0, -0.4]) wall_dist = 10 ** (torch.rand(n_query) * 5 - 8) interior = torch.stack([q_xy[:, 0], q_xy[:, 1], 1.0 + wall_dist], -1).float() - domain = _domain_with_interior(interior) + domain = _domain_with_interior(interior).to(device) transform = ComputeSDFFromBoundary( boundary_name="stl_geometry", @@ -103,29 +117,60 @@ def test_sdf_normals_near_wall_use_face_normal(): assert sdf.shape == (n_query, 1) assert normals.shape == (n_query, 3) # All points are outside (or, for the sub-float32-resolution wall - # distances, exactly on) the closed box. - assert torch.all(sdf >= 0) + # distances, exactly on) the closed box. The distance calculation can + # return a slightly negative value from float32 rounding at the wall; + # allow a few ulps at the coordinate scale on either CPU or CUDA. + sign_tolerance = ( + 8.0 + * torch.finfo(torch.float32).eps + * domain.interior.points.abs().amax(dim=-1, keepdim=True).clamp_min(1.0) + ) + assert torch.all(sdf >= -sign_tolerance) # Unit vectors, and the normal must be +z at ALL wall distances. The old # centroid fallback returned the direction away from the box centroid # (z-component ~0.1, i.e. nearly tangent) for points below 1e-6. torch.testing.assert_close( - normals.norm(dim=-1), torch.ones(n_query), atol=1e-4, rtol=0.0 + normals.norm(dim=-1), torch.ones(n_query, device=device), atol=1e-4, rtol=0.0 ) assert torch.all(normals[:, 2] > 0.99) +@pytest.mark.parametrize("offset", [0.0, 1e4]) +def test_sdf_normals_ignore_subresolution_sign_noise(monkeypatch, offset): + """Uncertain SDF signs at the wall cannot reverse its outward normal.""" + domain = _domain_with_interior(torch.tensor([[9.5, 0.0, 1.0]]).repeat(3, 1)) + domain = domain.translate(torch.tensor([offset, 0.0, 0.0])) + noise = torch.tensor([-1.0, 0.0, 1.0]) * torch.finfo(torch.float32).eps + noise *= domain.interior.points.abs().amax() + + def noisy_distance(surface, points, **kwargs): + """Return subresolution distance noise for points on the top face.""" + return noise, points.clone(), torch.full((3,), 2, dtype=torch.int64) + + monkeypatch.setattr(sdf_module, "signed_distance_field", noisy_distance) + result = ComputeSDFFromBoundary().apply_to_domain(domain) + + torch.testing.assert_close( + result.interior.point_data["sdf_normals"], + torch.tensor([[0.0, 0.0, 1.0]]).expand(3, -1), + ) + torch.testing.assert_close(result.interior.point_data["sdf"].squeeze(-1), noise) + + def test_sdf_normals_far_points_keep_closest_point_direction(): """Far from the surface the normals stay the closest-point direction.""" # Off the +x end (closest point on the box rim/edge) and above the middle - # (closest point on the top face). - off = torch.tensor([[12.0, 0.0, 2.0], [5.0, 0.0, 3.0]], dtype=torch.float32) + # (closest point on the top face), plus a clearly inside point below it. + off = torch.tensor( + [[12.0, 0.0, 2.0], [5.0, 0.0, 3.0], [5.0, 0.0, 0.9]], dtype=torch.float32 + ) domain = _domain_with_interior(off) result = ComputeSDFFromBoundary().apply_to_domain(domain) normals = result.interior.point_data["sdf_normals"] - closest = torch.tensor([[10.0, 0.0, 1.0], [5.0, 0.0, 1.0]]) + closest = torch.tensor([[10.0, 0.0, 1.0], [5.0, 0.0, 1.0], [5.0, 0.0, 1.0]]) expected = torch.nn.functional.normalize(off - closest, dim=-1) torch.testing.assert_close(normals, expected, atol=1e-5, rtol=1e-5) From b449804e255affc3ca1d27fda98b4c9370464448 Mon Sep 17 00:00:00 2001 From: Peter Sharpe Date: Sat, 19 Sep 2026 19:08:25 -0400 Subject: [PATCH 2/2] Preserve resolved interior normals in the near-wall fallback --- CHANGELOG.md | 3 +- .../unified_external_aero_recipe/src/sdf.py | 22 +++--- .../tests/test_sdf_transform.py | 79 ++++++++++++++----- 3 files changed, 76 insertions(+), 28 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f62d85a9ec..804bc8574a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -63,7 +63,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed - Unified external aero recipe: near-wall SDF normals ignore subresolution sign - noise, preventing tiny negative distances from reversing the outward face normal. + noise without reversing resolved interior normals, including on translated + geometry. Stored signed distances are unchanged. - Fixes mesh dtype handling: preserves integer-coordinate precision, normalizes connectivity safely, and rejects integer `.to()` casts. Floating/complex casts preserve the source mesh. diff --git a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/sdf.py b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/sdf.py index aec111dd99..1dff608b56 100644 --- a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/sdf.py +++ b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/sdf.py @@ -51,17 +51,18 @@ class ComputeSDFFromBoundary(MeshTransform): Reads the surface mesh from ``domain.boundaries[boundary_name]`` and evaluates the signed distance field at every interior point using :func:`physicsnemo.mesh.spatial.sdf.signed_distance_field`, - a mesh-native, pure-PyTorch implementation backed by a torch BVH. + a mesh-native wrapper around Warp mesh queries. The computed SDF is stored as a scalar field ``(N, 1)`` in ``interior.point_data[sdf_field]``. If ``normals_field`` is set, approximate surface normals ``(N, 3)`` are also stored, computed as - the normalized direction from each query point to its closest point - on the surface. Points essentially *on* the surface (boundary-layer + the normalized direction from the closest surface point to each query + point. Points essentially *on* the surface (boundary-layer points at sub-micron wall distances) instead use the oriented normal of the hit face, since at those distances the closest-point direction - and SDF sign can be float32 rounding noise. Within this uncertainty - band, the fallback uses the outward face normal. + and SDF sign can be float32 rounding noise. A smaller sign-uncertainty + band keeps on-wall normals outward while resolved interior points + retain inward normals throughout the direction-fallback band. Parameters ---------- @@ -164,10 +165,13 @@ def apply_to_domain(self, domain: DomainMesh) -> DomainMesh: # ``cell_normals`` returns a zero vector for it. Keep the raw # closest-point direction there instead of substituting zeros. face_normal_ok = (face_normals * face_normals).sum(-1) > 0.5 - # A tiny negative SDF at the wall can be rounding noise. Only - # reverse the fallback when the signed distance resolves a point - # inside the body beyond the same uncertainty band. - genuinely_inside = sdf_values < -surface_tolerance + # Resolving the side needs only a few ulps, whereas recovering + # a direction needs the wider band above. Reusing that band for + # the sign would force every fallback outward: |sdf| == dist. + # Ignore only small negative distances consistent with float32 + # roundoff, retaining inward normals for resolved inside points. + sign_tolerance = 8.0 * torch.finfo(torch.float32).eps * coord_scale + genuinely_inside = sdf_values < -sign_tolerance oriented = torch.where( genuinely_inside.unsqueeze(-1), -face_normals, face_normals ) diff --git a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_sdf_transform.py b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_sdf_transform.py index b6a460f765..69fc5343c2 100644 --- a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_sdf_transform.py +++ b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_sdf_transform.py @@ -70,6 +70,16 @@ dtype=torch.int64, ) +_DEVICES = [ + "cpu", + pytest.param( + "cuda:0", + marks=pytest.mark.skipif( + not torch.cuda.is_available(), reason="CUDA not available" + ), + ), +] + def _domain_with_interior(interior_points: torch.Tensor) -> DomainMesh: """Wrap query points and the box surface into a DomainMesh.""" @@ -80,18 +90,7 @@ def _domain_with_interior(interior_points: torch.Tensor) -> DomainMesh: ) -@pytest.mark.parametrize( - "device", - [ - "cpu", - pytest.param( - "cuda:0", - marks=pytest.mark.skipif( - not torch.cuda.is_available(), reason="CUDA not available" - ), - ), - ], -) +@pytest.mark.parametrize("device", _DEVICES) def test_sdf_normals_near_wall_use_face_normal(device): """Near-wall normals equal the oriented face normal, not tangent junk.""" torch.manual_seed(0) @@ -137,25 +136,69 @@ def test_sdf_normals_near_wall_use_face_normal(device): @pytest.mark.parametrize("offset", [0.0, 1e4]) -def test_sdf_normals_ignore_subresolution_sign_noise(monkeypatch, offset): +@pytest.mark.parametrize("device", _DEVICES) +def test_sdf_normals_ignore_subresolution_sign_noise(monkeypatch, offset, device): """Uncertain SDF signs at the wall cannot reverse its outward normal.""" domain = _domain_with_interior(torch.tensor([[9.5, 0.0, 1.0]]).repeat(3, 1)) - domain = domain.translate(torch.tensor([offset, 0.0, 0.0])) - noise = torch.tensor([-1.0, 0.0, 1.0]) * torch.finfo(torch.float32).eps + domain = domain.translate(torch.tensor([offset, 0.0, 0.0])).to(device) + noise = ( + torch.tensor([-1.0, 0.0, 1.0], device=device) * torch.finfo(torch.float32).eps + ) noise *= domain.interior.points.abs().amax() def noisy_distance(surface, points, **kwargs): """Return subresolution distance noise for points on the top face.""" - return noise, points.clone(), torch.full((3,), 2, dtype=torch.int64) + # Keep the distance magnitude consistent with the closest-point + # displacement, as the real SDF kernel does. Tangential roundoff on + # a wall can acquire either sign from the winding-number query. + closest = points.clone() + closest[:, 0] += noise.abs() + distances = (points - closest).norm(dim=-1) * noise.sign() + return distances, closest, torch.full((3,), 2, dtype=torch.int64, device=device) monkeypatch.setattr(sdf_module, "signed_distance_field", noisy_distance) result = ComputeSDFFromBoundary().apply_to_domain(domain) torch.testing.assert_close( result.interior.point_data["sdf_normals"], - torch.tensor([[0.0, 0.0, 1.0]]).expand(3, -1), + torch.tensor([[0.0, 0.0, 1.0]], device=device).expand(3, -1), + ) + expected_sdf, _, _ = noisy_distance(None, domain.interior.points) + torch.testing.assert_close( + result.interior.point_data["sdf"].squeeze(-1), expected_sdf + ) + + +@pytest.mark.parametrize("device", _DEVICES) +@pytest.mark.parametrize("use_winding_number", [False, True]) +@pytest.mark.parametrize("offset", [(0.0, 0.0, 0.0), (1e4, 0.0, 0.0), (1e4, 1e4, 1e4)]) +def test_sdf_normals_preserve_resolved_sides_after_translation( + device, use_winding_number, offset +): + """Resolved interior normals stay inward inside the direction-fallback band.""" + # At coordinate scale 1e4, the direction-fallback band is about 0.153. + # Distances 0.01 and 0.1 are inside that band but well above sign noise; + # 0.2 exercises the ordinary closest-point direction. Include the wall + # and exterior points to preserve the original near-wall correction. + points = torch.tensor([[5.0, 0.0, z] for z in [0.8, 0.9, 0.99, 1.0, 1.01, 1.1]]) + domain = _domain_with_interior(points).translate(torch.tensor(offset)).to(device) + transform = ComputeSDFFromBoundary(use_winding_number=use_winding_number) + + result = transform.apply_to_domain(domain) + + expected_normals = torch.tensor( + [[0.0, 0.0, z] for z in [-1.0, -1.0, -1.0, 1.0, 1.0, 1.0]], device=device + ) + torch.testing.assert_close( + result.interior.point_data["sdf_normals"], expected_normals + ) + # Compare against the represented geometry, including float32 rounding + # after translation, and ensure the normal correction never edits SDF. + wall_z = domain.boundaries["stl_geometry"].points[4, 2] + expected_sdf = domain.interior.points[:, 2] - wall_z + torch.testing.assert_close( + result.interior.point_data["sdf"].squeeze(-1), expected_sdf ) - torch.testing.assert_close(result.interior.point_data["sdf"].squeeze(-1), noise) def test_sdf_normals_far_points_keep_closest_point_direction():