openmc icon indicating copy to clipboard operation
openmc copied to clipboard

Adding support for NCrystal materials in OpenMC

Open marquezj opened this issue 2 years ago • 3 comments

After some discussion with @paulromano we cleaned up this version of OpenMC with support for NCrystal. NCrystal is a software library that provides low energy neutron physics to Monte Carlo codes. It was initially designed with Geant4 in mind, but it currently has interfaces with many other codes, and it can be called from Python, C and C++. NCrystal is developed mainly at the European Spallation Source. You can read more about NCrystal in these papers:

X.-X. Cai and T. Kittelmann, NCrystal: A library for thermal neutron transport, Computer Physics Communications 246 (2020) 106851, https://doi.org/10.1016/j.cpc.2019.07.015

T. Kittelmann and X.-X. Cai, Elastic neutron scattering models for NCrystal, Computer Physics Communications 267 (2021) 108082, https://doi.org/10.1016/j.cpc.2021.108082

X.-X. Cai, T. Kittelmann, et. al., "Rejection-based sampling of inelastic neutron scattering", Journal of Computational Physics 380 (2019) 400-407, https://doi.org/10.1016/j.jcp.2018.11.043

or watch the presentation from @tkittel, at ND2022:

https://indico.frib.msu.edu/event/52/contributions/954/

We have used NCrystal together with NJOY to generate a large library of new materials in ACE format (the NJOY-NCrystal library), but calling NCrystal directly allows to introduce in OpenMC additional physics that is not supported by the ACE format.

In particular, NCrystal allows to generate and sample on the fly scattering kernels for solid polycrystalline materials at any temperature. This is done at initialization, based on a human readable description stored in text form in the NCMAT format. For existing ENDF-6 scattering kernels, it allows the direct sampling of S(alpha, beta), which allows for continuous angular distributions (which are still not supported by the ACE format). NCrystal also support oriented materials (single crystals, highly-oriented pyrolytic graphite) and has plug-in support for other physics, currently being used to include small angle neutron scattering and inelastic magnetic scattering. All this would make OpenMC more useful to the neutron scattering community, while adding advanced nuclear data capabilities for the reactor community.

This pull request is planned as an optional feature that could be activated with cmake, using the OPENMC_USE_NCRYSTAL flag after installing and initializing NCrystal. A new method is added to the Material class in Python, to create materials from NCrystal, which takes the composition and density, creates the OpenMC material and associates it with NCrystal.

After reading the xml files, an NCrystal object is created in memory, which provides thermal scattering cross sections and scattering events (for the calculation of the scattering events, the direction of the neutron is a parameter to be able to consider oriented materials). When the thermal scattering cross section is computed, the nuclear elastic scattering component is removed, similarly to what is done in the standard S(alpha, beta) treatment.

There a couple of open topics that we could not decide, and are marked with "fixme" in the code. For instance, NCrystal uses a cache that could be implemented in OpenMC to reduce the computation time. Also, new tags for NCrystal events could be added for tallies. But, overall, we think it is ready for merging upstream.

marquezj avatar Sep 14 '22 14:09 marquezj

I can add to what @marquezj wrote, that since last week NCrystal is also available on conda-forge (since @paulromano suggested this during our meeting).

tkittel avatar Sep 14 '22 15:09 tkittel

Thanks for the PR @marquezj! Really excited to see this. Just wanted to give you a heads up I probably won't get to reviewing this until early next week.

paulromano avatar Sep 15 '22 12:09 paulromano

@paulromano I added comments concerning the two non-obvious fixmes in the PR. Let us know if you have any specific wishes concerning this (or anything else), and we will try to update the PR.

tkittel avatar Sep 16 '22 07:09 tkittel

So @paulromano, unless I overlooked something I believe @marquezj took care of all the changes requested so far (including applying clang-format).

tkittel avatar Sep 26 '22 11:09 tkittel

FYI we haven't forgotten this, but some other work came in the way. We hope to spend some time on this on Friday.

tkittel avatar Oct 19 '22 12:10 tkittel

