crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Improve radial crossing checks in SphericalMesh and CylindricalMesh (#3792)

Browse source 
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
GuySten 2026-02-13 07:10:20 +02:00 • committed by GitHub
parent 4e5deed0cc
commit b145fdd999
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
1 changed files with 10 additions and 8 deletions
Download patch file Download diff file Expand all files Collapse all files

18
src/mesh.cpp
View file

@ -1893,18 +1893,19 @@ double CylindricalMesh::find_r_crossing(
const double inv_denominator = 1.0 / denominator;
const double p = (u.x * r.x + u.y * r.y) * inv_denominator;
double c = r.x * r.x + r.y * r.y - r0 * r0;
double D = p * p - c * inv_denominator;
double R = std::sqrt(r.x * r.x + r.y * r.y);
double D = p * p - (R - r0) * (R + r0) * inv_denominator;
if (D < 0.0)
return INFTY;
D = std::sqrt(D);
// the solution -p - D is always smaller as -p + D : Check this one first
if (std::abs(c) <= RADIAL_MESH_TOL)
// Particle is already on the shell surface; avoid spurious crossing
if (std::abs(R - r0) <= RADIAL_MESH_TOL * (1.0 + std::abs(r0)))
return INFTY;
// Check -p - D first because it is always smaller as -p + D
if (-p - D > l)
return -p - D;
if (-p + D > l)
@ -2178,15 +2179,16 @@ double SphericalMesh::find_r_crossing(
if (r0 == 0.0)
return INFTY;
const double p = r.dot(u);
double c = r.dot(r) - r0 * r0;
double D = p * p - c;
double R = r.norm();
double D = p * p - (R - r0) * (R + r0);
if (std::abs(c) <= RADIAL_MESH_TOL)
// Particle is already on the shell surface; avoid spurious crossing
if (std::abs(R - r0) <= RADIAL_MESH_TOL * (1.0 + std::abs(r0)))
return INFTY;
if (D >= 0.0) {
D = std::sqrt(D);
// the solution -p - D is always smaller as -p + D : Check this one first
// Check -p - D first because it is always smaller as -p + D
if (-p - D > l)
return -p - D;
if (-p + D > l)