Add neutron elastic scattering physics

[whokion]

Add neutron elastic scattering physics

• Implement neutron-nucleus elastic scattering process, model, interactor and associated classes [enhancement, physics]
• Add a test for neutron-nucleus elastic scattering model and interactor
• Update import, process build, etc. for the newly added neutron elastic process and model

Description

The MR includes the first implementation of the neutron physics (neutron elastic scattering) based on the Geant4's CHIPS model as a part of add neutron physics in the medium-high energy range (up to 20GeV). The implementation is a minimal
without generating secondary ions, but features should be good enough for testing neutron transport in typical HEP user cases. A minimal verification (cross section, scattered energy spectrum of neutrons and direction) has been done with the default Geant4 physics list, but more rigorous validations and verification are required once a test transport application is added.

[whokion]

Sorry for a big set of commits! Please take your time as this MR should be totally decoupled from other tasks. Thanks.

[sethrj]

A few initial comments... this looks great overall!

[amandalund]

This looks great @whokion! Still working through it, but here are a few initial comments.

[amandalund]

Since this is essentially identical to the Livermore macro xs calculator, I wonder if there's a way we could use the same calculator for both and reduce the duplicate code...

[whokion]

Good idea as I will also add a similar set of calculators for NeutronInelastic. I think that a templated class, something like,

template<class MicroXsT, typename MacroXsUnitsT>
class MacroXsCalculator

may be used as a common method. As a quick look, this requires changes in the associated MicroXsCalculator (argument for constructor and operator) and a creation of a specific MicroXsCalculator before hand which imposes an explicit rule for MacroXsCalculator. Please let me know if you have the better idea or if a potential consolidation may be a separated MR.

[sethrj]

I like that idea. I think it would be better to create a new PR that refactors LivermorePEMacroXsCalculator to something like:

template<class F>
class MacroXsCalculator
{
    using ParamsRef = typename F::ParamsRef;
    using Energy = typename F::Energy;
    using XsUnits = units::Native;

    MacroXsCalculator(ParamsRef const&, MaterialView const& material);
    real_type operator()(Energy) const;
};

I think also LivermorePEMicroXsCalculator should be returning either a units::Barn quantity or a native (len^2) unit.

[amandalund]

Suggested change

	auto calc_xs = EPlusGGMacroXsCalculator(
	params_.hardwired.eplusgg_data, material);
	auto calc_xs = NeutronElasticMacroXsCalculator(
	params_.hardwired.chips_data, material);

[sethrj]

Maybe it's time to add kernels to replace these hard coded xs bits?

[sethrj]

@whokion I'm pausing my review because I think there's a lot of naming updates that need to be done, and it seems like you might have addressed some commits but not pushed.

[sethrj]

Forgot to add that close is redundant with the destructor (we're at the end of the block so it's going out of scope. The only reason to call it independently is to be able to check bool(infile) afterward to make sure it closed correctly, which frankly doesn't matter.

Suggested change

infile.close();

[whokion]

Removed!

[sethrj]

The action IDs no longer need to be stored in device data like this (they could be instead stored as member data on the Model class itself). The only place it's used is in the Model where the action ID is returned (and indirectly where make_action_track_executor is used).

[whokion]

As all other models/data have this ActionId consistently, It may be better to refactor all of them together with a separate MR (?).

[sethrj]

You should be able to remove this and just use the typedef in geocel/Types.hh.

[whokion]

Done.

[sethrj]

Suggested change

	CELER_ENSURE(v2 < real_type{1});
	CELER_EXPECT(v2 < real_type{1});

also vsq or v_sq ;)

[whokion]

Changed as suggested.

[sethrj]

The outer call to make unit vector is unnecessary since rotating a unit vector will result in a unit vector

[whokion]

Okay. Thanks - forget to remove the outer one when I added the inner one. Done.

[sethrj]

Is the result MevMomentum? It should definitely be a Quantity then. I think you can also use quantity to handle the transfer to and from GeV internally...

[whokion]

