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.
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
Initial portable implementation of ORANGE: oak ridge geometry engine #291
Conversation
The most important question — what is the "AN" in "ORANGE"?! |
Bwaha this was a discussion topic for a while! I want it to be "Arbitrarily Named" or "Acronym-Neutral", but the vote was to call it "Adaptable Nested"... |
Type-deleted collections-based surfaces are working! Next is to add a GPU test :) and then I think this phase will be ready to merge (unless I should keep going on this branch...) |
Pre-merge question probably for @tmdelellis : how should we namespace this package? SCALE's version uses an |
1de0101
to
34093d7
Compare
@tmdelellis et al: I think this is ready for review as the first phase of geometry implementation. I'll have another PR for transporting among volumes, then another one for Celeritas integration, then another for multiple universes etc. |
Before: ``` /home/s3j/.local/src/celeritas/src/orange/surfaces/detail/SurfaceAction.hh(112): warning: calling a __host__ function("__builtin_unreachable") from a __host__ __device__ function("celeritas::detail::SurfaceAction< ::celeritas_test::CalcSenseDistance> ::operator ()") is not allowed ``` Using `__trap` didn't reduce the number of instructions, but simply omitting the warning works.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Awesome @sethrj! Looks great from my end.
This is a nonstandard extension that lets us more easily interrogate array sizes.
CI failure is due to disk space issues because ... of course ... vecgeom. |
Reviewing now! |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Looks good @sethrj
*/ | ||
CELER_CONSTEXPR_FUNCTION real_type no_intersection() | ||
{ | ||
return numeric_limits<real_type>::infinity(); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I may be misreading this, but this will not be able to be called on device.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sneaky here: numeric_limits
is the one in the celeritas namespace. You imply a good point which is that maybe we should change such class names to "resemble" the standard library without actually using the same names: i.e. class NumericLimits
rather than class numeric_limits
so there's less ambiguity. That can be deferred to another merge but I'm curious of your and the group's thoughts. We already did this with Span and Array, for example, so I'd say we should probably do the same here.
const real_type x = pos[0]; | ||
const real_type y = pos[1]; | ||
const real_type z = pos[2]; | ||
|
||
Real3 norm; | ||
norm[0] = 2 * a_ * x + d_ * y + f_ * z + g_; | ||
norm[1] = 2 * b_ * y + d_ * x + e_ * z + h_; | ||
norm[2] = 2 * c_ * z + e_ * y + f_ * x + i_; | ||
|
||
normalize_direction(&norm); | ||
return norm; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hi @sethrj, while looking around in the ORANGE code, I came across this method and how you use it to calculate the safety, ie the smallest distance from a point to the surface. Do you have (public) notes how this is derived mathematically?
I'm wondering about the cases of ellipsoids (or the special case of spheres), cones and cylinders where there are only quadratic terms. If I read the code correctly, it would then produce a null vector for a point at the origin. This likely will give no meaningful answer for the safety, which is trying to intersect a line following the normal vector, will it?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@hahnjo The normal vector at a point on a surface is just the gradient of the surface equation evaluated at that point (x,y,z). I'm not sure what meaning it holds if you evaluate it for a point that's not on the surface (as it appears to be used for the safety distance) so maybe @sethrj can enlighten us both 😄
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes thank you Paul. The code highlighted here is indeed the gradient of the quadratic expression q(x,y,z)
, which is the vector normal on the surface but also (modulo the sign) the steepest-descent search direction to solve for a point on the surface q(x,y,z)=0
. There are of course points (for ellipsoids), lines (for elliptic cylinders), and I think planes (for hyperbola/parabola/cone) where multiple points are equally close and so components of the gradient go to zero -- so the CalcSafetyDistance
algorithm is incomplete for certain degenerate regions.
Anyway, once the steepest descent direction is calculated, the CalcSafetyDistance
algorithm solves the quadratic equation that results from evaluating the quadric equation along that direction, resulting in the analytic closest point along that vector.
The comment "I need to check the math on this" on the functor is because this is all rather hand-wavy and I haven't formally evaluated that the nearest surface point is along that gradient... 😬 Nor have I worked out how to calculate the distance for degenerate points. I haven't made safety distance a priority, so this functor is not yet used anywhere: I meant it primarily as an example of how to operate generically on a surface.
One less obvious important point about a complete safety distance algorithm in ORANGE is that the minimum distance to any surface is only the correct safety in the case of a "simple" volume, where crossing any surface exits the volume. In the case of a volume with internal surfaces (e.g. a sphere with a box excluded) the distance to the nearest surface may be far smaller than the distance to the nearest surface crossing that actually exits the volume.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No, I'm fairly convinced this doesn't hold mathematically. I've played with tubes on GeoGebra and it works for a = b = 1
. For a = 2
you can already see that the safety sphere intrudes into the tube, and for a = 3
the method doesn't find any intersection at all.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Well, you're definitely right! Sorry for the misleading class and incorrect algorithm. Maybe the safety could be approximated with a few iterations with steepest descent?
Thanks for digging into this, and for the link to the very nifty visualization tool!
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I created an issue #352 so we can have further discussions there rather than buried away inside this code review comment 😄
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The most promising literature that I could find to date is ON THE DISTANCE FROM A POINT TO A QUADRIC SURFACE. What we need for safety is a sort of "conservative" estimate, ie the returned value can be too low but it must never be too large, I think.
Thanks for digging into this, and for the link to the very nifty visualization tool!
Hehe, I never liked it in high school, but it's one of the few free tools to easily visualize objects and surfaces in 3 dimensions...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What we need for safety is a sort of "conservative" estimate, ie the returned value can be too low but it must never be too large, I think
Exactly. Anything that can give us a conservative estimate (especially if the accuracy increases near the surface itself) would be good.
This will be the first set of functionality for an in-house geometry (see #267 ). It will include: