Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix SimplePore distance function #5016

Merged
merged 2 commits into from
Dec 6, 2024
Merged

Fix SimplePore distance function #5016

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f83d98d
Fix SimplePore distance function
jngrad Dec 4, 2024
310d993
CI: Drop flaky hdf5-mpi package from macOS
jngrad Dec 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 2 additions & 3 deletions .github/actions/build_and_check/action.yml
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading