diff --git a/CHANGELOG.md b/CHANGELOG.md index 02cb3aeedc..804bc8574a 100644 --- a/CHANGELOG.md +++ b/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. 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..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,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 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..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 @@ -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)