CRAM clip negative atoms

[yrrepy]

Description

CRAM does not guarantee positivity of nuclide concentrations.
When atom densities get very low, the solver can start to oscillate about zero.

This PR addresses that by clipping to zero any negative value that the solver produces. This then propagates to the next step and the depletion_results.h5 file.
No impact is expected on Bateman solutions as this behaviour arises when concentrations are already infinitesimally low.

The new behaviour is made opt-in with clip_negative_densities=True

The attached sample script demonstrates the negative solutions.
demo.py

No clipping:

[openmc.deplete] t=0.0 s, dt=1.0 s, source=100000000000000.0
[openmc.deplete] t=1.0 s, dt=43200.0 s, source=0.0
/home/perry/py313/lib64/python3.13/site-packages/uncertainties/core.py:1024: UserWarning: Using UFloat objects with std_dev==0 may give unexpected results.
  warn("Using UFloat objects with std_dev==0 may give unexpected results.")
[openmc.deplete] t=43201.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=86401.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=129601.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=172801.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=216001.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=259201.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=302401.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=345601.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=388801.0 s, dt=43200.0 s, source=0.0

WITHOUT clipping — Co62 oscillates around zero:
  step  1 (t=   0h):   5.79e+12
  step  2 (t=  12h):  -1.25e-34  *** NEGATIVE ***
  step  3 (t=  24h):   2.70e-81
  step  4 (t=  36h): -5.85e-128  *** NEGATIVE ***
  step  5 (t=  48h):  1.26e-174
  step  6 (t=  60h): -2.73e-221  *** NEGATIVE ***
  step  7 (t=  72h):  5.91e-268
  step  8 (t=  84h): -1.28e-314  *** NEGATIVE ***
  step  9 (t=  96h):   0.00e+00
  step 10 (t= 108h):   0.00e+00
  step 11 (t= 120h):   0.00e+00

4 negative values out of 12 steps (without clipping)

With Clipping

[openmc.deplete] t=0.0 s, dt=1.0 s, source=100000000000000.0
[openmc.deplete] t=1.0 s, dt=43200.0 s, source=0.0
/home/perry/py313/lib64/python3.13/site-packages/openmc/deplete/abc.py:754: UserWarning: Clipping negative atom densities (worst: Co62 = -1.2509e-34)
  warn(f"Clipping negative atom densities "
[openmc.deplete] t=43201.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=86401.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=129601.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=172801.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=216001.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=259201.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=302401.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=345601.0 s, dt=43200.0 s, source=0.0
[openmc.deplete] t=388801.0 s, dt=43200.0 s, source=0.0

WITH clipping — negatives clipped to zero:
  step  1 (t=   0h):   5.79e+12
  step  2 (t=  12h):   0.00e+00
  step  3 (t=  24h):   0.00e+00
  step  4 (t=  36h):   0.00e+00
  step  5 (t=  48h):   0.00e+00
  step  6 (t=  60h):   0.00e+00
  step  7 (t=  72h):   0.00e+00
  step  8 (t=  84h):   0.00e+00
  step  9 (t=  96h):   0.00e+00
  step 10 (t= 108h):   0.00e+00
  step 11 (t= 120h):   0.00e+00

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) 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)

[yrrepy]

On further reflection, this may be the ideal location to implement a minimum number density filter (clip).

Here is an augmented method
https://github.com/yrrepy/openmc/tree/CRAM_clip-min-atom-density

  integrator = openmc.deplete.PredictorIntegrator(op, timesteps, source_rates,
                                                  clip_min_atom_density=1e-20)

    def _timed_deplete(self, n, rates, dt, i=None, matrix_func=None):
        start = time.time()
        results = deplete(
            self._solver, self.chain, n, rates, dt, i, matrix_func,
            self.transfer_rates, self.external_source_rates)
        if self.clip_min_atom_density is not None:
            volumes = self.operator.number.volume
            for j, r in enumerate(results):
                threshold = self.clip_min_atom_density * 1e24 * volumes[j]
                if not self._warned_negative_density:
                    neg_mask = r < 0
                    if neg_mask.any():
                        idx = np.argmin(r)
                        nuc_list = list(
                            self.operator.number.burnable_nuclides)
                        nuc = nuc_list[idx]
                        density = r[idx] / (1e24 * volumes[j])
                        warn(f"Clipping {neg_mask.sum()} negative atom "
                             f"densities (worst: {nuc} = "
                             f"{density:.4e} atoms/b-cm)")
                        self._warned_negative_density = True
                if not self._warned_sub_threshold and threshold > 0:
                    sub_mask = (r >= 0) & (r < threshold)
                    if sub_mask.any():
                        warn(f"Clipping {sub_mask.sum()} sub-threshold "
                             f"atom densities (threshold: "
                             f"{self.clip_min_atom_density:.4e} "
                             f"atoms/b-cm)")
                        self._warned_sub_threshold = True
                results[j] = np.where(r < threshold, 0.0, r)
        return time.time() - start, results

demo.py

The other alternative is to perform the clipping after integration is complete (after openmc deplete finalizes all time-steps). I already have a workflow (processing script) that does this but it is rather cumbersome as you have to batch-wise sweep through the depletion_results.h5 file to do it effectively.

[drewejohnson]

Trimming negatives makes a lot of sense. They're non-physical, but they can be an artifact of some CRAM workflows (decay or low flux, insufficient order)

[yrrepy]

Another examples of CRAM clipping:
https://github.com/jlanversin/ONIX/blob/master/onix/salameche/cram.py

in-case it inspires

another possible location for clipping would be on step_result.py

openmc.deplete.PredictorIntegrator(op, dt, source_rates, clip_min_atom_density=1e-20)

  # Add data — zero negligible densities before archival                                                                                                                   
          inds = [self.mat_to_hdf5_ind[mat] for mat in self.index_mat]                                                                                                   
          low = min(inds)                                                                                                                                                  
          high = max(inds)
          for i, mat in enumerate(self.index_mat):                                                                                                                         
              threshold = clip_min_atom_density * self.volume[mat] * 1e24                                                                                                                  
              self.data[i, self.data[i] < threshold] = 0.0                                                                                                                 
          number_dset[index, low:high+1] = self.data                                                                                                                       
          if has_reactions:                                                                                                                                              
              rxn_dset[index, low:high+1] = self.rates                                                                                                                     
          if comm.rank == 0:                                                                                                                                             
              eigenvalues_dset[index] = self.k                                                                                                                             
              time_dset[index] = self.time
              source_rate_dset[index] = self.source_rate                                                                                                                   
              if self.proc_time is not None:                                                                                                                             
                  proc_time_dset[index] = (                                                                                                                                
                      self.proc_time / (comm.size * self.n_hdf5_mats)                                                                                                      
                  )   
                  

This has the advantage of allowing CRAM to do it's numerical thing undisturbed.
But in rare situations you could get oscillations making it through to the depletion_results.h5 file, with the positives "blinking" in and out between zeroes, if clip_min_atom_density=0, for instance in the above examples