Fix caching issue when using NCrystal materials (#3538)

Co-authored-by: Jose Ignacio Marquez Damian <22483345+marquezj@users.noreply.github.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

Author: 3 people

include/openmc/particle_data.h

@@ -154,9 +154,10 @@ struct NuclideMicroXS {
 
   // Energy and temperature last used to evaluate these cross sections.  If
   // these values have changed, then the cross sections must be re-evaluated.
-  double last_E {0.0};      //!< Last evaluated energy
-  double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant
-                            //!< * temperature (eV))
+  double last_E {0.0};       //!< Last evaluated energy
+  double last_sqrtkT {0.0};  //!< Last temperature in sqrt(Boltzmann constant
+                             //!< * temperature (eV))
+  double ncrystal_xs {-1.0}; //!< NCrystal cross section
 };
 
 //==============================================================================

openmc/plotter.py

@@ -501,7 +501,7 @@ def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None,
                     elif ncrystal_cfg:
                         import NCrystal
                         nc_scatter = NCrystal.createScatter(ncrystal_cfg)
-                        nc_func = nc_scatter.crossSectionNonOriented
+                        nc_func = nc_scatter.xsect
                         nc_emax = 5 # eV # this should be obtained from NCRYSTAL_MAX_ENERGY
                         energy_grid = np.union1d(np.geomspace(min(energy_grid),
                                                               1.1*nc_emax,

src/particle.cpp

@@ -855,10 +855,12 @@ void Particle::update_neutron_xs(
 
   // 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) {
+      i_sab != micro.index_sab || sab_frac != micro.sab_frac ||
+      ncrystal_xs != micro.ncrystal_xs) {
     data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this);
 
     // If NCrystal is being used, update micro cross section cache
+    micro.ncrystal_xs = ncrystal_xs;
     if (ncrystal_xs >= 0.0) {
       data::nuclides[i_nuclide]->calculate_elastic_xs(*this);
       ncrystal_update_micro(ncrystal_xs, micro);