Storing surface source points using a cell ID

[JoffreyDorville]

Description

Currently, the ‘surf_source_write’ setting stores information from particles that are crossing surfaces of interest into an HDF5 file. These surfaces of interest are defined by the user with the ‘surface_ids’ parameter. The main objective of this development is to also provide the user with a ‘cell’ parameter to easily select surfaces from a cell ID without defining any surface IDs. This new ‘cell’ parameter can work on its own or can be used in combination with surface IDs. Additionally, two new arguments, ‘cellfrom’ and ‘cellto’, have been created to filter particles based on their direction.

Motivation:

Let’s consider the following simplified model (Fig. 1) containing a cylindrical source (in green), a first cube (in purple) and a second cube (in red) defined with a vacuum boundary condition:

Figure 1. Simplified model

I was interested in obtaining the particles that were crossing the outer surface of the purple box. I then used the available feature by declaring the corresponding surface IDs. Using the generated source file, 2D projections of the points distribution show that OpenMC effectively stores points at every declared surface (Fig 2.).

Figure 2. Correct results from the current feature

However, I only wanted to store points that were crossing the physical outer surface of the first box, not every individual surface used to construct the outer boundary of the cell. Compared to the previous 2D projections (Fig 2.), here is what I was looking for (Fig. 3):

Figure 3. Results I was looking for

New feature:

The proposed feature is to add a ‘cell’ argument to the ‘surf_source_write’ setting so that OpenMC only stores surface source points when particles are actually leaving or entering the cell of interest.

Using the previous example and knowing that the cell ID of the first box is 2, we can use ‘surf_source_write’ as follow to obtain the expected results (Fig. 3):

settings.surf_source_write = {
	'max_particles': 5000,
	'surface_ids': [4, 5, 6, 7, 8, 9],
	'cell': 2
}

With the new feature, if no surface indices are declared, OpenMC will store particles crossing any surfaces of the model. This means that the following call will store any particles crossing any surfaces of the model:

settings.surf_source_write = {
	'max_particles': 5000
}

This also means that the following call will store particles that are crossing any surfaces only if they are coming from or going to the user-defined cell:

settings.surf_source_write = {
	'max_particles': 5000,
	'cell': 2
}

It is important to note that with this call for the case of our first box, it will score points on the outer surface but also on the inner surface delimited by the cylinder.

Additionally, two new arguments, ‘cellfrom’ and ‘cellto’, have been created to filter particles based on their direction. The call to these arguments is syntactically similar to the call to the ‘cell’ argument:

settings.surf_source_write = {
	'max_particles': 5000,
	'cellfrom': 2
}

settings.surf_source_write = {
	'max_particles': 5000,
	'cellto': 2
}

The three new arguments ‘cell’, ‘cellfrom’ and ‘cellto’ cannot be used simultaneously and only one cell ID can be declared at a time for these arguments.

A regression test containing 20 use cases (non exhaustive) has been added to demonstrate the feature. A ‘_visualize.py” script has been added to this regression test folder to visualize the source distribution in 3D. Particular attention has been given to models with multiple level of coordinates but this new feature has not been specifically tested on lattices.

A current limitation with the new feature is that ‘surf_source_write’ cannot scores points at a surface with ‘periodic’ or ‘reflective’ boundary condition while using any of the ‘cell’, ‘cellfrom’ or ‘cellto’ arguments. If such surfaces were to be selected the points associated with these surfaces are simply not stored. For backward compatibility, storing points crossing a ‘periodic’ or ‘reflective’ surface when only surface IDs are selected is still available and is still done before handling the boundary conditions. Using ‘periodic’ and ‘reflective’ boundary conditions with the cell ID could be addressed as a further enhancement in a following PR.

Additional corrections:

During the development of this feature, I had to adjust particle.cpp to populate the ‘cell_last’ attribute of the ‘Particle’ class before the level of coordinates (n_coord) was updated. This previous behavior could potentially hide previously visited cells ID to the CellFrom tally filter if the particle was crossing a surface between two cells with different number of level of coordinates.

Additionally, the ‘cell_last’ attribute was not correctly initialized and was pointing to a cell index of zero. It could potentially lead to scoring contribution from the cell with index zero even if it does not make sense physically. In this PR, I propose to initialize the values at ‘C_NONE’ and update the content of ‘cell_last’ at birth so that it reflects the current cells (for all level of coordinates).

A regression test has been added for the CellFrom filter to check that the contributions are correctly calculated in a model with multiple level of coordinates. However, it has not been specifically tested on lattices.

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)

[gridley]

I thought the purpose of surface sources was for moving from something like a reactor problem to a reactor shielding problem. Why not just set the surfaces you're banking to to be vacuum type? I was under the impression this is always how people use surface_source_write.

[gridley]

Come to think of it--for writing surface sources, if you're not stopping the particle at the surface you write particles to the bank at, do you not run the risk of double-counting particles? It would consequently give the wrong source intensities for shielding problems. Admittedly I may not be an expert on shielding and may be 100% wrong, but this superficially seems like a problem to be solved by changing the overarching simulation technique, not by adding 4000 lines of code.

[gridley]

