Restricted file source #2916

ebknudsen · 2024-03-14T15:00:15Z

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:

"resample" a new particle is picked at random, until a valid particle has been accepted. N.b. this may potentially lead to biasing.
"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

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)

it appears that before that geometry is not initialized

add an option to, if a particle is rejected, kill that one and move to next particle instead of resampling. This is less efficient, but avoids renormalization

…into restricted_file_source

also apply a bit of formatting

…into restricted_file_source

gridley · 2024-03-16T04:03:57Z

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 · 2024-03-19T09:18:58Z

@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.

correct wrong condition

instead of a full coordinate vector instead use lower_left upper_right for spatial coordinates, and energy_bounds and time_bounds for energy and time. This should be more workable

ebknudsen · 2024-03-20T08:40:56Z

@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 · 2024-03-20T12:59:09Z

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 · 2024-03-20T13:10:14Z

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!

The source may now be restricted if it inherits from RestrictedSource instead of directly from Source

paulromano · 2024-05-09T02:48:24Z

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 · 2024-05-09T09:05:49Z

@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.

gridley

Alright, I've reviewed the code! Here is my overarching concern.

The Source class should probably be 100% responsible for site rejection. We have left IndependentSource with its own ability to reject based on if the location is in fissionable material, but that should also go into Source if you ask me!

Also, it appears that the rejection constraints that are implemented in Source are getting bypassed by IndependentSource. Ideally putting in this work to generalize rejection sampling should also be enjoyed by any child classes. Moreover, child classes should not have the option to override sample. It's a violation of the open-closed principle: we are allowing the class behavior to be modified rather than extended.

On top of that we're not clearly separating responsibilities at this point. The Source class should be responsible for site rejection, but we've left IndependentSource responsible for doing non-fissionable region rejections.

I know we're not hard and fast on developer documentation at this point, but going into the future it'd be nice to write a paragraph about the general structure of the Source system in the C++ code. It right now is a little opaque that child classes should be responsible for proposing source sites, while the parent class potentially rejects them.

gridley · 2024-05-09T14:55:38Z

include/openmc/source.h


  static unique_ptr<Source> create(pugi::xml_node node);
+
+protected:
+  // Methods that can be overridden


This comment is kind of extraneous. I don't think basic language features should be documented, but will leave this suggestion up to your discretion.

Instead, what is this function supposed to do? What would a base class have to implement? A comment should be there saying it I think.

gridley · 2024-05-09T14:58:43Z

include/openmc/source.h

+  void check_for_constraints(pugi::xml_node node);
+  bool satisfies_constraints(SourceSite& s) const;
+  bool satisfies_spatial_constraints(Position r) const;
+  bool satisfies_energy_constraints(const double E) const;


Marking values we pass in as const doesn't hurt; just as an observation I don't think there's much of anywhere else we do this.

For types longer than 8 bytes or so, such as the Position there, it might be good to pass as a const reference instead.

include/openmc/source.h

gridley · 2024-05-09T15:05:07Z

include/openmc/source.h

+protected:
+  // Note that this method is not used since IndependentSource directly
+  // overrides the sample method
+  SourceSite sample_no_rejection(uint64_t* seed) const override { return {}; }


Oh OK! I see what sample no rejection does. This is the proposal distribution. The first thing that came to mind for me reading "no_rejection" is that this implements some other sampling method that avoids using rejection. Instead it's the proposal distribution inside a rejection loop. I'd suggest renaming to "sample_proposal" or similar.

This should definitely be documented in the base class with a brief explainer on how it interacts with child classes.

gridley · 2024-05-09T15:07:20Z

include/openmc/source.h

+protected:
+  // Note that this method is not used since IndependentSource directly
+  // overrides the sample method
+  SourceSite sample_no_rejection(uint64_t* seed) const override { return {}; }


Moreover, I find it funny that we've taken the base class route here, so ideally IndependentSource should NOT override sample. That way it retains the ability to do source site rejection.

It appears to me that this approach does not allow source site rejection with IndependentSource. That's no good! The sample method in the base class should NOT be virtual. Rather, base classes should only fill out the proposal distribution. This way we guarantee rejection functionality on any source class.

gridley · 2024-05-09T15:40:35Z

src/source.cpp

+      }
+    }
+  }
+  // If at this point found is false, part. is outside all cells/materials


Let's not abbreviate particle. I was thinking this was a typo period and referring to a CAD part!

gridley · 2024-05-09T15:42:26Z

src/source.cpp

