Restricted file source

[ebknudsen]

Description

This PR adds functionality to spatially restrict emission of particles from a FileSource based on either domains (much like IndependentSource) or on the definition of a hypercube (strictly a hyperrectangle) in spatial, directional, energy, and time dimensions. This can for instance be useful when the particle-file was created by an external tool and some particles fall outside some important boundary in the problem.

Two rejection schemes are implemented:

1. "resample" a new particle is picked at random, until a valid particle has been accepted. N.b. this may potentially lead to biasing.
2. "kill" the particle is accepted but given wgt=0 to flag it for garbage collection.

I have not yet added regression tests.

Fixes # (issue) None defined.

Checklist

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

[gridley]

Very nice. But as for implementation, it looks like this may duplicate the other rejection code introduced in #2235. In fact there is no reason this should strictly pertain to a file source; it applies to any sort of source. This is very nice to have.

Did you consider just hoisting that code into the source base class rather than separately implementing it for the file source class?

Secondly, I'm not sold on the kill strategy being different from just resampling, but with some extra steps in between. Can you expand on how you thought this would bias the results?

Lastly, as for the "hypercube" terminology... IMO it just adds a layer of opacity that obscures what the intention of the code is in favor of using some fun mathematical terminology. It'd be much nicer to just include time_bounds etc. in addition to the spatial bounds included in #2235.

[ebknudsen]

@gridley Thanks for the comments - and yes, there is some code duplication from #2235, i.e. the Independent source. It's not entirely identical but close. I did consider elevating to the source base class, but as it is not completely identical, I deemed it better to leave it like this, to avoid creating a set of special cases that need to be dealt with. I could reconsider that.

The difference btw 'kill' and 'resample' ends up in normalization. If resampling you are always guaranteed that your source will emit the full source strength, albeit in a restricted phase space, which means that phase space is being oversampled. It is of course not all that difficult to downweight all your results by some factor to reflect this, but if you can live with somewhat inefficient runtime, kill is perhaps an easier option. Then the thing behaves as if you were to put some perfect neutron/photon absorber in that phase space. This is for "fixed source" runs of course.

Lastly, agreed - hypercube (or in fact hyperrectangle to be precise) may be an overly fancy term. I have nothing against using something more prosaic. It may be that it is much better to have a set of coordinate specific bounds, as you suggest, or even the option of a boolean function.

I did consider allowing a filter to be passed in as a restriction to piggy-back on existing infrastructure, but rejected that idea, since it requires a connection to the tally-infrastructures. I was a little worried about odd side-effects.

[ebknudsen]

@gridley Upon rereading your comment, perhaps you are suggesting hoisting the entire piece of code into the SourceBase class? I first read it as only the domain related bit, but this may not have been what you meant?

[gridley]

Hey Erik, yes, I think if it doesn't specifically pertain to sources from a file this should be in the base class. The capability added here seems to pertain to any source, so it would be nice to allow it as well for instance with a box source or similar.

Now there is a slight overlap with that other rejection sampling code, but finding their intersection shouldn't be too bad. Do you see any barrier to that?

[ebknudsen]

Hi @gridley I see what you mean, and I agree. I don't see any real obstacles as such, so I'll get to it. Thanks!

[ebknudsen]

I have now elevated the source restriction code to an intermediate class RestrictedSource from which IndependentSource and FileSource now inherit. I am sure t here are more elegant ways of doing this by none occur to me at present.

Interesting point: The default restriction strategy (which remains unchanged) by SourceIndependent is a mixture of "resample" and "kill". In terms of spatial sampling the modus operandi is "kill" whereas for energy it is "resample"

[gridley]

Interesting observation about kill vs. resample! Hm, I guess we really should have either one be possible, so I'll wait for someone else to weigh in on that. I will review the new code now, though!

[ebknudsen]

Interesting observation about kill vs. resample! Hm, I guess we really should have either one be possible, so I'll wait for someone else to weigh in on that. I will review the new code now, though!

Well - I suppose it is kind of uncommon for it to happen in terms of energy - it only checks if the sampled energy is "out of bounds" for cross-sections etc. I'd hazard a guess it is pretty insignificant.

[ebknudsen]

The asymmetric_lattice regression test fails. I don't know yet why or even if it is related to this, but will check.

[ebknudsen]

@gridley So - rereading the code again - I find that in terms of resample_strategy the original behavior is to always "resample" both for spatial and energy (and drop a fatal error if resample freq. >95%). I.e. my earlier claim was wrong. Further, the current implementation of "kill" turns out in our testcases to cause MPI-related problems. At present I think the way forward is to mimic the original behavior, and leave "kill" to future implementation. Would you agree?

[gridley]

Yes I completely agree on that, although the kill on reject feature would indeed be helpful!

[ebknudsen]

... and to turn things around yet again I'm pretty sure that https://github.com/openmc-dev/openmc/pull/2916/commits/aa8628575623b7e7ac6209154afacfbb63b40e15 fixed the kill-problem. Will run some more tests and ping when I'm positive.

[ebknudsen]

@gridley I am now positive this is right - will proceed to setup some kind of unit test following the procedure in test_source.py

[paulromano]