Lastly: this functionality could be achieved without any added C++ as a postprocessing step. You'd just open the surface source h5 file with a python script, loop over the sites, check the cell of each location via the C API, and reject sites that fall outside the region you're interested in.

I apologize to perhaps come across as hostile towards what appears to be quite a large amount of work you've done here.

[paulromano]

@JoffreyDorville Very nice work on this enhancement to our surface source writing capability, which gives users more flexibility in specifying how particles should be banked. I would disagree with @gridley's suggestion that this could be done just as effectively as a post-processing step; that would entail collecting a lot more data than the user really needs if they are interested in a small subset of the particles that are crossing a given surface.

The overall changes to the C++ code are modest (<300 lines of code in the src/directory) and the testing you've added here is very extensive. I've made a few requests below, the most important of which is trying to eliminate the runtime building of a set. I've also made some direct changes on your branch to simplify/improve the code a bit, so let me know if you have any questions about that!

[gridley]

@paulromano Well sure, it would produce more data than required if treated as a postprocessing step, true.

But the more fundamental reason underlying my complaint is that this is a misuse of writing to a surface source, and can in fact cause double-counting for downstream simulations. I can't imagine any scenario where is's preferable to not just set the surface we're writing to as a vacuum BC, then have the second calculation pick up from there.

Why would you not just set it as a vacuum BC if generating a surface source in this case? For particles that cross out of the surface, scatter back in, then exit again, this causes double counting to the bank. It's one history in that case, whereas in the case of terminating the particle after writing it to a bank, you know for sure it's not being double counted.

So to be explicit: in the example provided you could have just set the surfaces you're writing banks at to be vacuum BCs and achieved exactly the same result @JoffreyDorville. Sorry if this has come across as hostile, and maybe I should have opened with more questions because perhaps I'm missing something about why this type of simulation approach is necessary.

[paulromano]

But the more fundamental reason underlying my complaint is that this is a misuse of writing to a surface source, and can in fact cause double-counting for downstream simulations.

Your complaint is not about the feature itself as implemented here. The example @JoffreyDorville gave is exactly that—an example. I can imagine lots of ways the feature here is useful beyond the particular motivating example. The main questions a reviewer should be answering are:

• Is the feature generally useful?
• Does it affect maintainability of the code? Is it an intrusive change?
• Does it have a negative impact on performance?
• etc.

If you have comments or suggestions about the implementation of the feature that has been added here, those would be most welcome. If you want discuss the applicability of this feature to the particular use case @JoffreyDorville raised in his description, I would encourage you to do so separately from this PR (discourse?) so that we can focus on the code changes.

[gridley]

Alright I don't mean to be stubborn, but I feel my comments have been addressing exactly those questions.

The question I am answering with nearly every comment I've made is pertaining to "is the feature useful?"

And every point raised is substantiating why I believe the answer to that question is no. It is not useful because 1) this can be achieved through other nonintrusive means as it appears to be an extremely niche situation and 2) appears to not even give correct answers as it would result in double-counting particles. It degrades source code maintainability because it introduces a substantial delta of code (cognitive burden on future developers) which will be touched infrequently, and lastly indeed negatively impacts performance in the current implementation due to allocating memory on the fly. And this is all the basic substance of what my previous comments convey.

[paulromano]

this can be achieved through other nonintrusive means as it appears to be an extremely niche situation

If a surface bounds multiple cells on a given side, there is no way to differentiate particles entering/leaving a specific cell. The additions in this PR make it possible to filter by neighboring cell in addition to a surface ID. I don't believe this to be a niche situation. As a point of reference, MCNP also allows specifying entering/exiting cell on the SSW card.

appears to not even give correct answers as it would result in double-counting particles

"Correct" is up to the user. The feature simply banks particles on a surface and gives that information to the user. How they use that information is not our problem. We are not endorsing any specific methodology here, only supplying the capability to obtain information on particles passing through a surface.

It degrades source code maintainability because it introduces a substantial delta of code

Believe me when I tell you that I am sensitive to accepting changes to the code that would impact maintainability because—let's be honest—that is a burden that would fall on me. The changes in this PR are really not concerning to me. Most of the lines of change are simply tests. There are 233 lines changed in .cpp files, and that can probably be reduced further with some refactoring. Most of those lines are also contained in a single function that handles banking particles to the surface source.

negatively impacts performance in the current implementation due to allocating memory on the fly

Agreed -- I also noted this in my review, but I think that can be fixed easily and won't be a showstopper.

[gridley]

OK Paul you make some good points. You've swayed my position with the MCNP point. To clarify my position on the first point in your reply I'm claiming that surface source writing should be carried out with vacuum BCs on the surface each time, which immediately solves the problem here.

But if MCNP supports it, I'm sure there's a good reason to allow this. I just have no idea what that could possibly be.

On "it'll fall on me..." 😆

[gridley]

I was reflecting on this some more this morning, and still can't think of a situation in which this is a preferable or more correct approach over applying vacuum boundary conditions to the surface particles are being saved at.

Am I misinterpreting the purpose of surface sources? I was under the impression this is used to simulate leakage through some complicated region with a source inside, perhaps a reactor of some sort, then move to a separate calculation for simulating the externals.

I apologize for the insistence on debating this, but after Paul's comment I have been racking my mind to figure out why this might be necessary or useful.

[paulromano]

@gridley You are not misinterpreting the purpose of surface sources. The main use case is capturing the particle current from some system (where it is assumed it is "difficult" for particles to get from the original source to the surface) and then using that as a starting source for simulating some external geometry.

As @JoffreyDorville and I (and @shikhar413) have been iterating on a particular use case for this functionality, we ourselves have learned more about the limitations and I feel more confident saying that just imposing vacuum boundaries is generally not the right thing to do when there is any influence of the external geometry on the starting geometry. That is, if you have backscatter / reflection of particles coming back into your original geometry, imposing a vacuum boundary exactly where you are collecting a surface source will result in an incorrect current on the surface (also if you are doing a k-eigenvalue calculation on your original geometry, that backscattering can also affect k, which is another issue). To model things correctly, the surface at which you are collecting particles should be a bit inward from wherever particles are getting terminated. So, you might have a setup something like this:

where the external geometry starts at the surface source but the original geometry penetrates a bit into that external geometry to capture backscattering. Validating this point is the following sentence from the MCNP 6.3 manual (section 5.8.8):

Include enough geometry beyond the specified surfaces to account for albedo effects.

[gridley]

@paulromano thank you very much for clarifying that. I suppose that either approach is imperfect: with what I've suggested some small eigenvalue error is incurred, and with the one suggested in the MCNP manual it still seems to me like you can double count particles that leak, reflect back, then pass the surface again.

At least using the cell-to-cell marking method should cut down the amount of double-counted particles substantially.

If we want to correctly normalize the surface source against the source rate in the internal region, don't we need a method to prevent double-counting? For example in a reactor shielding problem you would have a total power tally, multiply the source to your desired thermal power, then for follow-on shielding calcs (in a full power problem for example) you'd need to know that normalizing factor: the total number of surface source points per source neutron.

As for the vacuum approach: if you tally surface sources on a vacuum BC, you're guaranteed to stop the history there and prevent double-counting. For the follow-on calculations, you'd turn off fission in the internal region, then allow particles that reflect back in to do that as they please. They might reflect back out, but if you've turned off fission and normalize to a certain power you'll be using the correct source neutron rate for the surface source.

I agree at least that if MCNP has this feature, at least, it's worth having.

Lastly, I disagree with the MCNP manual here. A vacuum BC should be used, and the surface on which surface sources are tallied should be far enough out so that albedo is correctly accounted for. Maybe MCNP accounts for double-counting by marking the history ID of a particle banked on the surface.

[gridley]

@paulromano in your example why wouldn't you just make the surface to collect particles at the outer one and not the inner one?

[gridley]

[paulromano]

@gridley Thank you for your careful consideration and for again coming back to the issue of double-counting. I think I've come to a clear mental picture of this now. To avoid double-counting you either need to have a vacuum boundary in the first calculation or in the second calculation but not in both. So, the options are either:

1. Run the original geometry with a vacuum BC between original and external (i.e., external geometry is excluded). Then run a calculation with the surface source into external but include the original geometry so that particles can scatter back and forth.
2. Run the original geometry and include the external geometry, collecting the surface source at the interface between original and external. Then run a calculation with the surface source into external but exclude the original geometry by imposing a vacuum BC on particles that would have gone back into the original geometry.

The picture I drew before is a slight variation on option 2 where the first calculation includes enough of the external geometry to capture the albedo but not the full thing.

[gridley]

Ah! Well that makes perfect sense to me.

On top of that, I thought of another case where this new approach certainly makes more sense than using a vacuum BC. It would be useful if you knew ahead of time that the dominant contribution to shield originates from some penetration into a shield, and wanted to only tally what's leaving through that penetration.

Sorry to be stubborn about this but I was genuinely convinced until now that this approach was not only unnecessary, but also incorrect. After digesting it for a while, it makes sense to me.

[paulromano]

I'm glad that you were stubborn because I think we all came to a better understanding of it in the end!

[paulromano]

Thanks for all the updates on this @JoffreyDorville. I made a few small changes myself and am happy with where it is now so I'll go ahead and approve.

@gridley Please let us know if you have any further concerns / requests before merging.

[paulromano]

@gridley Just pinging you once again to make sure you have no further requests before we merge.

[gridley]

Nope, I'll go ahead and hit the button. Sorry to miss that last notification!

[yrrepy]

I'm just seeing this now.
This feature will be extremely useful to me.
Thank you @JoffreyDorville

Some comments from my perspective:
Often in SSR, it is up to the user to ensure they properly setup the simulation to avoid double-counting particles. Post-processing is laborious and these files can become very large.

I often have complex surfaces over which I desire the SSW, a plane with cone protruding, and for maximum acceleration I want the particles starting on the plane surface and on the complex cone surface.

I agree with Paul's final two bullets.
In my application, I have to use 1), as the external geometry changes from case to case.
IIRC, 2) was best for reactor geometries where I am performing deep penetration and want to avoid re-simulating transport/albedo back into the core.