-
Notifications
You must be signed in to change notification settings - Fork 526
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Pulse height tally for photons #2452
Pulse height tally for photons #2452
Conversation
I'm impressed by how few lines of code this is! Sweet! |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for the updated PR here @cfichtlscherer. Like @gridley I'm impressed by how small the implementation is! There are some design issues we'll need to work through but I'm hopeful we can get to the point of having something that can be merged.
Also, FYI it looks like you have tests failing in CI which you'll need to look into. |
Yes, I think this should be solved by 09de63d |
Thanks a lot @gridley @paulromano Yes we dont need the special estimator anymore which is much better for the general code structure I think. Currently I am still fighting with another bug but I am confident to show you the next version soon :-D |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@cfichtlscherer Thanks for the updates. Nice to see that the filters are essentially unchanged now! A few more requests but this is getting really close.
One thing I would like to see before we merge: I know that you've previously done comparisons with MCNP. Would it be easy enough to repeat this comparison using the current version in this PR just to make sure things still look good?
src/particle.cpp
Outdated
if (cell_born() == C_NONE) | ||
cell_born() = coord(n_coord() - 1).cell; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ah, sorry for my misunderstanding. Now I see that this code is actually exactly what's being called in event_calculate_xs
, so no behavior is changing. That being said, it still feels awkward to me that this block of code is called in pht_secondary_particles
. Can you move it up into event_revive_from_secondary
?
src/tallies/tally.cpp
Outdated
@@ -381,6 +393,10 @@ void Tally::add_filter(Filter* filter) | |||
energyout_filter_ = filters_.size(); | |||
} else if (dynamic_cast<const DelayedGroupFilter*>(filter)) { | |||
delayedgroup_filter_ = filters_.size(); | |||
} else if (dynamic_cast<const CellFilter*>(filter)) { | |||
cell_filter_ = filters_.size(); | |||
} else if (dynamic_cast<const EnergyFilter*>(filter)) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Still need updates here on the type checks
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for your patience @cfichtlscherer. Hopefully the last round of comments here!
src/tallies/tally.cpp
Outdated
@@ -381,6 +393,10 @@ void Tally::add_filter(Filter* filter) | |||
energyout_filter_ = filters_.size(); | |||
} else if (dynamic_cast<const DelayedGroupFilter*>(filter)) { | |||
delayedgroup_filter_ = filters_.size(); | |||
} else if (dynamic_cast<const CellFilter*>(filter)) { | |||
cell_filter_ = filters_.size(); | |||
} else if (dynamic_cast<const EnergyFilter*>(filter)) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, I would recommend we change all of them here to not use dynamic_cast
and instead rely on filt->type()
checks.
All done, think we can merge. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for all the work getting this to the finish line!
This pull request represents a second attempt to provide OpenMC with the capability of a pulse-height tally (see
#1881). This non-Boltzmann tally allows the simulation of all kinds of gamma detector measurements.
Thanks to @gridley and @paulromano for their help and input during the implementation. The feature is extensively validated and produces the same results as MCNP.
This tally only makes sense if photon transport is activated and works only in combination with an
EnergyFilter
and aCellFilter
- checked at the beginning of the simulation intally.cpp
.Furthermore, the variable
vector<int> active_pulse_height_tallies;
is introduced, which is used to check if the simulation tallies the pulse-height.The basic idea for storing the information of the delivered energy over the entire history of a particle is to give the particle class a new attribute,
pht_storage
(a vector with the length of the number of cells in the system). In this vector, the pulse-height contributions are updated whenever the photon performs a collision or secondary particles are started. At the end of the particle's history, the functionscore_pulse_height_tally
is called. While scoring the pulse-height, we must always iterate over all selected cells. Even if the pulse-height value is zero, this result must be scored. For this reason also theFilterCell
class was changed.I am very happy and thankful for all feedback and hope we can have this feature in the official OpenMC code soon!
When this pull-request is merged, I would happily create a sample page (jupyter notebook) that shows how to quickly and directly simulate detector responses in your system.
Update July 7, 2023
The results of OpenMC pulse-height tally simulations fit very well with the corresponding MCNP6.2 simulations (F8 tally).
For the following plot, photons with the energy of 1MeV were started from a disk source into a sodium-iodine scintillation detector. The identical model was implemented in OpenMC and MCNP6.2. The MCNP simulation used 1e8 particles and ENDFB_VI_8 crosssections, the OpenMC simulation 1e9 particles, and the Lanl_based/mcnp_endfb71 crosssections.
I am happy to share input files and more details with anyone who is interested.
![Screenshot from 2023-07-07 18-49-56](https://private-user-images.githubusercontent.com/29277544/251806125-d9253fe2-cbb6-429f-bc35-a6898d9d2358.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MzkwODA2MTMsIm5iZiI6MTczOTA4MDMxMywicGF0aCI6Ii8yOTI3NzU0NC8yNTE4MDYxMjUtZDkyNTNmZTItY2JiNi00MjlmLWJjMzUtYTY4OThkOWQyMzU4LnBuZz9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNTAyMDklMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjUwMjA5VDA1NTE1M1omWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPTI5YjE3NjNmYmQ3NTg1ZDIwYTVkYzUyODNjMzBhZDhjNzRkZTA5NjdjNjM3M2U5ZjExOTUzNmJkNGRmZmRmMTAmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0In0.dWUywcXi0xIDwwqXjm9EcSKPBVNz4OE9NSgUvV4_uOk)