Fixed a bug in cell last update

include/openmc/particle_data.h

@@ -210,6 +210,8 @@ struct MacroXS {
   double incoherent;      //!< macroscopic incoherent xs
   double photoelectric;   //!< macroscopic photoelectric xs
   double pair_production; //!< macroscopic pair production xs
+
+  double last_E {0.0}; //!< last evaluated energy in [eV]
 };
 
 //==============================================================================

openmc/data/decay.py

@@ -564,7 +564,7 @@ def sources(self):
 _DECAY_PHOTON_ENERGY = {}
 
 
-def decay_photon_energy(nuclide: str) -> Univariate | None:
+def decay_photon_energy(nuclide: str, chain = None) -> Univariate | None:
     """Get photon energy distribution resulting from the decay of a nuclide
 
     This function relies on data stored in a depletion chain. Before calling it
@@ -577,14 +577,25 @@ def decay_photon_energy(nuclide: str) -> Univariate | None:
     ----------
     nuclide : str
         Name of nuclide, e.g., 'Co58'
-
+    chain : Chain, optional.
+            Chain file to get decay energy from
+
     Returns
     -------
     openmc.stats.Univariate or None
         Distribution of energies in [eV] of photons emitted from decay, or None
         if no photon source exists. Note that the probabilities represent
         intensities, given as [Bq/atom] (in other words, decay constants).
     """
+    from openmc.deplete.chain import Chain
+    if chain is not None:
+        cv.check_type('chain', chain, Chain)
+        for nuc in chain.nuclides:
+            if nuc.name == nuclide:
+                return nuc.sources.get('photon')
+        else:
+            return
+
     if not _DECAY_PHOTON_ENERGY:
         chain_file = openmc.config.get('chain_file')
         if chain_file is None:
@@ -593,7 +604,6 @@ def decay_photon_energy(nuclide: str) -> Univariate | None:
                 "openmc.config['chain_file'] in order to load decay data."
             )
 
-        from openmc.deplete import Chain
         chain = Chain.from_xml(chain_file)
         for nuc in chain.nuclides:
             if 'photon' in nuc.sources:
@@ -610,7 +620,7 @@ def decay_photon_energy(nuclide: str) -> Univariate | None:
 _DECAY_ENERGY = {}
 
 
-def decay_energy(nuclide: str):
+def decay_energy(nuclide: str, chain = None):
     """Get decay energy value resulting from the decay of a nuclide
 
     This function relies on data stored in a depletion chain. Before calling it
@@ -623,13 +633,24 @@ def decay_energy(nuclide: str):
     ----------
     nuclide : str
         Name of nuclide, e.g., 'H3'
+    chain : Chain, optional.
+        Chain file to get decay energy from
 
     Returns
     -------
     float
         Decay energy of nuclide in [eV]. If the nuclide is stable, a value of
         0.0 is returned.
     """
+    from openmc.deplete.chain import Chain
+    if chain is not None:
+        cv.check_type('chain', chain, Chain)
+        for nuclide in chain.nuclides:
+            if nuc.name == nuclide and nuc.decay_energy:
+                return nuc.decay_energy
+        else:
+            return 0.0
+
     if not _DECAY_ENERGY:
         chain_file = openmc.config.get('chain_file')
         if chain_file is None:
@@ -638,7 +659,6 @@ def decay_energy(nuclide: str):
                 "openmc.config['chain_file'] in order to load decay data."
             )
 
-        from openmc.deplete import Chain
         chain = Chain.from_xml(chain_file)
         for nuc in chain.nuclides:
             if nuc.decay_energy:

openmc/deplete/r2s.py

@@ -230,7 +230,7 @@ def run(
         )
         self.step3_photon_transport(
             photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
-            mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs
+            mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs, chain_file = chain_file
         )
 
         return output_dir
@@ -405,6 +405,7 @@ def step3_photon_transport(
         output_dir: PathLike = 'photon_transport',
         mat_vol_kwargs: dict | None = None,
         run_kwargs: dict | None = None,
+        chain_file: PathLike | None = None,
     ):
         """Run the photon transport step.
 
@@ -437,12 +438,17 @@ def step3_photon_transport(
         run_kwargs : dict, optional
             Additional keyword arguments passed to :meth:`openmc.Model.run`
             during the photon transport step. By default, output is disabled.
+        chain_file : PathLike, optional
+                    Path to the depletion chain XML file to use during activation. If
+                    not provided, the default configured chain file will be used.
         """
 
         # TODO: Automatically determine bounding box for each cell
         if bounding_boxes is None and self.method == 'cell-based':
             raise ValueError("bounding_boxes must be provided for cell-based "
                              "R2S calculations.")
+        from openmc.deplete import Chain
+        chain = Chain.from_xml(chain_file) if chain_file is not None else None
 
         # Set default run arguments if not provided
         if run_kwargs is None:
