Skip to content

s_model_levelset does not subtract centroid_offset, so a moving STL's image points are measured against a phantom body #1896

Description

@sbryngelson

Symptom

For an STL/model immersed boundary whose working centroid has been moved to the measured centre of mass (moving_ibm with geometry 12), every ghost point's image point is placed against a body displaced by 2 * centroid_offset from the one the marker field describes. On our case that was 0.3 chords, enough to put image points on the wrong side of the body and off their own rank.

Cause

The marker test and the levelset disagree about where the body is.

s_apply_ib_patches subtracts centroid_offset after rotating into the body frame. s_airfoil_levelset (src/simulation/m_compute_levelset.fpp:110) does the same:

offset(:) = patch_ib(ib_patch_id)%centroid_offset(:)
...
xy_local = matmul(inverse_rotation, xy_local)

s_model_levelset (m_compute_levelset.fpp:606) does not:

xyz_local = matmul(inverse_rotation, xyz_local)

! 3D models
if (p > 0) then
    call s_distance_normals_3D(...)

so it measures distance to a body sitting centroid_offset away from the one that was voxelised.

Why nothing upstream has seen it

centroid_offset is zero unless the working centroid was moved to the measured centre of mass, which happens only for a moving body of these geometries. Every static STL in the suite has a zero offset, so the two agree and the bug is invisible.

Fix

One line, the same as the airfoil levelset already has:

xyz_local = matmul(inverse_rotation, xyz_local)
xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset

Check that the STL path is otherwise sound

Worth recording because it was the obvious next worry: a 5120-triangle icosphere and the analytic sphere at Re 100 on the same 20-cells-per-diameter grid give C_D = 0.9917 and 0.9913, the same drift and the same 1e-5 transverse forces. The STL path adds nothing to the error.

Part of the moving-IB-across-ranks set filed today; see also #1895.

Found with Claude Code.

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions