Fix caching issue when using NCrystal materials

[tkittel]

Description

This small but important PR fixes a caching issue affecting NCrystal materials when any of the isotopes of NCrystal materials also appear in other materials with the same temperature. It seems to be simply that the "is the cached xs still valid" check never actually took the NCrystal contribution into account.

This issue was observed by two different people recently, and was initially reported in the ncrystal issue tracker at mctools/ncrystal#243.

Fixes #3540

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)

[marquezj]

I think the update of the NCrystal XS could be taken out of the if condition. Then the ACE XS is not recalculated if only the NCrystal XS changes, but the XS is updated anyway:

void Particle::update_neutron_xs(
  int i_nuclide, int i_grid, int i_sab, double sab_frac, double ncrystal_xs)
{
  // Get microscopic cross section cache
  auto& micro = this->neutron_xs(i_nuclide);

  // If the cache doesn't match, recalculate micro xs
  bool recalculate_xs = ( this->E() != micro.last_E || 
                          this->sqrtkT() != micro.last_sqrtkT ||
                          i_sab != micro.index_sab || 
                          sab_frac != micro.sab_frac );
  if ( recalculate_xs ) {
    data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this);
  }
  // If NCrystal is being used, update micro cross section cache
  if ((recalculate_xs || ncrystal_xs != micro.last_ncrystal_xs ) 
       && (ncrystal_xs >= 0.0) ) {
    ncrystal_update_micro(ncrystal_xs, micro);
    micro.last_ncrystal_xs = ncrystal_xs;
  }
}

[tkittel]

That is getting a bit too complicated for my liking, I can't really wrap my head around whether or not that will always end up correct.

[marquezj]

Ok, let's keep it your way.

[marquezj]

The crossSectionNonOriented() method is now deprecated in NCrystal.

[marquezj]

This fixes a bug where a neutron flying crossed from a cell filled with a material A to a NCrystal material B, but the XS was not updated because the last interaction was with the same nuclide, energy and temperature. Now the XS is updated if the energy changes, the temperature changes or the NCrystal XS changes. Note that the NCrystal XS is always updated for NCrystal materials if the neutron energy is below NCRYSTAL_MAX_ENERGY.

[tkittel]

All looks fine for merging from my POV.

[paulromano]

Just putting a placeholder review here as I'd like to take a look at this before we merge

[marquezj]

Hi @paulromano. This is a MWE that shows the issue with two materials (one provided by NCrystal) that share a nuclide at the same temperature. We can open an issue here if you consider it necessary to keep track of which bug we are fixing.

mwe.txt

[paulromano]

Thanks @marquezj. The issue makes sense to me -- I'm just slightly wary of adding new members on the NuclideMicroXS struct because performance can be quite sensitive to it so I'm trying to figure out if there's a way we can avoid that.

[marquezj]

@paulromano: it might not be necessary to add another member to NuclideMicroXS. The last NCrystal XS will be stored in micro.elastic, and then we could use the following condition:

  // If the cache doesn't match, recalculate micro xs
  if (this->E() != micro.last_E || this->sqrtkT() != micro.last_sqrtkT ||
      i_sab != micro.index_sab || sab_frac != micro.sab_frac ||
      ncrystal_xs != micro.elastic) {
    data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this);

    // If NCrystal is being used, update micro cross section cache
    if (ncrystal_xs >= 0.0) {
      data::nuclides[i_nuclide]->calculate_elastic_xs(*this);
      ncrystal_update_micro(ncrystal_xs, micro);
    }
  }

[tkittel]

@marquezj But how would you know if the last stored elastic xs is coming from NCrystal, and not an elastic process in another material?

[marquezj]

@tkittel I think it does not matter. Let's say micro.elastic and ncrystal_xs have the same value.

• The material has to be from NCrystal, otherwise ncrystal_xs will be -1.0. And this will never be equal to micro.elastic.
• Then two things could happen:

1. One of the ACE recalculation conditions is triggered and the XS is recomputed and updated, or
2. None is triggered, which means that all other XS are also the same and then not need to be updated.

[tkittel]

@marquezj Hmm... maybe I am confused, but I am not sure I get your first bullet point. If ncrystal_xs is 5.0 and micro.elastic is also 5.0, then you know that the current material is an NCrystal material, but you can't know that the 5.0 in micro elastic was not due to some free-gas nuclei in a previous material (i.e. non-NCrystal). Or what am I missing?

[marquezj]

Why do we need to know if the previous material was from NCrystal or not?

All we need to know is if we need to update the XS or not. If the previous material was non-NCrystal and the elastic XS was 5 b, it still is 5 b (because the other ACE conditions are not satisfied). And calling ncrystal_update_micro() to add and substract 5 b will be unnecessary.

Am I missing something?

[tkittel]

@marquezj probably not, it all just seems a tad too fragile for my liking :-)

But if you and @paulromano are convinced it makes sense, then ok with me, go ahead and do the changes you envision.

[tkittel]

So I have a feeling @paulromano is waiting for the updates from @marquezj while @marquezj is waiting for comments from @paulromano ;-)

Could you guys perhaps clarify who is waiting for who and get this moving again? :-)

[paulromano]

Thanks for the nudge @tkittel. I'll take a shot at it!

[paulromano]

Well, after trying out different ideas, I couldn't come up with a good solution that didn't involve the addition of a new data member on NuclideMicroXS. Even though ncrystal_xs is the same as elastic when an NCrystal cross section is present, the problem is that it always forces a recalculation of the cross section when ncrystal_xs < 0.0 (even when it wouldn't otherwise be necessary). So, we'll go ahead with the original approach here and can always re-evaluate later.