@@ -490,7 +496,7 @@ def step3_photon_transport(
             # Create decay photon source
             if self.method == 'mesh-based':
                 self.photon_model.settings.source = \
-                    self.get_decay_photon_source_mesh(time_index)
+                    self.get_decay_photon_source_mesh(time_index, chain_file=chain_file)
             else:
                 sources = []
                 results = self.results['depletion_results']
@@ -510,7 +516,7 @@ def step3_photon_transport(
 
                     # Create decay photon source source
                     space = openmc.stats.Box(*bounding_box)
-                    energy = activated_mat.get_decay_photon_energy()
+                    energy = activated_mat.get_decay_photon_energy(chain=chain)
                     strength = energy.integral() if energy is not None else 0.0
                     source = openmc.IndependentSource(
                         space=space,
@@ -539,7 +545,8 @@ def step3_photon_transport(
 
     def get_decay_photon_source_mesh(
         self,
-        time_index: int = -1
+        time_index: int = -1,
+        chain_file: PathLike | None = None,
     ) -> list[openmc.IndependentSource]:
         """Create decay photon source for a mesh-based calculation.
 
@@ -558,6 +565,9 @@ def get_decay_photon_source_mesh(
         ----------
         time_index : int, optional
             Time index for the decay photon source. Default is -1 (last time).
+        chain_file : PathLike, optional
+                    Path to the depletion chain XML file to use during activation. If
+                    not provided, the default configured chain file will be used.
 
         Returns
         -------
@@ -566,6 +576,9 @@ def get_decay_photon_source_mesh(
             each mesh element-material combination with non-zero source strength.
 
         """
+        from openmc.deplete import Chain
+        chain = Chain.from_xml(chain_file) if chain_file is not None else None
+
         mat_dict = self.neutron_model._get_all_materials()
 
         # List to hold all sources
@@ -607,7 +620,7 @@ def get_decay_photon_source_mesh(
                 activated_mat = results[time_index].get_material(str(original_mat.id))
 
                 # Create decay photon source
-                energy = activated_mat.get_decay_photon_energy()
+                energy = activated_mat.get_decay_photon_energy(chain=chain)
                 if energy is not None:
                     strength = energy.integral()
                     space = openmc.stats.Box(*bbox)

openmc/material.py

@@ -28,7 +28,6 @@
 from openmc.data.function import Tabulated1D
 from openmc.data import mass_energy_absorption_coefficient, dose_coefficients
 
-
 # Units for density supported by OpenMC
 DENSITY_UNITS = ('g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3', 'sum',
                  'macro')
@@ -338,7 +337,8 @@ def get_decay_photon_energy(
         units: str = 'Bq',
         volume: float | None = None,
         exclude_nuclides: list[str] | None = None,
-        include_nuclides: list[str] | None = None
+        include_nuclides: list[str] | None = None,
+        chain: Chain | None = None,
     ) -> Univariate | None:
         r"""Return energy distribution of decay photons from unstable nuclides.
 
@@ -359,6 +359,8 @@ def get_decay_photon_energy(
         include_nuclides : list of str, optional
             Nuclides to include in the photon source calculation. If specified,
             only these nuclides are used.
+        chain : Chain, optional.
+            Chain file to get decay energy from
 
         Returns
         -------
@@ -367,7 +369,9 @@ def get_decay_photon_energy(
             the total intensity of the photon source in the requested units.
 
         """
+        from openmc.deplete import Chain
         cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3'})
+        cv.check_type('chain', chain, Chain, none_ok=True)
 
         if exclude_nuclides is not None and include_nuclides is not None:
             raise ValueError("Cannot specify both exclude_nuclides and include_nuclides")
@@ -391,7 +395,7 @@ def get_decay_photon_energy(
             if include_nuclides is not None and nuc not in include_nuclides:
                 continue
 
-            source_per_atom = openmc.data.decay_photon_energy(nuc)
+            source_per_atom = openmc.data.decay_photon_energy(nuc, chain=chain)
             if source_per_atom is not None and atoms_per_bcm > 0.0:
                 dists.append(source_per_atom)
                 probs.append(1e24 * atoms_per_bcm * multiplier)

src/geometry.cpp

@@ -173,12 +173,21 @@ bool find_cell_inner(
       }
 
       // Set the material, temperature and density multiplier.
-      p.material_last() = p.material();
-      p.material() = c.material(p.cell_instance());
-      p.sqrtkT_last() = p.sqrtkT();
-      p.sqrtkT() = c.sqrtkT(p.cell_instance());
-      p.density_mult_last() = p.density_mult();
-      p.density_mult() = c.density_mult(p.cell_instance());
+      if (p.sqrtkT() != -1) {
+        // if current material is set, set last material
+        p.material_last() = p.material();
+        p.material() = c.material(p.cell_instance());
+        p.sqrtkT_last() = p.sqrtkT();
+        p.sqrtkT() = c.sqrtkT(p.cell_instance());
+        p.density_mult_last() = p.density_mult();
+        p.density_mult() = c.density_mult(p.cell_instance());
+      } else {
+        // set both current data and last data to current data
+        p.material_last() = p.material() = c.material(p.cell_instance());
+        p.sqrtkT_last() = p.sqrtkT() = c.sqrtkT(p.cell_instance());
+        p.density_mult_last() = p.density_mult() =
+          c.density_mult(p.cell_instance());
+      }
 
       return true;
 

src/material.cpp

@@ -817,6 +817,7 @@ void Material::calculate_xs(Particle& p) const
   p.macro_xs().absorption = 0.0;
   p.macro_xs().fission = 0.0;
   p.macro_xs().nu_fission = 0.0;
+  p.macro_xs().last_E = p.E();
 
   if (p.type().is_neutron()) {
     this->calculate_neutron_xs(p);

src/mgxs.cpp

@@ -615,6 +615,7 @@ void Mgxs::calculate_xs(Particle& p)
   p.macro_xs().nu_fission =
     fissionable ? xs[temperature].nu_fission(angle, p.g()) * p.density_mult()
                 : 0.;
+  p.macro_xs().last_E = p.E();
 }
 
 //==============================================================================