Fixed a bug when combining TimeFilter, MeshFilter, and tracklength estimator

[GuySten]

Description

This PR change the particle advance method to stop at time boundaries so the combination of time filter, mesh filter and tracklength estimator works correctly.

This work is inspired by PR #3107 and solves the same issues in a different way.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 15) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[GuySten]

@paulromano, as a reviewer of #3107, are you interested in reviewing this PR as well?

[ilhamv]

Hi @GuySten - Thanks for directing me to this PR!

I agree that this will address the issue in the use of the combination of mesh+time filters and tracklength estimators. I would recommend running a verification test, like using an analytical problem in #3107.

This PR and #3107 effectively give equivalent effect and address the same issue. I think the decision on which one should be merged is really a trade-off that weighs in on the planning of the overall code design.

While resulting in fewer code changes, this PR fixes the tally filtering issue by introducing operations into the main particle transport logic. On the other hand, #3107 isolates the equivalent operations only into the most relevant parts of the code (tally filter), hiding them from the main particle transport logic. (I think, to an extent, this is also part of the reasoning of the creation of a specialized MeshMaterialFilter in #3406.) Keeping the main transport logic as clean/light as possible is useful, as there will be fewer complications when more complex operations inevitably need to be added to the main transport logic. Also, I think the addition of TimedMeshFilter, along with MeshMaterialFilter and others, gives emphasis to the user that individual filters perform their respective tasks individually, and if different behavior is expected from a combination of multiple filters, the appropriate, derived filter should be used instead.

[GuySten]

Thanks @ilhamv for your thoughts about this PR!
I will run the analytical test when I have some free time.

Your comment summarize the main differences between the two aproaches to the same issue.

There are also some other benefits to my approach:

• Support for unstructured mesh out of the box
• From a users' perspective the default behavior of time filter + mesh filter + tracklength is correct
• My implementation can be easily generalized to bank all particles at a time boundary

I also tried to minimize additional overhead when no time boundaries are needed.

I think the ultimate decision should be at the core devs hands.

[ilhamv]

Thanks @GuySten for emphasizing the benefits of the approach. Note that

• While not tested yet, the approach in Add TimedMeshFilter #3107 should also work with other meshes derived from MeshBase;
• Similar to other derived, combination filters like MeshMaterialFilter, we can add warning/recommendation messages whenever the individual corresponding filters are used in combination; and
• Time grid resolutions needed for tallying and time-banking may be different. Separating tally and time-banking time grids is recommended for a flexible, general-purpose MC transport code.

[paulromano]

@ilhamv @GuySten Many thanks to you both for the two contrasting approaches and for the friendly dialogue around the pros and cons of each approach. My gut instinct is that I lean toward the TimedMeshFilter approach for the reasons elucidated by @ilhamv but I'm not yet set on this and I'd like to take a closer look at both PRs before making a final decision.

I'll note that in the case of MeshMaterialFilter there is an additional motivation to minimize the total number of bins. Using a combination of MeshFilter and MaterialFilter would result in the Cartesian product of all mesh and material bins when most combinations are irrelevant for common use cases. Although the same idea may hold here (the Cartesian product of time and mesh bins is very sparse), I don't think there's a way we can exploit that since, unlike the mesh-material case, there is no a priori knowledge on what combinations would get hit (this is the point where @pshriwise chimes in to tell us that we ought to have sparse tallies 😄).

[GuySten]

@ilhamv and @paulromano, after reading your concerns about touching the transport logic, I added another commit that restructured the code in a way that only the tallying section is affected and the main transport logic remain unchanged.

Under the hood the particle is moved back and forth in the tallying routines but in the end it returns to its correct place.

I think this change address the main concerns you had so I am happy to hear your thoughts on the new code structure.

[ilhamv]

Thanks, @GuySten, for continuing to look for alternatives for this!

• I really like the new approach. I think adding the complexity into the score_tracklength_tally is effective as it hits right at the source of the issue: tracklength estimator.
• This approach would address any possible similar future issues coming from a combination of tracklength estimator and time filter.
• The careful book-keeping of relevant time filter grids is a key piece in this approach. Consider sorting the time grid when new time grids are merged in. Imagine the case when one wants to have several overlapping mesh tallies with different time resolution.
• Some particle attributes are modified (and recovered in the end) during the tally scoring. It may not be the most ideal, but since the function is well isolated, it should be safe. (Some of the attribute changes may not be necessary, though.)
• I'm not sure having two functions named score_tracklength_tally (despite the different signatures) is OK. Consider renaming it.

The only downside that I can see relative to the TimedMeshFilter is that in this approach, the particle needs to "stop-and-go" at each time point in the collected tally time grid, even if the time point is not relevant to a tracklength tally.

Howver, I now noticed a downside of the TimedMeshFilter relative to this approach: the new filter TimedMeshFilter is oly useful for tracklength tallies. That's why dealing with the issue at score_tracklength_tally, like this new approach, may be the most optimal.

So, @paulromano, I think I like this approach better. If you agree, I'll go ahead and close #3107.

Nevertheless, we should still verify this approach against a test problem (like the AZURV1 in #3107). Feel free to let me know if you need help with that, @GuySten.

[GuySten]

@ilhamv, Thank you for your thoughts about the new approach!

Regarding the time grid, the add_to_time_grid function already support sorting and removing repeating items from the time grid.

I tried to run your AZURV1 problem according to your verification script but your script does not generate the reference solution. Maybe you can include the scripts needed to generate your nice time dependent plots or at least the script to generate the reference solution.

[ilhamv]

I found the script to generate the reference solution:
Archive.zip
However, I couldn't find the plotting script yet (it has been a while...)

[ilhamv]

Found it!
Archive_plot.zip
plot.py is to generate the animation, and process.py is to generate the error convergence.

[GuySten]

Thanks @ilhamv for the scripts!

I have run the script, and it looks great:

collision estimate:

tracklength estimate:

[GuySten]

@paulromano, could you take some time to review the new approach?

[paulromano]

Many thanks to both @GuySten and @ilhamv for the work on this issue! This solution is indeed quite elegant and I'm happy to merge it. I went through and made a few stylistic updates. I also updated the test to use a chi-squared test instead of comparing values within 3σ (which would result in more random test failures). @GuySten if you're OK with the changes, please feel free to hit the auto-merge button.