Aero recipe: prevent near-wall SDF sign noise from reversing normals

CHANGELOG.md

@@ -62,6 +62,9 @@ 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 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.

examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/sdf.py

@@ -51,16 +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
-    is float32 rounding noise.
+    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
     ----------
@@ -147,27 +149,31 @@ 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
+            # 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(
-                (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

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
 
@@ -68,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."""
@@ -78,7 +90,8 @@ def _domain_with_interior(interior_points: torch.Tensor) -> DomainMesh:
     )
 
 
-def test_sdf_normals_near_wall_use_face_normal():
+@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)
     n_query = 500
@@ -88,7 +101,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 +116,104 @@ 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])
+@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])).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."""
+        # 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]], 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
+    )
+
+
 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)