Initial implementation of fixed-source activation calculations

[paulromano]

This pull request implements initial support for fixed-source activation calculations. The Integrator classes now accept a source_rates argument that corresponds to the absolute source rates in neutrons/sec that can be used in lieu of the power or power_density arguments. Additionally, the energy_mode argument of Operator() was changed to normalization_mode and a new "source-rate" mode was introduced for this purpose. All internal function calls that accepted power argument now accept a source_rate argument that represents either a power in [W] or a source rate in [neutrons/sec].

A few implementation notes:

• The abstract EnergyHelper class is now NormalizationHelper. There is a new intermediate EnergyNormalizationHelper class that is the parent of ChainFissionHelper and EnergyScoreHelper.
• To avoid all the repetitive docstrings for integrator classes (which are easy to mess up if you need to change), I've done a little metaprogramming hack to ensure they are consistent. Each integrator class now has an @add_params class decorator that modifies the docstring, adding the "Parameters" section that is shared among the classes.
• The Nuclide class that is used in Chain was enhanced to have add_reaction and add_decay_mode methods. This makes it easy to construct chains "by hand", if you wish.
• Before, the deplete module accessed internal tally results directly, which are not normalized by the number of tally realizations. This is actually ok in k-eigenvalue depletion calculations because the absolute numbers are just renormalized based on the power, but in the case of fixed-source activation calculations, it's important that we get the properly batch-normalized values. Thus, all references to results now instead rely on mean which will divide by the number of realizations.

Lastly, I've added a new "unit" test for the activation capability. It uses a simple one-nuclide model, and what it does is 1) run a transport calculation to determine the (n,gamma) reaction rate, 2) calculate the necessary source rate to deplete the nuclide concentration by half over a given time (chosen at random), 3) runs the activation calculation, and 4) confirms that the nuclide was indeed depleted by half.

Some caveats to note:

• Activation steps with a zero source rate will still result in an unnecessary transport calculation being run;
• When using the integrator classes, a transport calculation is always run at the end of a step, which may be undesirable for activation calculations.

[makeclean]

Exciting @paulromano I look forward to giving this a go :)

[pshriwise]

What @makeclean said! I'm planning on trying it out soon as well and will review the code too.

[drewejohnson]

I am a huge fan of the add_params decorator and the new chain and nuclide methods. There are some minor in-line comments, and I have some follow on questions.

From the new source_rates documentation on the Operator and how you added this in the new unit test, it appears that we prefer source_rates be passed as an iterable of floats. The uses for power and power density support passing a single float, and it looks like that would also work for source_rates. I don't have experience with activation analysis, so it very well could be that having varied source rates like pulses is the most common situation, as opposed to a constant power / power density. Do we want to explicitly support a constant source?

And it looks like the depletion reference file has been ever so slightly updated (~300 bytes 😅) More of a curiosity than anything, but do you know what caused it? h5diff gives me zero differences between the files, so it's nothing major

[drewejohnson]

If you're using ravel to just get the very first entry from an array, you can also use

Suggested change

	capture_rate = tally.get_values().ravel()[0]
	capture_rate = tally.get_values().flat[0]

Also, there's no way we get a float from tally.get_values? We have a score, so we'll get back an array, but still only containing a single value?

[paulromano]

Good suggestion. To be clearer, I'll change this to tally.mean.flat[0]. And yes, I don't think get_values ever returns a scalar.

[drewejohnson]

Makes sense to me. Just wanted to point out that the v 0.12.0 docs indicate that a scalar may be returned

[drewejohnson]

Thanks for this awesome feature @paulromano! Not sure if @makeclean, @shimwell, or @pshriwise want to chime in here, but I'm good to merge as is 👍

[paulromano]

Thanks @drewejohnson! Yes, I'd like to get input from @shimwell or @makeclean before we merge since this work was done under a contract with them.

[pshriwise]

This is a really nice extension of the depletion module for this capability. Can't wait to try it out on some fusion problems someday!

