Improve TrackLength tally behavior

[gonuke]

Control flow and logic of TrackLength tally leads to behaviors that undermine accuracy and/or confidence, in particular, the negative track length warning.

Study and resolve the internal flow of the TrackLength tally to better understand and avoid such scenarios, and other numerical cases.

Depends on #93

[makeclean]

I have a proposal for this, the current implementation has a few flaws which we have noticed, in order to improve confidence in the routine I suggest the following,

1. fire a ray from the start point and determine the distance to all tet faces that we intersect, by doing so record the handle to the nodes that belong to that tet
2. Determine if the start point is inside a tet or not

• If inside then the score for the first tet is the first element of the distance vector
• not inside then then the first element of the distance vector is the distance to the first intersection

3. Take the vertex data and determine the tet centroid

• If inside a tet then the score is taken from the distance vector for that tet
• If not inside then this corresponds to a non score
• If we scored for that tet add that tracklength score to the array

[makeclean]

Let me know if this isn't clear, I will try and improve it.

[kldunn]

Wouldn't it be better to determine if the start point is in a tet or not first? Not much point getting the distances if we aren't even in a tet to begin with. There are methods in TrackLengthMeshTally that already do this step separately from the compute method, but we want to make sure it can be extended to hex elements as well.

[makeclean]

We of course can be given a track that starts outside of a mesh and then enters it later on. I.e. a track that starts at 0,0,0 but where the face of an element lies in the plane 10,0,0, checking for not being in a tet will not help us. Instead, knowing the number of face crossings and their distances, we can determine if in tet or not.

[kldunn]

True, but that was not quite what I meant. The current implementation looks for the first tet that is intersected along the track before it does any computations. So it starts off in a tet that is in the mesh, if one exists along the path. If it doesn't exist along the path, then the compute method doesn't do anything else for that track as it is not needed.

[makeclean]

This implementation does the same, it collects all the hits along the ray with the track length and returns them, if 0, no computation.

[makeclean]

I should clarify, I have retasked one of the original functions, which retains all information rather than binning everything but the first hit

[kldunn]

Ah ok, I understand what you mean now. How are you determining whether you are in a tet or not?

[makeclean]

So from the function I get back a list of intersection distances, e.g 0, 10,10.5,20,21,31, since the ray has intersected at these points we either are in or not in a tet, so firstly if the length of the array is 0, then there were no intersections, exit.

If the length is greater than 0, then for each step I compute, pos = origin + direction*distance, take the average, i.e. (pos[1]-pos[0])/2 giving the centroid of the interaction, and then testing to see if this point belongs to any tets, if not, its not score and carry on through the list, until we get a match in a tet, then the track length is given by distance[i+1] - distance[i]

[kldunn]

This sounds very similar to how I think the current implementation already works, with the exception of having all the intersection points stored up front instead of calculating them as we go through tets. Although it sounds like you are planning on using point in volume for every tet instead of just the first one that is intersected. Wouldn't that be less efficient than the current implementation?

[makeclean]

It probably is very similar, however hopefully this/a/some new implementation will avoid the problems of negative track lengths. I would initially settle for slower as long as it is correct.

[gonuke]

I think one of the keys here is to avoid decrementing some length variable with individual subtractions.

One clarification: in step 1 - we should only get the tets that are crossed within the distance of the track we are currently trying to score. It's not clear what the best algorithm is for this. No matter how we do it, we'll need to do ray tracing across as many facets as the number of elements we encounter.

[makeclean]

Indeed, the first query gets all facets crossed from the point p along a vector r of distance sub track.

[kldunn]

When I was looking into this issue I narrowed the negative track length issue down to either the calculations of "t" that were incorrect, and/or the next tetrahedron was not being set correctly. In many cases (not all) "t" was returned as zero, which meant that:

track_length = t - last_t < 0

The value of "t" is set by the ray_tri_intersect method in GeomUtil. Another possibility was in the find_next_tet_by_ray_fire method in TrackLengthMeshTally. Perhaps the next tet was not being set correctly.

It would be interesting to see if this approach makes a difference, and is definitely worth trying.

[kldunn]

For clarification purposes I think that "t" is the distance from the start point to the next surface that was intersected by the ray, which would make the above equation basically the same thing that you are calculating with distance[i+1] - distance[i].

[makeclean]

next_tet_by_ray_fire uses, kdtree->intersect trianlges produces for reason unbeknownst to tim tautges an unorderd list of intersections, there is always a 0 in the list, I think it is this function that needs some attention.

[makeclean]

This was taken care of with commits regarding the MeshTally refactor, ending with commit 1b1f3ab