Skip to content

Commit

Permalink
Fix SimplePore distance function (#5016)
Browse files Browse the repository at this point in the history 
The `SimplePore` shape used trigonometric functions incorrectly and generated NaN values in the core that would not propagate, but would still be a source of undefined behavior in if/else statements and would cause fatal errors when floating-point exceptions trapping was enabled.

The new distance calculation doesn't rely on trigonometric functions, is 5% faster, and is fully tested.
  • Loading branch information
@kodiakhq
kodiakhq[bot] authored Dec 6, 2024
2 parents f6d12a5 + 310d993 commit 6dbf05f
Show file tree
Hide file tree
Showing 3 changed files with 214 additions and 34 deletions.
5 changes: 2 additions & 3 deletions .github/actions/build_and_check/action.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,12 +5,11 @@ runs:
steps:
- run: |
brew install boost boost-mpi fftw
brew install hdf5-mpi
pip3 install -c requirements.txt "cython<3.0" numpy scipy h5py packaging
pip3 install -c requirements.txt "cython<3.0" numpy scipy packaging
shell: bash
if: runner.os == 'macOS'
- run: |
export myconfig=maxset with_cuda=false with_gsl=false test_timeout=800 check_skip_long=true
export myconfig=maxset with_cuda=false with_gsl=false with_hdf5=false test_timeout=800 check_skip_long=true
bash maintainer/CI/build_cmake.sh
shell: bash
# This is a workaround for the unfortunate interaction of MacOS and OpenMPI 4
Expand Down
40 changes: 22 additions & 18 deletions src/shapes/src/SimplePore.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include <shapes/SimplePore.hpp>

#include <utils/Vector.hpp>
#include <utils/math/sqr.hpp>

#include <cassert>
#include <cmath>
Expand All @@ -39,28 +40,33 @@ std::pair<double, double> SimplePore::dist_half_pore(double r, double z) const {
assert(r >= 0.0);

/*
* We have to find the line that splits area 1 (r determines distance) from
* area 2 (z determines distance) inside pore. In area 3 we have to consider
* z and r to determine the distance.
* We have to find the line that splits area (1), where r determines
* the distance, from area (2), where z determines the distance.
* In area (3) we have to consider z and r to determine the distance.
* The line that separates area (1) from area (2) has parametric equation
* @f$ r = c_z + c_r - z @f$.
*
* | x
* | 2 x
* | x
* _|_ x 1
* \| ^ r
* 3 |-------|
* | z <-
* | . ^ r
* | (2) . |
* | . |
* ........|.. (1) |
* ^ \ : |
* | \_________|
* c_r | (3) : |
* | :<------>|
* v : c_z |
* z <-----------------+
*/

if ((z <= c_z) && (r <= (c_z + c_r - z))) {
/* Cylinder section, inner */
/* Cylinder section, inner, area (1) */
return {m_rad - r, 0};
}
if (((z >= c_z) && (r >= c_r)) || ((z <= c_z) && (r > (c_z + c_r - z)))) {
/* Wall section and outer cylinder */
/* Wall section and outer cylinder, area (2) */
return {0, m_half_length - z};
}
/* Smoothing area */
/* Smoothing area (3) */
/* Vector to center of torus segment */
auto const dr = c_r - r;
auto const dz = c_z - z;
Expand All @@ -76,7 +82,7 @@ void SimplePore::calculate_dist(const Utils::Vector3d &pos, double &dist,
Utils::Vector3d &vec) const {
/* Coordinate transform to cylinder coords
with origin at m_center. */
Utils::Vector3d const c_dist = pos - m_center;
auto const c_dist = pos - m_center;
auto const z = e_z * c_dist;
auto const r_vec = c_dist - z * e_z;
auto const r = r_vec.norm();
Expand All @@ -97,10 +103,8 @@ void SimplePore::calculate_dist(const Utils::Vector3d &pos, double &dist,
} else {
// smoothing area
if (std::abs(z) >= c_z) {
auto const angle = std::asin((std::abs(z) - c_z) / m_smoothing_rad);
auto const dist_offset =
m_smoothing_rad - (std::cos(angle) * m_smoothing_rad);
if (m_half_length < std::abs(z) || r <= (m_rad + dist_offset)) {
auto const d_sq = Utils::sqr(r - c_r) + Utils::sqr(z - c_z);
if (d_sq > Utils::sqr(m_smoothing_rad)) {
side = 1;
}
}
Expand Down
203 changes: 190 additions & 13 deletions testsuite/python/simple_pore.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,20 +19,24 @@

import unittest as ut
import unittest_decorators as utx
import numpy as np
import espressomd
import espressomd.shapes

# Integration test for simple pore
# The rationale is to hit the pore everywhere with particles
# and check that it does not blow up. The cylinder is needed
# because the pore is tilted with respect to the box, without
# it particles could enter the constraint over the periodic boundaries,
# leading to force jumps.


@utx.skipIfMissingFeatures(["LENNARD_JONES"])
class SimplePoreConstraint(ut.TestCase):

box_yz = 15.
box_x = 20.
system = espressomd.System(box_l=[box_x, box_yz, box_yz])
system.time_step = 0.01
system.cell_system.skin = 0.4

def tearDown(self):
self.system.constraints.clear()
self.system.part.clear()
self.system.non_bonded_inter.reset()

def test_orientation(self):
pore = espressomd.shapes.SimplePore(axis=[1., 0., 0.], radius=2., smoothing_radius=.1,
length=2., center=[5., 5., 5.])
Expand All @@ -43,13 +47,186 @@ def test_orientation(self):
d, _ = pore.calc_distance(position=[5., 5., .0])
self.assertLess(d, 0.)

def test_distance_calculation(self):
"""
Check the distance calculation function for all 3 sections of
a simple pore. The test works in cylindrical coordinates.
"""
twopi = 2. * np.pi
shape = espressomd.shapes.SimplePore(
axis=[0., 0., 1.], radius=3., smoothing_radius=.1, length=6.,
center=self.system.box_l / 2.)

def calc_cartesian_coord(body, pos):
"""
Convert cylindrical coordinates in the body frame
to Cartesian coordinates in the lab frame.
"""
origin = body.center
r, phi, z = pos
x = r * np.cos(phi) + origin[0]
y = r * np.sin(phi) + origin[1]
z = z + origin[2]
return [x, y, z]

def calc_distance(shape, pos):
"""Compute the distance vector from a shape surface."""
return shape.calc_distance(
position=calc_cartesian_coord(shape, pos))

def angledist(val, ref):
"""Compute the absolute distance between two angles."""
twopi = 2. * np.pi
d1 = np.abs((ref - val) % twopi)
d2 = np.abs((ref - val) % twopi - twopi)
return min(d1, d2)

tols = {"atol": 1e-10, "rtol": 1e-7}

# scan points along the cylinder main axis
for h in np.linspace(-1., 1., 21):
z = h * (shape.length / 2. - shape.smoothing_radius)
for phi in np.linspace(0., twopi, 31):
dist, vec = calc_distance(shape, (0., phi, z))
np.testing.assert_allclose(dist, shape.radius, **tols)
np.testing.assert_allclose(vec[0], -shape.radius, **tols)
np.testing.assert_allclose(vec[2], 0., **tols)

# scan cylinder section
for ref_dist in np.linspace(-0.99, 0.99, 21) * shape.smoothing_radius:
r = shape.radius - ref_dist
for h in np.linspace(-0.99, 0.99, 21):
z = h * (shape.length / 2. - shape.smoothing_radius)
for phi in np.linspace(0., twopi, 31):
pos = calc_cartesian_coord(shape, (r, phi, z))
dist, vec = shape.calc_distance(position=pos)
is_inside = shape.is_inside(position=pos)
cur_r = np.linalg.norm(vec[:2])
np.testing.assert_allclose(dist, ref_dist, **tols)
np.testing.assert_allclose(cur_r, np.abs(ref_dist), **tols)
np.testing.assert_allclose(vec[2], 0., **tols)
if np.abs(ref_dist) > 1e-4:
ref_phi = phi if ref_dist < 0. else phi + np.pi
cur_phi = np.arctan2(vec[1], vec[0])
angle_diff = angledist(cur_phi, ref_phi)
np.testing.assert_allclose(angle_diff, 0., **tols)
self.assertEqual(is_inside, ref_dist < 0.)

# scan wall section
for r in np.linspace(2.01, 3.01, 21) + shape.radius:
for ref_dist in np.linspace(-2.0, 2.0, 21):
for sgn in [+1, -1]:
z = ref_dist + shape.length / 2.
for phi in np.linspace(0., twopi, 31):
pos = calc_cartesian_coord(shape, (r, phi, sgn * z))
dist, vec = shape.calc_distance(position=pos)
is_inside = shape.is_inside(position=pos)
ref_z = sgn * ref_dist
np.testing.assert_allclose(dist, ref_dist, **tols)
np.testing.assert_allclose(vec[0], 0., **tols)
np.testing.assert_allclose(vec[1], 0., **tols)
np.testing.assert_allclose(vec[2], ref_z, **tols)
if np.abs(ref_dist) > 1e-4:
self.assertEqual(is_inside, ref_dist < 0.)

# scan torus section
for phi in np.linspace(0., twopi, 21):
r = shape.radius + shape.smoothing_radius
z = shape.length / 2. - shape.smoothing_radius
origin = np.array(calc_cartesian_coord(shape, (r, phi, z)))
for rho in np.linspace(0.01, 2.01, 21) * shape.smoothing_radius:
ref_dist = rho - shape.smoothing_radius
for theta in np.linspace(np.pi / 2. + 0.01, np.pi - 0.01, 31):
dr = rho * np.cos(theta)
dz = rho * np.sin(theta)
rdr = ref_dist * np.cos(theta)
rdz = ref_dist * np.sin(theta)
pos_cyl = (r + dr, phi, z + dz)
normal = (r + rdr, phi, z + rdz)
pos = calc_cartesian_coord(shape, pos_cyl)
dist, vec = shape.calc_distance(position=pos)
is_inside = shape.is_inside(position=pos)
ref_vec = calc_cartesian_coord(shape, normal) - origin
np.testing.assert_allclose(dist, ref_dist, **tols)
np.testing.assert_allclose(np.copy(vec), ref_vec, **tols)
if np.abs(ref_dist) > 1e-4:
self.assertEqual(is_inside, ref_dist < 0.)

@utx.skipIfMissingFeatures(["LENNARD_JONES"])
def test_scattering(self):
"""
Create a line of particles along the pore main axis with a velocity
vector pointing up. Particles inside the pore and outside the pore
will bounce up and down (force vector is perpendicular to main axis).
Particles at the pore mouth will scatter along the main axis due
to reflection by the torus. Two peaks will appear in the trajectory
of the particles x-position, which roughly follow a sine wave
(we'll use a polynomial fit). Since we don't use a thermostat,
the two peaks are perfectly symmetric.
"""
system = self.system
system.non_bonded_inter[0, 1].lennard_jones.set_params(
epsilon=1., sigma=1., cutoff=2**(1. / 6.), shift="auto")
xpos_init = np.arange(800) * (self.box_x / 800.)

def get_diffusion():
return np.copy(self.system.part.all().pos[:, 0]) - xpos_init

def get_series():
ydata = get_diffusion()
xdata = np.arange(250, 340)
xdata = xdata[np.nonzero(ydata[250:340])[0]]
ydata = ydata[xdata[0]:xdata[-1] + 1]
xdata = np.array(xdata, dtype=float) / system.box_l[0]
return (xdata, ydata)

system.constraints.add(
particle_type=0, penetrable=False, only_positive=False,
shape=espressomd.shapes.SimplePore(
axis=[1., 0., 0.], radius=3., smoothing_radius=.1, length=5.,
center=system.box_l / 2.))

for x0 in xpos_init:
rpos = [x0, 0.5 * self.box_yz, 0.5 * self.box_yz]
system.part.add(pos=rpos, type=1, v=[0., 1., 0.])

# integrate until particles interact with the simple pore
system.integrator.run(300)
system.time = 0.
xdata, ydata = get_series()
coefs_init = np.polyfit(xdata, ydata, 4)[::-1]
coefs_growth = np.array([5398.954, -1483.69746, 152.9315,
-7.0115546, 0.1207058])

# verify particles trajectory
for _ in range(10):
system.integrator.run(20)
# check symmetry
ydata = get_diffusion()[1:]
np.testing.assert_allclose(ydata, -ydata[::-1], rtol=0., atol=1e-9)
# check one peak
coefs = coefs_init + system.time * coefs_growth
xdata, ydata = get_series()
ydata_hat = np.sum([coefs[j] * xdata**j for j in range(5)], axis=0)
magnitude = -np.min(ydata_hat)
deviation = (ydata - ydata_hat)[2:-2] # remove outliers
rmsd = np.sqrt(np.mean(deviation**2))
self.assertLess(rmsd / magnitude, 0.01)

@utx.skipIfMissingFeatures(["LENNARD_JONES"])
def test_stability(self):
box_yz = 15.
box_x = 20.
system = espressomd.System(box_l=[box_x, box_yz, box_yz])
system.time_step = 0.01
system.cell_system.skin = 0.4
"""
Stability test for a simple pore tilted along the box diagonal.
The rationale is to hit the pore everywhere with particles
and check that no singularity occurs. The cylinder is needed
because the pore is tilted with respect to the box, without
it particles could enter the constraint over the periodic
boundaries, leading to force jumps.
"""
box_yz = self.box_yz
box_x = self.box_x
system = self.system
lj_eps = 1.0
lj_sig = 1.0
lj_cut = lj_sig * 2**(1. / 6.)
Expand Down

0 comments on commit 6dbf05f

Please sign in to comment.