One thing I noticed is that we should probably check to make sure a normalization mode in the Operator isn't being combined with an inappropriate power/source rate setting when creating an Integrator. For example, an Operator's SourceRateHelper could be operating on an Integrator's source_rate created using a power or power density.

I gave this a try by modifying the test_depletion_activation case and right now the integrator seems to do this happily. The easiest check I can think of off-hand is to check the type of the Operator's _normalization_helper when creating the Integrator. Maybe I'm missing something here though.

Other than that, just a question about the depletion_file specification now that the source rate can have different units and a couple of super small line comments.

[pshriwise]

Now that the units of the source rate can vary, should we also write an indicator of whether this value is in Watts or neutron/s?

[paulromano]

Right now, the Results class doesn't really have that information. It would have to either be passed from the Integrator (which also doesn't keep track of the units currently) or Results would have to dig into the Operator class to determine what units are being used. Either way doesn't seem great to me at present. I would probably lean toward just leaving it as is, but let me know if you see a cleaner way to do this.

[pshriwise]

Yeah I'm good with leaving it as it is. As you mentioned in our offline convo, it isn't really possible to know everything from the current depletion_results.h5 spec as it stands now.

[shimwell]

Looks good from what I can see with the smartphone, but I'm without a computer for a couple of weeks but will take a proper look when I get back from holiday.

[makeclean]

If you're happy I can drive while you're gone, I'm on leave too - back in a week?

…

________________________________ From: Jonathan Shimwell <notifications@github.com> Sent: Monday, August 10, 2020 1:45:26 PM To: openmc-dev/openmc <openmc@noreply.github.com> Cc: Davis, Andrew <andrew.davis@ukaea.uk>; Mention <mention@noreply.github.com> Subject: Re: [openmc-dev/openmc] Initial implementation of fixed-source activation calculations (#1628) Looks good from what I can see with the smartphone, but I'm without a computer for a couple of weeks but will take a proper look when I get back from holiday. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<#1628 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AASTUSR7STOOO5QQBOWP6TLR77TWNANCNFSM4PURTEPA>.

[shimwell]

If you're happy I can drive while you're gone, I'm on leave too - back in a week?
…
________________________________ From: Jonathan Shimwell notifications@github.com Sent: Monday, August 10, 2020 1:45:26 PM To: openmc-dev/openmc openmc@noreply.github.com Cc: Davis, Andrew andrew.davis@ukaea.uk; Mention mention@noreply.github.com Subject: Re: [openmc-dev/openmc] Initial implementation of fixed-source activation calculations (#1628) Looks good from what I can see with the smartphone, but I'm without a computer for a couple of weeks but will take a proper look when I get back from holiday. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<#1628 (comment)>, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AASTUSR7STOOO5QQBOWP6TLR77TWNANCNFSM4PURTEPA.

Thanks @makeclean that sounds good. Sorry for both being away just when there is stuff to do @paulromano

[paulromano]

No problem! This is vacation season after all.

I have another branch dealing with the extra transport solves for decay-only timesteps (as well as the transport solve at the very last time point). It's only about 10 lines of changes though, so I could add it here. @drewejohnson @pshriwise What do you think?

[drewejohnson]

@paulromano that's excellent to hear on the decay-only solves. I've been kicking it around in my head for a bit

My preference would be to have a unique PR, if anything to have it merged sooner 😅

[shimwell]

Just out of interest, can users accidentally put in negative time steps be put in to do reverse decay

[paulromano]

@shimwell I've never heard of such a thing. There's currently a value check on the timesteps to ensure they are real and positive. I'm not sure what would happen if you allowed a negative timestep. CRAM is designed on the premise that the eigenvalues of the burnup matrix are confined to a small region around the negative real axis. If you were to allow negative timesteps, that may not be true anymore. So, perhaps it can be done but probably not with CRAM is my guess.

[shimwell]

Thanks for this nice feature.

We have been using it a bit and it looks great.

We plan to keep testing and will feedback if needed

Perhaps on day we should chat about how to add gamma information from the decay, but that is for the future

[pshriwise]

Thank you @paulromano! And thanks to @shimwell and @makeclean for doing some of their own testing on these changes.