I just pushed a few updates, mainly to use a new py API method I have just added in NCrystal 3.4.1, which helps getting the basic material composition in Python even for complicated (multiphase/enriched/...) materials. In the rare and somewhat academic case of materials where the same element now appears both as a natural element and as a explicit isotope, we are now using the openmc natural abundancies to break apart the natural element. This is like what we have been doing all along in Geant4 (with the C++ API), where we instead of course were asking Geant4 to provide the natural abundancies.

I also added material_id and name optional parameters on the Material.from_ncrystal method, which then simply pass these two parameters on to the Material constructor.

tkittel avatar Oct 20 '22 11:10 tkittel

Outstanding issues are tests, documentation and CI as requested. We will have a look at those tomorrow.

tkittel avatar Oct 20 '22 11:10 tkittel

@paulromano I added a tools/ci/gha-install-ncrystal.sh script now (can most likely be improved a bit in the future). Should we also update the gha-install.sh and gha-install.py files? Also, will it be possible to enable the workflows on the current PR, so we can try to get the new tests for NCrystal to work?

Update: I guess there will be a few more files to update, like perhaps gha-script.sh, ci.yml, ...). Do you prefer to take care of that, or should we try?

tkittel avatar Oct 24 '22 09:10 tkittel

Hi, sorry for the delay. I added documentation and tests. @tkittel check the documentation: is brief but more can be found in the NCrystal repo and papers (I put links). I added one simple test to check it works. It will be easy to add tests with more complex config strings.

marquezj avatar Nov 17 '22 15:11 marquezj

@marquezj it looks fine to me. Once we hopefully have some even better NCrystal documentation, we can go back to the openmc documentation and see if we need to update various links, or add new links. But with the current pile of items to do, it might not happen for quite some time :grin:

tkittel avatar Nov 17 '22 18:11 tkittel

So the tests are unsurprisingly failing with:

 import NCrystal
 ModuleNotFoundError: No module named 'NCrystal'

@paulromano could you kindly suggest where you would like changes, to enable NCrystal being built? (btw. NCrystal is on conda-forge now, but I guess this CI setup is not conda-based?).

tkittel avatar Nov 18 '22 10:11 tkittel

@paulromano : Is it normal that some of the checks do not start? All the started ones have now passed. Are there any remaining things we need to look at, or are you happy to merge this now?

Cheers, Thomas

tkittel avatar Dec 01 '22 10:12 tkittel

@tkittel No, you don't need to worry about that. GitHub is expecting a set of test configurations that does not include the new ncrystal option, so that's why those extra ones are shown. Once this is merged, we'll just change the list of expected test configurations to include the ones with ncrystal included.

paulromano avatar Dec 06 '22 16:12 paulromano

Thanks @paulromano , all changes committed as requested.

Anything else? Let me know if you want me to do a rebase of the branch on main now, or do any squashing or whatever :-)

tkittel avatar Dec 07 '22 11:12 tkittel

Hi @paulromano

Despite the status saying "1 change requested", I think we did everything. Let us know if you think more is needed :smile:

Cheers, Thomas

tkittel avatar Jan 06 '23 12:01 tkittel

@tkittel Apologies this fell off my radar -- will take another look soon!

paulromano avatar Jan 09 '23 01:01 paulromano

@tkittel One final question -- should ncrystal_max_energy be configurable? Right now it defaults to 5 eV and it doesn't appear that there's any way to change this.

paulromano avatar Jan 10 '23 12:01 paulromano

I would postpone that for another time, since there is a common issue here with e.g. the Geant4 hooks. Once we eventually look into it, it might anyway be that we go for another solution altogether for this threshold (e.g. a per-material threshold rather than a global one).

tkittel avatar Jan 10 '23 13:01 tkittel

Ok, in that case let's change it to be a global constant as opposed to a global variable in the settings:: namespace. I'll send a quick PR to your branch updating that momentarily.

paulromano avatar Jan 10 '23 18:01 paulromano

Very exciting to us as well, thanks a lot @paulromano for the support and shepherding! :smile:

tkittel avatar Jan 12 '23 13:01 tkittel