-    // Now search to see if location exists in geometry
-    found = exhaustive_find_cell(p);
+    // Check if sampled position satisfies spatial constraints
+    found = satisfies_spatial_constraints(site.r);

    // Check if spatial site is in fissionable material


Let's move this fissionable check to the outer Source code. I really deeply believe that's the right thing to do here.

Why should an IndependentSource enjoy fissionable rejection but not a FileSource?

gridley · 2024-05-09T15:46:15Z

src/source.cpp

+
+  // Sample position and apply rejection on spatial domains
+  Position r;
+  while (true) {


Rather than a while/break, I'd suggest hearkening back to our Fortran roots here with the do/while construct more commonly used in that language.

gridley · 2024-05-09T15:47:14Z

src/source.cpp

+  Position r;
+  while (true) {
+    r = space_->mesh()->sample_element(element, seed);
+    if (this->satisfies_spatial_constraints(r)) {


OK, I see why you'd want to separate this out as its own function now. It's nice to pre-screen out some rejection sites before passing back to the Source class's rejection loop.

gridley · 2024-05-09T15:47:37Z

src/source.cpp

-  // Replace spatial position with the one already sampled
-  site.r = mesh_location.second;
+  SourceSite site;
+  while (true) {


Again do/while is a natural, beautiful way to handle this type of loop.

It is, but when I tried this, clang-format put the closing semicolon on a line by itself, which looks both ugly and confusing to me. Sooo, I'm going to leave this one as is (I will use do/while above though)

paulromano · 2024-05-09T16:21:54Z

@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 · 2024-05-09T16:28:05Z

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 · 2024-05-13T17:32:52Z

@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!

gridley

I like how you handled potentially mixed constraint handling! It looks 99% of the way good to go but I have a few small questions before approving.

include/openmc/source.h

src/source.cpp

gridley · 2024-05-14T13:30:35Z

src/source.cpp

-        accepted = true;
-        site.wgt = 0.0;
+      if (!accepted) {
+        ++n_reject;


I believe this is not thread safe use of a static variable, although the probability of two threads simultaneously hitting it is low. I suppose that having these be 32 bit types means you'll never get unreasonable values. Maybe the counters should be confined to a single thread? I know it doesn't influence correctness but worth a little thought.

I also had this realization looking through the code and thought about how much we should care about it (the same thing applies in the existing logic in IndependentSource). The important thing is to trigger the condition when too many sites are getting rejected; even if there is a race condition on incrementing the number of accepts/rejects, that will still happen, just in a semi-fickle way. I was going to defer fixing it for now but if you think we should take care of it as part of this PR, I'm happy to do so. Let me know how I should proceed.

Sounds good! I'm just going to raise a separate issue to handle it.

ebknudsen · 2024-05-15T08:02:20Z

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

gridley

LGTM, but maybe we can #include <cassert> and add assert on to the pointers we're getting from dynamic casts? cassert is a small header, so there's no tangible overhead from adding it. Alternatively dynamic casting to a reference rather than pointer would suffice.

It's better than the GSL Expects approach because they're removed in any compile that's not debug:

https://stackoverflow.com/questions/34302265/does-cmake-build-type-release-imply-dndebug

Co-authored-by: church89 <l.chierici89@gmail.com> Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

ebknudsen and others added 20 commits March 11, 2024 01:18

add code to reject on domains also for file sources

3df948b

also add the domain restriction to python layer

9e850f6

move rejection to the sampling step

35e160c

it appears that before that geometry is not initialized

improve commenting and formatting

0fbc48d

fix flipped logic

471a4aa

add option to kill non-conforming particles

972a281

add an option to, if a particle is rejected, kill that one and move to next particle instead of resampling. This is less efficient, but avoids renormalization

add option to reject by hypercube of particles

b1537a7

allow hypercube description in python layer

e12732b

fix typos

03e122c

fix typos

1baf072

yet another misprint

916315e

Merge branch 'restricted_file_source' of github.com:ebknudsen/openmc …

248353d

…into restricted_file_source

hypercube and rejection_strategy added to python

672f4e6

also apply a bit of formatting

Merge branch 'openmc-dev:develop' into restricted_file_source

11da3c5

adjust whitespace acc. to clang-format

ed46945

Merge branch 'restricted_file_source' of github.com:ebknudsen/openmc …

183779f

…into restricted_file_source

remove unimplemented FILTER option

5a892f8

fix var name

eed4aaf

initialize attributes to None

bb8e892

initialize the rej. stragtegy and add doc string

23a29b4

church89 and others added 3 commits March 19, 2024 16:18

correct wrong condition

20b27f8

Merge pull request openmc-dev#29 from openmsr/correct_iterator

c188fab

correct wrong condition

move from a single vector to discrete bounds

d93ac83

instead of a full coordinate vector instead use lower_left upper_right for spatial coordinates, and energy_bounds and time_bounds for energy and time. This should be more workable

ebknudsen added 2 commits April 8, 2024 13:50

source restriction is now an intermediate class

86ae064

The source may now be restricted if it inherits from RestrictedSource instead of directly from Source

apply clang-format

78a8716

gridley requested changes May 9, 2024

View reviewed changes

paulromano added 7 commits May 10, 2024 07:22

Respond to many @gridley comments

0e2f747

Add fissionable source constraint

cd04ecb

Rename to sample / sample_with_constraints

91cfdcc

Make sample_with_constraints non-virtual

a60ff90

Add user's guide section on source constraints

e3f56e7

Put source constraints into a <constraints> XML element

e468dd4

Add comment describing Source base/derived relationship

42d5cad

gridley reviewed May 13, 2024

View reviewed changes

include/openmc/source.h Outdated Show resolved Hide resolved

src/source.cpp Show resolved Hide resolved

paulromano added 3 commits May 13, 2024 21:20

Update reference inputs to reflect contraints format

c9d130d

Fix one more reference input file

bad3a01

Implement rejection fraction check in sample_with_constraints

42898cc

gridley reviewed May 14, 2024

View reviewed changes

Fix rejection threshold check in sample_with_constraints

ca3ba80

gridley approved these changes May 14, 2024

View reviewed changes

gridley mentioned this pull request May 14, 2024

(potential bug) Use atomic add operations for incrementing source site rejection counters #3005

Open

Skip rejection limit test for MPI config

da08d77

paulromano merged commit 1702b45 into openmc-dev:develop May 15, 2024
18 checks passed

paulromano mentioned this pull request Jun 2, 2024

Make sure direction is set when checking source spatial constraints #3022

Merged

5 tasks

paulromano mentioned this pull request Jun 11, 2024

Fixed Source Random Ray #2988

Merged

5 tasks

gridley reviewed Jun 12, 2024

View reviewed changes

This was referenced Jul 5, 2024

Update openmc neams-th-coe/cardinal#903

Closed

Fission site bounding box reporting when threshold fissionable site rejection sampling error is triggered #3073

Open

church89 added a commit to openmsr/openmc that referenced this pull request Jul 18, 2024

Restricted file source (openmc-dev#2916)

86d7fe6

Co-authored-by: church89 <l.chierici89@gmail.com> Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Restricted file source #2916

Restricted file source #2916

ebknudsen commented Mar 14, 2024 •

edited by paulromano

Loading

gridley commented Mar 16, 2024 •

edited

Loading

ebknudsen commented Mar 19, 2024 •

edited

Loading

ebknudsen commented Mar 20, 2024

gridley commented Mar 20, 2024

ebknudsen commented Mar 20, 2024

paulromano commented May 9, 2024

ebknudsen commented May 9, 2024 •

edited

Loading

gridley left a comment

gridley May 9, 2024

gridley May 9, 2024

gridley May 9, 2024

gridley May 9, 2024

gridley May 9, 2024

gridley May 9, 2024

gridley May 9, 2024

gridley May 9, 2024

gridley May 9, 2024

paulromano May 9, 2024

paulromano commented May 9, 2024

gridley commented May 9, 2024

paulromano commented May 13, 2024

gridley left a comment

gridley May 14, 2024

paulromano May 14, 2024

gridley May 14, 2024

ebknudsen commented May 15, 2024

gridley left a comment •

edited

Loading

Restricted file source #2916

Restricted file source #2916

Conversation

ebknudsen commented Mar 14, 2024 • edited by paulromano Loading

Description

Checklist

gridley commented Mar 16, 2024 • edited Loading

ebknudsen commented Mar 19, 2024 • edited Loading

ebknudsen commented Mar 20, 2024

gridley commented Mar 20, 2024

ebknudsen commented Mar 20, 2024

paulromano commented May 9, 2024

ebknudsen commented May 9, 2024 • edited Loading

gridley left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

paulromano commented May 9, 2024

gridley commented May 9, 2024

paulromano commented May 13, 2024

gridley left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ebknudsen commented May 15, 2024

gridley left a comment • edited Loading

Choose a reason for hiding this comment

ebknudsen commented Mar 14, 2024 •

edited by paulromano

Loading

gridley commented Mar 16, 2024 •

edited

Loading

ebknudsen commented Mar 19, 2024 •

edited

Loading

ebknudsen commented May 9, 2024 •

edited

Loading

gridley left a comment •

edited

Loading