openmc-dev / openmc

OpenMC Monte Carlo Code
https://docs.openmc.org
Other
699 stars 444 forks source link

Allow pure decay IndependentOperator #2966

Closed gridley closed 1 week ago

gridley commented 2 weeks ago

Description

Fixes #2963

Fixes #2968

EDIT this also fixes #2968 with the second commit it adds.

Autopep8'd the independent operator file, and adds handling for decay only steps as we have in the coupled operator. This prevents an error where the heating is calculated as zero and an error is thrown erroneously. After taking this PR, the following code can produce this nice little plot in which we track the first two decay progeny of U239. It has a half life around 1400s and we can see that behavior nicely in this plot. Notably I need to turn off multiprocessing to have this work, described below.

It's so nice to have standalone decay like this, thanks @paulromano for implementing the IndependentOperator 🙂

import numpy as np
import openmc
import openmc.deplete
import matplotlib.pyplot as plt

# TODO I need help on why this does not work
openmc.deplete.USE_MULTIPROCESSING = False

chain_file = '/Users/gavin/Code/chain_endfb80_pwr.xml'
openmc.config['chain_file'] = chain_file
openmc.config['cross_sections'] = '/Users/gavin/Code/endfb-viii.0-hdf5/cross_sections.xml'

decay_times = np.ones(100) * 50.0

# Create one material to deplete. It is ten grams of U239.
mat = openmc.Material()
mat.name = 'I decay!'
mat.add_nuclide('U239', 1.0, 'ao')
mat.volume = 10.0
mat.set_density('g/cc', 1.0)

mats = openmc.Materials([mat])

empty_micro = openmc.deplete.MicroXS(np.array([]).reshape((0, 0)), [], [])
op = openmc.deplete.IndependentOperator(mats, [1.0], [empty_micro], chain_file)
integ = openmc.deplete.PredictorIntegrator(op, decay_times, power=0.0, timestep_units='s')
integ.integrate()

results = openmc.deplete.Results()

atoms_u239 = results.get_mass(mat, 'U239')
atoms_pu239 = results.get_mass(mat, 'Pu239')
atoms_np239 = results.get_mass(mat, 'Np239')
plt.semilogy(atoms_u239[0], atoms_u239[1])
plt.semilogy(atoms_np239[0], atoms_np239[1])
plt.semilogy(atoms_pu239[0], atoms_pu239[1])
plt.legend(('U239', 'Np239', 'Pu239'))
plt.xlabel('Time (s)')
plt.ylabel('Mass (g)')
plt.ylim([0.01, 11.0])
plt.show()
image

HOWEVER I am running into issues with multiprocessing set to true and may need some help. This error recurrently appears if depletion multiprocessing is turned on:

    super().__init__(process_obj)
  File "/opt/homebrew/Cellar/python@3.12/3.12.2_1/Frameworks/Python.framework/Versions/3.12/lib/python3.12/multiprocessing/popen_fork.py", line 19, in __init__
    self._launch(process_obj)
  File "/opt/homebrew/Cellar/python@3.12/3.12.2_1/Frameworks/Python.framework/Versions/3.12/lib/python3.12/multiprocessing/popen_spawn_posix.py", line 42, in _launch
    prep_data = spawn.get_preparation_data(process_obj._name)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/Cellar/python@3.12/3.12.2_1/Frameworks/Python.framework/Versions/3.12/lib/python3.12/multiprocessing/spawn.py", line 164, in get_preparation_data
    _check_not_importing_main()
  File "/opt/homebrew/Cellar/python@3.12/3.12.2_1/Frameworks/Python.framework/Versions/3.12/lib/python3.12/multiprocessing/spawn.py", line 140, in _check_not_importing_main
    raise RuntimeError('''
RuntimeError: 
        An attempt has been made to start a new process before the
        current process has finished its bootstrapping phase.

        This probably means that you are not using fork to start your
        child processes and you have forgotten to use the proper idiom
        in the main module:

            if __name__ == '__main__':
                freeze_support()
                ...

        The "freeze_support()" line can be omitted if the program
        is not going to be frozen to produce an executable.

        To fix this issue, refer to the "Safe importing of main module"
        section in https://docs.python.org/3/library/multiprocessing.html

I guess we've missed some idiom here and haven't hit this race condition before because this case runs the depletion solver without much latency? If anyone knows where this freeze_support might reasonably go, I would very much appreciate your input.

Checklist

gridley commented 2 weeks ago

Note: the multiprocessing error only seems to happen on mac. It's working fine on my ubuntu desktop.

gridley commented 1 week ago

Oh sick, thank you Olek! My apologies, unsure why I thought that. It is a lovely feature to have for analyzing reactors.

Will get the test in fosho.

gridley commented 1 week ago

OK @yardasol, I found that the test I've added here was able to reproduce the original problem, and passes with the introduction of these commits.

paulromano commented 1 week ago

@yardasol are you satisfied with the changes in this PR?

gridley commented 1 week ago

I would really prefer to not do that. The only potential problem with zero power is hitting that branch between calculating a normalizing factor and not, so I can't think of any difference that would take place with a filled MicroXS object.

yardasol commented 1 week ago

If @paulromano thinks it's fine then I'm satisfied.