Skip to content

Increase numerical robustness of cylindrical and spherical distance to boundary checks #3923

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

Open
GuySten wants to merge 4 commits into
base: develop
Choose a base branch
from
Open

Increase numerical robustness of cylindrical and spherical distance to boundary checks #3923

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
edae136
increase robustness of cylindrical and surface robustness checks
GuySten Apr 15, 2026
d9d37e9
update regressions
GuySten Apr 15, 2026
9fc69e9
add test that was previously failing
GuySten Apr 15, 2026
69ec752
add a test with a small cylinder
GuySten May 4, 2026
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
22 changes: 12 additions & 10 deletions src/surface.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -391,7 +391,7 @@ double axis_aligned_cylinder_evaluate(
{
const double r1 = r.get<i1>() - offset1;
const double r2 = r.get<i2>() - offset2;
return r1 * r1 + r2 * r2 - radius * radius;
return std::sqrt(r1 * r1 + r2 * r2) - radius;
Copy link
Copy Markdown
Contributor

@paulromano paulromano May 1, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I understand the motivation but this changes the meaning of evaluate ("evaluate the surface equation and return the value"), which I'm not a huge fan of even if currently it's only used to check for the sense

}

// The first template parameter indicates which axis the cylinder is aligned to.
Expand All @@ -408,14 +408,15 @@ double axis_aligned_cylinder_distance(Position r, Direction u, bool coincident,
const double r2 = r.get<i2>() - offset1;
const double r3 = r.get<i3>() - offset2;
const double k = r2 * u.get<i2>() + r3 * u.get<i3>();
const double c = r2 * r2 + r3 * r3 - radius * radius;
const double quad = k * k - a * c;
const double R = std::sqrt(r2 * r2 + r3 * r3);
const double quad = k * k - a * (R - radius) * (R + radius);

if (quad < 0.0) {
// No intersection with cylinder.
return INFTY;

} else if (coincident || std::abs(c) < FP_COINCIDENT) {
} else if (coincident ||
std::abs(R - radius) < FP_COINCIDENT * (1.0 + radius)) {
Copy link
Copy Markdown
Contributor

@paulromano paulromano May 1, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Like @gonuke this also makes me uneasy about the tolerance at large radii

Copy link
Copy Markdown
Contributor

@gonuke gonuke May 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that this form is already merged into the similar code for mesh tallies, upon which this PR was based

// Particle is on the cylinder, thus one distance is positive/negative
// and the other is zero. The sign of k determines if we are facing in or
// out.
Expand All @@ -425,7 +426,7 @@ double axis_aligned_cylinder_distance(Position r, Direction u, bool coincident,
return (-k + sqrt(quad)) / a;
}

} else if (c < 0.0) {
} else if (R < radius) {
// Particle is inside the cylinder, thus one distance must be negative
// and one must be positive. The positive distance will be the one with
// negative sign on sqrt(quad).
Expand Down Expand Up @@ -601,7 +602,7 @@ double SurfaceSphere::evaluate(Position r) const
const double x = r.x - x0_;
const double y = r.y - y0_;
const double z = r.z - z0_;
return x * x + y * y + z * z - radius_ * radius_;
return std::sqrt(x * x + y * y + z * z) - radius_;
}

double SurfaceSphere::distance(Position r, Direction u, bool coincident) const
Expand All @@ -610,14 +611,15 @@ double SurfaceSphere::distance(Position r, Direction u, bool coincident) const
const double y = r.y - y0_;
const double z = r.z - z0_;
const double k = x * u.x + y * u.y + z * u.z;
const double c = x * x + y * y + z * z - radius_ * radius_;
const double quad = k * k - c;
const double R = std::sqrt(x * x + y * y + z * z);
const double quad = k * k - (R - radius_) * (R + radius_);

if (quad < 0.0) {
// No intersection with sphere.
return INFTY;

} else if (coincident || std::abs(c) < FP_COINCIDENT) {
} else if (coincident ||
std::abs(R - radius_) < FP_COINCIDENT * (1.0 + radius_)) {
// Particle is on the sphere, thus one distance is positive/negative and
// the other is zero. The sign of k determines if we are facing in or out.
if (k >= 0.0) {
Expand All @@ -626,7 +628,7 @@ double SurfaceSphere::distance(Position r, Direction u, bool coincident) const
return -k + sqrt(quad);
}

} else if (c < 0.0) {
} else if (R < radius_) {
// Particle is inside the sphere, thus one distance must be negative and
// one must be positive. The positive distance will be the one with
// negative sign on sqrt(quad)
Expand Down
Loading
Loading