Actually the result is $({\rm MevMomentum})^{2}$.

[sethrj]

More comments now that I've had a chance to dig down into the guts of the momentum transfer sampler 😬

[sethrj]

How hard would it be to preprocess the data so that its units are consistent with the rest of the code? It seems a little odd to divide by 1000 when entering and then multiply by 1000 when leaving. If we don't have to do the extra division and multiplication, we can also make stronger guarantees about the resulting value and skip the extra clamping.

[whokion]

Totally agree! Not sure why the parameterization was done in GeV. Nevertheless, I am not sure whether it is easy to re-parameterize those coefficients by MeV in order to remove this extra conversion. For the moment, I would like to keep it as it is until all validations are done, then we may refactor it later.

[sethrj]

Yeah, it looks like there are a lot of odd exponents and factors...

[sethrj]

On GPU it might be better to not precalculate these: we should do a test compile of this kernel and see if the register pressure here is high. (That can be a follow-on "optimization" check, just something to keep in mind.)

[whokion]

@sethrj Trying to understand those failures 1) Are tests with float using different seeds or are those failures due to something else? 2) in the last test (clhep, ... ubuntu-cuda), values are off by a factor 10 which I do not understand. Have any idea? Thanks.

[sethrj]

@whokion the single-precision failures are because the RNG stream changes and everything is super sensitive. if you look at the test/celeritas/CMakeLists.txt you'll see that most of the EM tests have ${_needs_double} as part of the test call.. that'll disable the test for double precision.

For the failing clhep test, you need to do something like

        macro_xs.push_back(
            native_value_to<units::InvCmXs>(calc_macro_xs(MevEnergy{e})).value());

[whokion]

@sethrj, @amandalund Thank you for all your review comments! I hope that most of suggestions are reflected now and the MR may be good enough to be merged - I will continue to improve codes as I find more documentation related to the model and do more tests along the way. Please let me know any major things that I missed to address. Thanks again!

[sethrj]

I still feel icky about "here's some code that reproduces geant4 but we have no idea what it means, what it does, or where it's from", but I guess that's a good enough starting point for reference and comparison?

[whokion]

I still feel icky about "here's some code that reproduces geant4 but we have no idea what it means, what it does, or where it's from", but I guess that's a good enough starting point for reference and comparison?

I fully agree! On the other hand, we actually know what each component does and how things are related even though detail parameterizations (usually with respect to a certain data set, which sometime even have uncertainties or incompleteness) are not known or fully justified as the nature of NP would not reveal underlying physics cleanly. Anyway, the best standard candle may be physics validation with respect to data besides Geant4 and other MCs results. I will search more documentation on CHIPS or try to contact Kossov, if reachable, how the chips parameterizations were done and put more documentation into the code. After all, details may not be sometimes digestible, but still good to have solid background information on them. Otherwise, we should write a new simulation model based on the celeriton theory.

[amandalund]

Great work @whokion, thanks for our first neutron model! I agree there are still some parts that could be clarified/improved, but it makes sense to wait until after we've done some initial validation.

[amandalund]

Suggested change

	//! A-dependent coefficients for the momentum transfer) of the CHIPS model
	//! A-dependent coefficients for the momentum transfer of the CHIPS model

[whokion]

Fixed.

[amandalund]

Suggested change

// cutoff_params_ = {};

[whokion]

Fixed.

[amandalund]

Since these cross sections are so small, I think even EXPECT_SOFT_EQ(0, calc_micro_xs(el_id)) will pass. We could do something like convert these to barns before checking.

[whokion]

Okay. I will consider this when I consolidate MacroXsCalculator where we probably use a common unit convention for all XsCalculator (next MR).

[amandalund]

Not sure how much it will matter, but it looks like the limit in Geant4 is only approximately 14 MeV/c (actually 13.568559 MeV/c, or log of momentum in GeV/c = -4.3).

[whokion]

Okay. I changed it to 13.568559 even though it should be a rough value anyway.

[whokion]

@amandalund Thanks for additional review comments!