@ebknudsen If it's OK with you, I may do a bit of refactoring on your branch. One thing I wanted to bring up before doing so is the lower_left_ and upper_right_ members. One feature we are planning on adding in the future is "cookie cutter" regions that could be used for domain rejection but that otherwise aren't part of the transport geometry. This would allow a user to define whatever shape for rejection (including a box with a lower-left and upper-right coordinate) but would be more general. In light of that, are you OK with me removing lower_left_ and upper_right_ from this PR?

[ebknudsen]

That is entirely fine by me - and good that you bring it up now - I realized that I needed to refactor the python layer anyway, since it does not reflect the restricted source being in the base class. Is it better I do that now before you do the necessary refactoring for cookie cutter?

[ebknudsen]

out of curiosity: one notion from the beginning was to use a surface/set of surfaces as a restriction geometry, but it would require adding surfaces into the geometry model, no? (Or of course relying on already existing ones) I could not see a nice way of doing that at that point - I trust you can come up with a better idea than I could.

[paulromano]

@ebknudsen I can take care of the Python bit while I'm doing my refactoring. Thanks!

one notion from the beginning was to use a surface/set of surfaces as a restriction geometry, but it would require adding surfaces into the geometry model, no?

Correct, this is the idea that we will likely pursue eventually where you can define a CSG cell as normal that it used as a restriction for the source. As you point out, it means that geometry information would be need to stored somewhere so we'll have to think about that aspect.

[paulromano]

@ebknudsen I just pushed a major update on this branch. Let me know if you're alright with the changes I've made (happy to iterate with you if it doesn't meet your needs). The highlights are:

• From the Python side, everything is now specified by the user as a dictionary of "constraints". So this looks like:

  source = openmc.FileSource(path, constraints={'domains': [steel], 'time_bounds': [0.5, 2.0]})

  That same constraints argument can now be given for any of the source classes.

• I tried to unify the reading, writing, and use of the constraint information as much as possible for both Python and C++.
• Tests have been added for everything other than using constraints with CompiledSource

Since I have now written quite a bit of the code on this branch, we should have an independent reviewer. @gridley would you mind giving this another review?

[paulromano]

One more thought to add — with the generic constraints dictionary, we could also refactor the only_fissionable argument from openmc.stats.Box to be its own constraint, which would then allow a user to apply that constraint for any source distribution.

[ebknudsen]

@paulromano I've taken a quick look - looks good to me. I really like the dictionary idea. I think this does what we need/intend. I'll run a few tests tom. but I don't expect any problems. Totally agree with having a 2nd review btw.

[paulromano]

@gridley Thanks for the detailed comments. I pretty much agree with everything you've said. Ideally, the derived classes should not override the sample method. However, the problems I ran into with that are:

• For IndependentSource, as you know the constraint checking is interspersed with the sampling logic and there are "early outs", where e.g. if the position is rejected, there's no point in sampling an angle and energy. We could remove this but 1) it would be less efficient and 2) it would result in all our test results changing (not the biggest deal ever, but a consideration nonetheless).
• For MeshSource, there's a motivation for overriding sample that is a bit deeper than just it being convenient. There are two steps to the mesh source sampling: 1) sample a mesh element, and 2) sample a position within that element. If you apply domain rejection after both sampling an element and position within the element, you can end up biasing the the original distribution of mesh element probabilities. Take the following thought experiment: you have a mesh with two elements, one of which has probability 0.1 and the other of which has probability 0.9. However, within element 1, there is 0% chance of rejecting based on the domain whereas within element 2, there is a 99% probability of rejecting. In this case, if you sample element+position and then apply domain rejection, you'll end up grossly oversampling the source in the first mesh element.

If you have ideas for how to address those two points while having sample be non-virtual, I'm all ears. In the meantime, I'll try to address all the other comments while we stew on what the right solution is for the virtual sample method issue.

[gridley]

Hey Paul! OK, yes, that makes sense now. I agree that it'd be a bummer to not take advantage of some efficiency considerations in rejection in favor of the most separated responsibility code design.

Regardless, I think you could still make sample non-virtual. It's because nothing precludes you from using rejection in sample_no_rejection. Wouldn't that let the tests pass then, if rejection was applied in both places? Maybe that's a little icky code design though.

So are you saying the objective in mesh source sampling is to preserve the element sampling probabilities even after applying a rejection filter? I am not 100% understanding that one right now. I would have guessed it to be natural to have different element sampling probabilities in the end after applying a rejection filter, but I can imagine there are applications where it's desirable to preserve these after applying rejection.

[paulromano]

@gridley I've updated this branch in response to all your excellent comments. Notably, the interaction between the Source class and derived classes is now cleaner. Derived classes are required to implement the sample method (that provides a "proposal" sampled site) and then the source rejection is done in the non-virtual sample_with_constraints method from the base class. In order to get IndependentSource/MeshSource to work with that, they check for constraints directly in their sample method and then there's a flag that tells the base class not to do a redundant check against the constraints.

I've also added a generic "fissionable" constraint that can be used with any source class and that replaces the only_fissionable argument on openmc.stats.Box (which is hereby deprecated).

Let me know if you see any further issues!

[ebknudsen]

Just wanted to say - thanks guys for all the work you've put in on this one. Much appreciated.