Skip to content

core/immersed_boundary: fix Tribend concave reference angle (bug-sweep #5)#5385

Draft
RudolfWeeber wants to merge 1 commit into
espressomd:pythonfrom
RudolfWeeber:fix/bug-5-tribend-concave-angle
Draft

core/immersed_boundary: fix Tribend concave reference angle (bug-sweep #5)#5385
RudolfWeeber wants to merge 1 commit into
espressomd:pythonfrom
RudolfWeeber:fix/bug-5-tribend-concave-angle

Conversation

@RudolfWeeber

Copy link
Copy Markdown
Contributor

IBMTribend::calc_forces() builds the dynamic dihedral angle as
theta = acos(n1n2) * copysign(1, desc), confined to [-pi, pi], while
initialize() maps a concave reference dihedral (desc = dx1 . (n1 x n2) < 0)
onto theta0 = 2
pi - acos(...) in [pi, 2pi]. The two angles thus live on
different branches and the raw difference DTh = theta - theta0 was never
reduced modulo 2
pi. At a concave reference equilibrium DTh = -2pi instead
of 0, applying a spurious ~2
pi*kb bending force at the relaxed shape --
exactly the canonical biconcave-RBC / vesicle membrane use case.

Fix: reduce DTh to (-pi, pi] via std::remainder(theta - theta0, 2pi),
keeping the signed theta and the copysign(1, theta) orientation factor
unchanged. This zeroes the force at both convex and concave reference
shapes and yields correctly restoring forces on both sides of equilibrium.
(Note: the alternative of remapping the dynamic theta to [0, 2
pi] would
collapse copysign(1, theta) to +1 and make concave references unstable.)

Adds two regression tests to testsuite/python/ibm.py: a convex control
(passes even on buggy code) and the concave equilibrium case (force ~6 on
buggy HEAD, < 1e-6 after the fix).

Co-Authored-By: Claude Opus 4.8 noreply@anthropic.com

🤖 Generated with Claude Code

…#5)

IBMTribend::calc_forces() builds the dynamic dihedral angle as
theta = acos(n1*n2) * copysign(1, desc), confined to [-pi, pi], while
initialize() maps a concave reference dihedral (desc = dx1 . (n1 x n2) < 0)
onto theta0 = 2*pi - acos(...) in [pi, 2*pi]. The two angles thus live on
different branches and the raw difference DTh = theta - theta0 was never
reduced modulo 2*pi. At a concave reference equilibrium DTh = -2*pi instead
of 0, applying a spurious ~2*pi*kb bending force at the relaxed shape --
exactly the canonical biconcave-RBC / vesicle membrane use case.

Fix: reduce DTh to (-pi, pi] via std::remainder(theta - theta0, 2*pi),
keeping the signed theta and the copysign(1, theta) orientation factor
unchanged. This zeroes the force at both convex and concave reference
shapes and yields correctly restoring forces on both sides of equilibrium.
(Note: the alternative of remapping the dynamic theta to [0, 2*pi] would
collapse copysign(1, theta) to +1 and make concave references unstable.)

Adds two regression tests to testsuite/python/ibm.py: a convex control
(passes even on buggy code) and the concave equilibrium case (force ~6 on
buggy HEAD, < 1e-6 after the fix).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
@RudolfWeeber

Copy link
Copy Markdown
Contributor Author

not obvious to me. Needs more detailed analysis.

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

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant