Allow depletion by burnup

[drewejohnson]

Ideally, this could / should mimic the time steps, e.g.

>>> pred = PredictorIntegrator(op, burnup=[0.5, 1.0, 2.0, ...], ...)

where burnup is in MWd/kgHM [kg heavy metal].

Looking over the Integrator, it would actually be a pretty straight forward inclusion. Given some power (or power density) and mass of initial heavy metal (through operator.heavy_metal) we can back calculate the corresponding step sizes in days.

If we want to add a regression test, it would probably benefit us to use the same HDF result file and specify the burnup steps to exactly match the current time steps. But these depletion tests are not super short...

I'm open to suggestions on testing this

[rockfool]

Good idea. I agree to add this kind of input parameters for burnup steps.
Serpent also supports "bustep", "butot", "daystep", "daytot", "decstep", "dectot", which mean that:

1. burnup interval
2. cumulative burnup
3. time step interval
4. cumulative time
5. pure decay time interval
6. pure decay cumulative time
   BTW, no further comment on the regression test. It seems easier to use the same results file. ;-)

[drewejohnson]

Agree, it would be really nice to support those settings as well. I think ultimately, the Integrators should have two vectors: time [s] and power. Mainly because the Integrator.__call__(...) methods are already configured to accept a step size in seconds and step power. If we can compute those quantities during Integrator.__init__, then we don't have to compute those every step.

This also eases our testing, because we can throw a myriad of inputs at the __init__ methods, and check their time and power vectors rather than run a full regression test [OpenMC + depletion]. The assumption here is that all the solutions would be the same if you had the same time and power vectors.

One thing worth noting, is that some additional changes may be required for decay. Technically, we could simulate a decay step with a zero power step to compute zero particle-induced reaction rates [fission, (n,2n), etc.]. But this would require running a useless transport simulation only to throw the reaction rates away. Maybe modifying

openmc/openmc/deplete/abc.py

Lines 739 to 743 in 50a271b

	for i, (dt, p) in enumerate(self):
	if i > 0 or self.operator.prev_res is None:
	conc, res = self._get_bos_data_from_operator(i, p, conc)
	else:
	conc, res = self._get_bos_data_from_restart(i, p, conc)

if the step power is zero?

[paulromano]

I'm planning on implementing support for indicating steps by burnup rather than by time. My proposal is as follows:

• Have Integrator accept (operator, timesteps, ..., timestep_units='s')
• The user specifies either 's' or 'MWd/kg' for timestep_units.
• Decay timesteps are indicated by power=0
• Generally, timesteps is just a plain array/list with numeric values. However, in the case where a user wants to specify timesteps by burnup and have decay timesteps, we can support this by having timesteps by an iterable of (value, unit) tuples, e.g. timesteps=[(3.0, 'MWd/kg'), (3.0, 'MWd/kg'), (10*86400, 's')] would mean two 3 MWd/kg timesteps and then a 10 day decay period (assuming the power associated with the third step was zero)
• For convenience, we might want to support alternate time units (day, year, etc.) in timestep_units

@drewejohnson @rockfool @pshriwise Let me know if any of you have objections or can think of a better way of supporting all the different cases.

[rockfool]

I'm planning on implementing support for indicating steps by burnup rather than by time. My proposal is as follows:

• Have Integrator accept (operator, timesteps, ..., timestep_units='s')
• The user specifies either 's' or 'MWd/kg' for timestep_units.
• Decay timesteps are indicated by power=0
• Generally, timesteps is just a plain array/list with numeric values. However, in the case where a user wants to specify timesteps by burnup and have decay timesteps, we can support this by having timesteps by an iterable of (value, unit) tuples, e.g. timesteps=[(3.0, 'MWd/kg'), (3.0, 'MWd/kg'), (10*86400, 's')] would mean two 3 MWd/kg timesteps and then a 10 day decay period (assuming the power associated with the third step was zero)
• For convenience, we might want to support alternate time units (day, year, etc.) in timestep_units

@drewejohnson @rockfool @pshriwise Let me know if any of you have objections or can think of a better way of supporting all the different cases.

Agree. That would be much flexible. Just one idea to add kind of cumulative-timesteps to avoid the calc of time interval sometimes.

[paulromano]

No need to reinvent the wheel -- we should just point users to numpy.cumsum and numpy.diff if they need to convert between cumulative / non-cumulative.

[drewejohnson]

This sounds great; I really like the ability to support multiple units (days and years) 💯