diff --git a/.gitignore b/.gitignore index 88e59895d..3b2a24a4d 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,7 @@ *.o *.log *.out +*.pkl # Compiler python objects *.pyc diff --git a/docs/source/methods/depletion.rst b/docs/source/methods/depletion.rst index dce1f2503..1b131bed5 100644 --- a/docs/source/methods/depletion.rst +++ b/docs/source/methods/depletion.rst @@ -103,16 +103,17 @@ integrate over the entire timestep. Our aim here is not to exhaustively describe all integration methods but rather to give a few examples that elucidate the main considerations one must take into account when choosing a method. Generally, there is a tradeoff between the -accuracy of the method and its computational expense. The expense is driven -almost entirely by the time to compute a transport solution, i.e., to evaluate -:math:`\mathbf{A}` for a given :math:`\mathbf{n}`. Thus, the cost of a method -scales with the number of :math:`\mathbf{A}` evaluations that are performed per -timestep. On the other hand, methods that require more evaluations generally -achieve higher accuracy. The predictor method only requires one evaluation and -its error converges as :math:`\mathcal{O}(h)`. The CE/CM method requires two -evaluations and is thus twice as expensive as the predictor method, but achieves -an error of :math:`\mathcal{O}(h^2)`. An exhaustive description of time -integration methods and their merits can be found in the `thesis of Colin Josey +accuracy of the method and its computational expense. In the case of +transport-coupled depletion, the expense is driven almost entirely by the time +to compute a transport solution, i.e., to evaluate :math:`\mathbf{A}` for a +given :math:`\mathbf{n}`. Thus, the cost of a method scales with the number of +:math:`\mathbf{A}` evaluations that are performed per timestep. On the other +hand, methods that require more evaluations generally achieve higher accuracy. +The predictor method only requires one evaluation and its error converges as +:math:`\mathcal{O}(h)`. The CE/CM method requires two evaluations and is thus +twice as expensive as the predictor method, but achieves an error of +:math:`\mathcal{O}(h^2)`. An exhaustive description of time integration methods +and their merits can be found in the `thesis of Colin Josey `_. OpenMC does not rely on a single time integration method but rather has several @@ -169,12 +170,14 @@ Data Considerations In principle, solving Eq. :eq:`depletion-matrix` using CRAM is fairly simple: just construct the burnup matrix at various times and solve a set of sparse -linear systems. However, constructing the burnup matrix itself involves not only -solving the transport equation to estimate transmutation reaction rates but also -a series of choices about what data to include. In OpenMC, the burnup matrix is -constructed based on data inside of a *depletion chain* file, which includes -fundamental data gathered from ENDF incident neutron, decay, and fission product -yield sublibraries. For each nuclide, this file includes: +linear systems. However, constructing the burnup matrix itself involves not +only solving the transport equation to estimate transmutation reaction rates +(in the case of transport-coupled depletion) or to obtain microscopic cross +sections (in the case of transport-independent depletion), but also a series of +choices about what data to include. In OpenMC, the burnup matrix is constructed +based on data inside of a *depletion chain* file, which includes fundamental +data gathered from ENDF incident neutron, decay, and fission product yield +sublibraries. For each nuclide, this file includes: - What transmutation reactions are possible, their Q values, and their products; - If a nuclide is not stable, what decay modes are possible, their branching @@ -185,9 +188,12 @@ yield sublibraries. For each nuclide, this file includes: Transmutation Reactions ----------------------- -OpenMC will setup tallies in a problem based on what transmutation reactions are -available in a depletion chain file, so any arbitrary number of transmutation -reactions can be tracked. The pregenerated chain files that are available on +In transport-coupled depletion, OpenMC will setup tallies in a problem based on +what transmutation reactions are available in a depletion chain file, so any +arbitrary number of transmutation reactions can be tracked. In +transport-independent depletion, OpenMC will calculate reaction rates for every +reaction that is present in both the available cross sections and the depletion +chain file. The pregenerated chain files that are available on https://openmc.org include the following transmutation reactions: fission, (n,\ :math:`\gamma`\ ), (n,2n), (n,3n), (n,4n), (n,p), and (n,\ :math:`\alpha`\ ). @@ -202,11 +208,12 @@ accurately model the branching of the capture reaction in Am241. This is complicated by the fact that the branching ratio may depend on the incident neutron energy causing capture. -OpenMC does not currently allow energy-dependent capture branching ratios. -However, the depletion chain file does allow a transmutation reaction to be -listed multiple times with different branching ratios resulting in different -products. Spectrum-averaged capture branching ratios have been computed in LWR -and SFR spectra and are available at https://openmc.org/depletion-chains. +OpenMC's transport solver does not currently allow energy-dependent capture +branching ratios. However, the depletion chain file does allow a transmutation +reaction to be listed multiple times with different branching ratios resulting +in different products. Spectrum-averaged capture branching ratios have been +computed in LWR and SFR spectra and are available at +https://openmc.org/depletion-chains. Fission Product Yields ---------------------- @@ -217,26 +224,31 @@ energies. It is an open question as to what the best way to handle this energy dependence is. OpenMC includes three methods for treating the energy dependence of FPY: -1. Use FPY data corresponding to a specified energy. +1. Use FPY data corresponding to a specified energy. This is used by default in + both transport-coupled and transport-independent depletion. 2. Tally fission rates above and below a specified cutoff energy. Assume that all fissions below the cutoff energy correspond to thermal FPY data and all - fission above the cutoff energy correspond to fast FPY data. + fission above the cutoff energy correspond to fast FPY data. Only applicable + to transport-coupled depletion. 3. Compute the average energy at which fission events occur and use an effective FPY by linearly interpolating between FPY provided at neighboring energies. + Only applicable to transport-coupled depletion. -The method can be selected through the ``fission_yield_mode`` argument to the -:class:`openmc.deplete.Operator` constructor. +The method for transport-coupled depletion can be selected through the +``fission_yield_mode`` argument to the :class:`openmc.deplete.CoupledOperator` +constructor. Power Normalization ------------------- -The reaction rates provided OpenMC are given in units of reactions per source -particle. For depletion, it is necessary to compute an absolute reaction rate in -reactions per second. To do so, the reaction rates are normalized based on a -specified power. A complete description of how this normalization can be -performed is described in :ref:`usersguide_tally_normalization`. Here, we simply -note that the main depletion class, :class:`openmc.deplete.Operator`, allows the -user to choose one of two methods for estimating the heating rate, including: +In transport-coupled depletion, the reaction rates provided OpenMC are given in +units of reactions per source particle. For depletion, it is necessary to +compute an absolute reaction rate in reactions per second. To do so, the +reaction rates are normalized based on a specified power. A complete +description of how this normalization can be performed is described in +:ref:`usersguide_tally_normalization`. Here, we simply note that the main +depletion class, :class:`openmc.deplete.CoupledOperator`, allows the user to +choose one of two methods for estimating the heating rate, including: 1. Using fixed Q values from a depletion chain file (useful for comparisons to other codes that use fixed Q values), or @@ -244,4 +256,4 @@ user to choose one of two methods for estimating the heating rate, including: energy-dependent estimate of the true heating rate. The method for normalization can be chosen through the ``normalization_mode`` -argument to the :class:`openmc.deplete.Operator` class. +argument to the :class:`openmc.deplete.CoupledOperator` class. diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index c74769be6..c7c65e14a 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -61,11 +61,13 @@ Core Functions atomic_mass atomic_weight + combine_distributions decay_constant dose_coefficients gnd_name half_life isotopes + kalbach_slope linearize thin water_density diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 9f7d8c447..0b7faceb7 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -15,16 +15,18 @@ are: 1) A transport operator 2) A time-integration scheme -The former is responsible for executing a transport code, like OpenMC, -and retaining important information required for depletion. The most common examples -are reaction rates and power normalization data. The latter is responsible for -projecting reaction rates and compositions forward in calendar time across -some step size :math:`\Delta t`, and obtaining new compositions given a power -or power density. The :class:`Operator` is provided to handle communicating with -OpenMC. Several classes are provided that implement different time-integration -algorithms for depletion calculations, which are described in detail in Colin -Josey's thesis, `Development and analysis of high order neutron -transport-depletion coupling algorithms `_. +The former is responsible for calculating and retaining important information +required for depletion. The most common examples are reaction rates and power +normalization data. The latter is responsible for projecting reaction rates and +compositions forward in calendar time across some step size :math:`\Delta t`, +and obtaining new compositions given a power or power density. The +:class:`CoupledOperator` class is provided to obtain reaction rates via tallies +through OpenMC's transport solver, and the :class:`IndependentOperator` class is +provided to obtain reaction rates from cross-section data. Several classes are +provided that implement different time-integration algorithms for depletion +calculations, which are described in detail in Colin Josey's thesis, +`Development and analysis of high order neutron transport-depletion coupling +algorithms `_. .. autosummary:: :toctree: generated @@ -40,18 +42,20 @@ transport-depletion coupling algorithms `_. SICELIIntegrator SILEQIIntegrator -Each of these classes expects a "transport operator" to be passed. An operator -specific to OpenMC is available using the following class: +Each of these classes expects a "transport operator" to be passed. OpenMC +provides The following classes implementing transpor operators: .. autosummary:: :toctree: generated :nosignatures: :template: mycallable.rst - Operator + CoupledOperator + IndependentOperator -The :class:`Operator` must also have some knowledge of how nuclides transmute -and decay. This is handled by the :class:`Chain`. +The :class:`CoupledOperator` and :class:`IndependentOperator` classes must also +have some knowledge of how nuclides transmute and decay. This is handled by the +:class:`Chain`. Minimal Example --------------- @@ -64,11 +68,12 @@ A minimal example for performing depletion would be: >>> import openmc.deplete >>> geometry = openmc.Geometry.from_xml() >>> settings = openmc.Settings.from_xml() + >>> model = openmc.model.Model(geometry, settings) # Representation of a depletion chain >>> chain_file = "chain_casl.xml" - >>> operator = openmc.deplete.Operator( - ... geometry, settings, chain_file) + >>> operator = openmc.deplete.CoupledOperator( + ... model, chain_file) # Set up 5 time steps of one day each >>> dt = [24 * 60 * 60] * 5 @@ -131,6 +136,7 @@ data, such as number densities and reaction rates for each material. :template: myclass.rst AtomNumber + MicroXS OperatorResult ReactionRates Results @@ -172,7 +178,7 @@ with :func:`cram.CRAM48` being the default. :class:`multiprocessing.pool.Pool` class. If set to ``None`` (default), the number returned by :func:`os.cpu_count` is used. -The following classes are used to help the :class:`openmc.deplete.Operator` +The following classes are used to help the :class:`openmc.deplete.CoupledOperator` compute quantities like effective fission yields, reaction rates, and total system energy. @@ -189,14 +195,45 @@ total system energy. helpers.FissionYieldCutoffHelper helpers.FluxCollapseHelper +The :class:`openmc.deplete.IndependentOperator` uses inner classes subclassed +from those listed above to perform similar calculations. + +Intermediate Classes +-------------------- + +Specific implementations of abstract base classes may utilize some of +the same methods and data structures. These methods and data are stored +in intermediate classes. + +Methods common to tally-based implementation of :class:`FissionYieldHelper` +are stored in :class:`helpers.TalliedFissionYieldHelper` + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + helpers.TalliedFissionYieldHelper + +Methods common to OpenMC-specific implementations of :class:`TransportOperator` +are stored in :class:`openmc_operator.OpenMCOperator` + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: mycallable.rst + + openmc_operator.OpenMCOperator + + Abstract Base Classes --------------------- A good starting point for extending capabilities in :mod:`openmc.deplete` is to examine the following abstract base classes. Custom classes can inherit from :class:`abc.TransportOperator` to implement alternative -schemes for collecting reaction rates and other data from a transport code -prior to depleting materials +schemes for collecting reaction rates and other data prior to depleting +materials .. autosummary:: :toctree: generated @@ -206,7 +243,9 @@ prior to depleting materials abc.TransportOperator The following classes are abstract classes used to pass information from -OpenMC simulations back on to the :class:`abc.TransportOperator` +transport simulations (in the case of transport-coupled depletion) or to +simply calculate these quantities directly (in the case of +transport-independent depletion) back on to the :class:`abc.TransportOperator` .. autosummary:: :toctree: generated @@ -216,7 +255,6 @@ OpenMC simulations back on to the :class:`abc.TransportOperator` abc.NormalizationHelper abc.FissionYieldHelper abc.ReactionRateHelper - abc.TalliedFissionYieldHelper Custom integrators or depletion solvers can be developed by subclassing from the following abstract base classes: diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index 57619611f..29c8114d9 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -4,50 +4,30 @@ Depletion and Transmutation =========================== -OpenMC supports coupled depletion, or burnup, calculations through the -:mod:`openmc.deplete` Python module. OpenMC solves the transport equation to -obtain transmutation reaction rates, and then the reaction rates are used to -solve a set of transmutation equations that determine the evolution of nuclide -densities within a material. The nuclide densities predicted as some future time -are then used to determine updated reaction rates, and the process is repeated -for as many timesteps as are requested. +OpenMC supports transport-coupled and transport-independent depletion, or +burnup, calculations through the :mod:`openmc.deplete` Python module. OpenMC +uses transmutation reaction rates to solve a set of transmutation equations +that determine the evolution of nuclide densities within a material. The +nuclide densities predicted at some future time are then used to determine +updated reaction rates, and the process is repeated for as many timesteps as +are requested. -The depletion module is designed such that the flux/reaction rate solution (the +The depletion module is designed such that the reaction rate solution (the transport "operator") is completely isolated from the solution of the -transmutation equations and the method used for advancing time. At present, the -:mod:`openmc.deplete` module offers a single transport operator, -:class:`openmc.deplete.Operator` (which uses the OpenMC transport solver), but -in principle additional operator classes based on other transport codes could be -implemented and no changes to the depletion solver itself would be needed. The -operator class requires a :class:`openmc.Geometry` instance and a -:class:`openmc.Settings` instance:: - - geom = openmc.Geometry() - settings = openmc.Settings() - ... - - op = openmc.deplete.Operator(geom, settings) - -Any material that contains a fissionable nuclide is depleted by default, but -this can behavior can be changed with the :attr:`Material.depletable` attribute. - -.. important:: The volume must be specified for each material that is depleted by - setting the :attr:`Material.volume` attribute. This is necessary - in order to calculate the proper normalization of tally results - based on the source rate. +transmutation equations and the method used for advancing time. :mod:`openmc.deplete` supports multiple time-integration methods for determining material compositions over time. Each method appears as a different class. For example, :class:`openmc.deplete.CECMIntegrator` runs a depletion calculation using the CE/CM algorithm (deplete over a timestep using the middle-of-step -reaction rates). An instance of :class:`openmc.deplete.Operator` is passed to -one of these functions along with the timesteps and power level:: +reaction rates). An instance of :class:`~openmc.deplete.abc.TransportOperator` +is passed to one of these functions along with the timesteps and power level:: power = 1200.0e6 # watts timesteps = [10.0, 10.0, 10.0] # days openmc.deplete.CECMIntegrator(op, timesteps, power, timestep_units='d').integrate() -The coupled transport-depletion problem is executed, and once it is done a +The depletion problem is executed, and once it is done a ``depletion_results.h5`` file is written. The results can be analyzed using the :class:`openmc.deplete.Results` class. This class has methods that allow for easy retrieval of k-effective, nuclide concentrations, and reaction rates over @@ -56,11 +36,41 @@ time:: results = openmc.deplete.Results("depletion_results.h5") time, keff = results.get_keff() -Note that the coupling between the transport solver and the transmutation solver -happens in-memory rather than by reading/writing files on disk. +Note that the coupling between the reaction rate solver and the transmutation +solver happens in-memory rather than by reading/writing files on disk. OpenMC +has two categories of transport operators for obtaining transmutation reaction +rates. + +.. _coupled-depletion: + +Transport-coupled depletion +=========================== + +This category of operator solves the transport equation to obtain transmutation +reaction rates. At present, the :mod:`openmc.deplete` module offers a single +transport-coupled operator, :class:`openmc.deplete.CoupledOperator` (which uses +the OpenMC transport solver), but in principle additional transport-coupled +operator classes based on other transport codes could be implemented and no +changes to the depletion solver itself would be needed. The +:class:`openmc.deplete.CoupledOperator` class requires a :class:`~openmc.Model` +instance containing material, geometry, and settings information:: + + model = openmc.Model() + ... + + op = openmc.deplete.CoupledOperator(model) + +Any material that contains a fissionable nuclide is depleted by default, but +this can behavior can be changed with the :attr:`Material.depletable` attribute. + +.. important:: + + The volume must be specified for each material that is depleted by setting + the :attr:`Material.volume` attribute. This is necessary in order to + calculate the proper normalization of tally results based on the source rate. Fixed-Source Transmutation -========================== +-------------------------- When the ``power`` or ``power_density`` argument is used for one of the Integrator classes, it is assumed that OpenMC is running in k-eigenvalue mode, @@ -77,11 +87,11 @@ using the :attr:`Material.depletable` attribute:: mat = openmc.Material() mat.depletable = True -When constructing the :class:`~openmc.deplete.Operator`, you should indicate -that normalization of tally results will be done based on the source rate rather -than a power or power density:: +When constructing the :class:`~openmc.deplete.CoupledOperator`, you should +indicate that normalization of tally results will be done based on the source +rate rather than a power or power density:: - op = openmc.deplete.Operator(geometry, settings, normalization_mode='source-rate') + op = openmc.deplete.CoupledOperator(model, normalization_mode='source-rate') Finally, when creating a depletion integrator, use the ``source_rates`` argument:: @@ -92,19 +102,22 @@ timestep in the calculation. A zero source rate for a given timestep will result in a decay-only step, where all reaction rates are zero. Caveats -======= +------- + +.. _energy-deposition: Energy Deposition ------------------ +~~~~~~~~~~~~~~~~~ The default energy deposition mode, ``"fission-q"``, instructs the -:class:`openmc.deplete.Operator` to normalize reaction rates using the product -of fission reaction rates and fission Q values taken from the depletion chain. -This approach does not consider indirect contributions to energy deposition, -such as neutron heating and energy from secondary photons. In doing this, the -energy deposited during a transport calculation will be lower than expected. -This causes the reaction rates to be over-adjusted to hit the user-specific -power, or power density, leading to an over-depletion of burnable materials. +:class:`~openmc.deplete.CoupledOperator` to normalize reaction rates using the +product of fission reaction rates and fission Q values taken from the depletion +chain. This approach does not consider indirect contributions to energy +deposition, such as neutron heating and energy from secondary photons. In doing +this, the energy deposited during a transport calculation will be lower than +expected. This causes the reaction rates to be over-adjusted to hit the +user-specific power, or power density, leading to an over-depletion of burnable +materials. There are some remedies. First, the fission Q values can be directly set in a variety of ways. This requires knowing what the total fission energy release @@ -113,29 +126,32 @@ should be, including indirect components. Some examples are provided below:: # use a dictionary of fission_q values fission_q = {"U235": 202e+6} # energy in eV + # create a Model object + model = openmc.Model(geometry, settings) + # create a modified chain and write it to a new file chain = openmc.deplete.Chain.from_xml("chain.xml", fission_q) chain.export_to_xml("chain_mod_q.xml") - op = openmc.deplete.Operator(geometry, setting, "chain_mod_q.xml") + op = openmc.deplete.CoupledOperator(model, "chain_mod_q.xml") # alternatively, pass the modified fission Q directly to the operator - op = openmc.deplete.Operator(geometry, setting, "chain.xml", + op = openmc.deplete.CoupledOperator(model, "chain.xml", fission_q=fission_q) A more complete way to model the energy deposition is to use the modified -heating reactions described in :ref:`methods_heating`. These values can be used +heating reactions described in :ref:`methods_heating`. These values can be used to normalize reaction rates instead of using the fission reaction rates with:: - op = openmc.deplete.Operator(geometry, settings, "chain.xml", + op = openmc.deplete.CoupledOperator(model, "chain.xml", normalization_mode="energy-deposition") These modified heating libraries can be generated by running the latest version -of :meth:`openmc.data.IncidentNeutron.from_njoy`, and will eventually be bundled +of :meth:`openmc.data.IncidentNeutron.from_njoy()`, and will eventually be bundled into the distributed libraries. Local Spectra and Repeated Materials ------------------------------------- +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ It is not uncommon to explicitly create a single burnable material across many locations. From a pure transport perspective, there is nothing wrong with @@ -160,7 +176,7 @@ the next transport step. This can be countered by instructing the operator to treat repeated instances of the same material as a unique material definition with:: - op = openmc.deplete.Operator(geometry, settings, chain_file, + op = openmc.deplete.CoupledOperator(model, chain_file, diff_burnable_mats=True) For our example problem, this would deplete fuel on the outer region of the @@ -177,3 +193,164 @@ across all material instances. This will increase the total memory usage and run time due to an increased number of tallies and material definitions. +Transport-independent depletion +=============================== + +.. warning:: + + This feature is still under heavy development and has yet to be rigorously + verified. API changes and feature additions are possible and likely in + the near future. + +This category of operator uses pre-calculated one-group microscopic cross +sections to obtain transmutation reaction rates. OpenMC provides the +:class:`~openmc.deplete.IndependentOperator` for this method of calculation. +While the one-group microscopic cross sections can be calculated using a +transport solver, :class:`~openmc.deplete.IndependentOperator` is not directly +coupled to any transport solver. The +:class:`~openmc.deplete.IndependentOperator` class requires a +:class:`openmc.Materials` object, a :class:`~openmc.deplete.MicroXS` object, +and a path to a depletion chain file:: + + # load in the microscopic cross sections + materials = openmc.Materials() + ... + + micro_xs = openmc.deplete.MicroXS.from_csv(micro_xs_path) + op = openmc.deplete.IndependentOperator(materials, micro_xs, chain_file) + +.. note:: + + The same statements from :ref:`coupled-depletion` about which + materials are depleted and the requirement for depletable materials to have + a specified volume also apply here. + +An alternate constructor, +:meth:`~openmc.deplete.IndependentOperator.from_nuclides`, accepts a volume and +dictionary of nuclide concentrations in place of the :class:`openmc.Materials` +object:: + + nuclides = {'U234': 8.92e18, + 'U235': 9.98e20, + 'U238': 2.22e22, + 'U236': 4.57e18, + 'O16': 4.64e22, + 'O17': 1.76e19} + volume = 0.5 + op = openmc.deplete.IndependentOperator.from_nuclides(volume, + nuclides, + micro_xs, + chain_file, + nuc_units='atom/cm3') + +A user can then define an integrator class as they would for a coupled +transport-depletion calculation and follow the same steps from there. + +.. note:: + + Ideally, one-group cross section data should be available for every + reaction in the depletion chain. If a nuclide that has a reaction + associated with it in the depletion chain is present in the `nuclides` + parameter but not the cross section data, that reaction will not be + simulated. + +Generating Microscopic Cross Sections +------------------------------------- + +Users can generate the one-group microscopic cross sections needed by +:class:`~openmc.deplete.IndependentOperator` using the +:class:`~openmc.deplete.MicroXS` class:: + + import openmc + + model = openmc.Model.from_xml() + + micro_xs = openmc.deplete.MicroXS.from_model(model, + model.materials[0], + chain_file) + +The :meth:`~openmc.deplete.MicroXS.from_model()` method will produce a +:class:`~openmc.deplete.MicroXS` object with microscopic cross section data in +units of barns, which is what :class:`~openmc.deplete.IndependentOperator` +expects the units to be. The :class:`~openmc.deplete.MicroXS` class also +includes functions to read in cross section data directly from a ``.csv`` file +or from data arrays:: + + micro_xs = MicroXS.from_csv(micro_xs_path) + + nuclides = ['U234', 'U235', 'U238'] + reactions = ['fission', '(n,gamma)'] + data = np.array([[0.1, 0.2], + [0.3, 0.4], + [0.01, 0.5]]) + micro_xs = MicroXS.from_array(nuclides, reactions, data) + +.. important:: + + Both :meth:`~openmc.deplete.MicroXS.from_csv()` and + :meth:`~openmc.deplete.MicroXS.from_array()` assume the cross section values + provided are in barns by defualt, but have no way of verifying this. Make + sure your cross sections are in the correct units before passing to a + :class:`~openmc.deplete.IndependentOperator` object. + +Caveats +------- + +Reaction Rate Normalization +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The :class:`~openmc.deplete.IndependentOperator` class supports two methods for +normalizing reaction rates: + +.. important:: + + Make sure you set the correct parameter in the :class:`openmc.abc.Integrator` + class. Use the ``source_rates`` parameter when + ``normalization_mode == source-rate``, and use ``power`` or ``power_density`` + when ``normalization_mode == fission-q``. + +1. ``source-rate`` normalization, which assumes the ``source_rate`` provided by + the time integrator is a flux, and obtains the reaction rates by multiplying + the cross-sections by the ``source-rate``. +2. ``fission-q`` normalization, which uses the ``power`` or ``power_density`` + provided by the time integrator to obtain reaction rates by computing a value + for the flux based on this power. The general equation for the flux is + + .. math:: + + \phi = \frac{P}{V \cdot \sum_i (Q_i \cdot \Sigma^f_i \cdot \rho_i)} + + where :math:`\sum_i` is the sum over all nuclides :math:`i`. This equation + makes the same assumptions and issues as discussed in + :ref:`energy-deposition`. Unfortunately, the proposed solution in that + section does not apply here since we are decoupled from transport code. + However, there is a method to converge to a more accurate value for flux by + using substeps during time integration. + `This paper `_ provides a + good discussion of this method. + +.. warning:: + + The accuracy of results when using ``fission-q`` is entirely dependent on + your depletion chain. Make sure it has sufficient data to resolve the + dynamics of your particular scenario. + +Multiple Materials +~~~~~~~~~~~~~~~~~~ + +Running a depletion simulation with multiple materials using the +``source-rate`` normalization method treats each material as completely +separate with respect to reaction rates. This can be useful for running many +different cases of a particular scenario. However, running a depletion +simulation with multiple materials using the ``fission-q`` normalization method +treats each material as part of the same "reactor" due to how ``fission-q`` +normalization accumulates energy values from each material to a single value. +This behavior may change in the future. + +Time integration +~~~~~~~~~~~~~~~~ + +The one-group microscopic cross sections passed to +:class:`openmc.deplete.IndependentOperator` are fixed values for the entire +depletion simulation. This implicit assumption may produce inaccurate results +for certain scenarios. diff --git a/openmc/data/decay.py b/openmc/data/decay.py index 51a594bc8..2c2505bf5 100644 --- a/openmc/data/decay.py +++ b/openmc/data/decay.py @@ -9,7 +9,9 @@ from uncertainties import ufloat, UFloat import openmc.checkvalue as cv from openmc.mixin import EqualityMixin +from openmc.stats import Discrete, Tabular, combine_distributions from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER +from .function import INTERPOLATION_SCHEME from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record @@ -314,6 +316,12 @@ class Decay(EqualityMixin): 'excited_state', 'mass', 'stable', 'spin', and 'parity'. spectra : dict Resulting radiation spectra for each radiation type. + sources : dict + Radioactive decay source distributions represented as a dictionary + mapping particle types (e.g., 'photon') to instances of + :class:`openmc.stats.Univariate`. + + .. versionadded:: 0.13.1 """ def __init__(self, ev_or_filename): @@ -329,6 +337,7 @@ class Decay(EqualityMixin): self.modes = [] self.spectra = {} self.average_energies = {} + self._sources = None # Get head record items = get_head_record(file_obj) @@ -495,3 +504,69 @@ class Decay(EqualityMixin): """ return cls(ev_or_filename) + + @property + def sources(self): + """Radioactive decay source distributions""" + # If property has been computed already, return it + # TODO: Replace with functools.cached_property when support is Python 3.9+ + if self._sources is not None: + return self._sources + + sources = {} + name = self.nuclide['name'] + decay_constant = self.decay_constant.n + for particle, spectra in self.spectra.items(): + # Set particle type based on 'particle' above + particle_type = { + 'gamma': 'photon', + 'beta-': 'electron', + 'ec/beta+': 'positron', + 'alpha': 'alpha', + 'n': 'neutron', + 'sf': 'fragment', + 'p': 'proton', + 'e-': 'electron', + 'xray': 'photon', + 'anti-neutrino': 'anti-neutrino', + 'neutrino': 'neutrino', + }[particle] + + if particle_type not in sources: + sources[particle_type] = [] + + # Create distribution for discrete + if spectra['continuous_flag'] in ('discrete', 'both'): + energies = [] + intensities = [] + for discrete_data in spectra['discrete']: + energies.append(discrete_data['energy'].n) + intensities.append(discrete_data['intensity'].n) + energies = np.array(energies) + intensity = spectra['discrete_normalization'].n + rates = decay_constant * intensity * np.array(intensities) + dist_discrete = Discrete(energies, rates) + sources[particle_type].append(dist_discrete) + + # Create distribution for continuous + if spectra['continuous_flag'] in ('continuous', 'both'): + f = spectra['continuous']['probability'] + if len(f.interpolation) > 1: + raise NotImplementedError("Multiple interpolation regions: {name}, {particle}") + interpolation = INTERPOLATION_SCHEME[f.interpolation[0]] + if interpolation not in ('histogram', 'linear-linear'): + raise NotImplementedError("Continuous spectra with {interpolation} interpolation ({name}, {particle}) not supported") + + intensity = spectra['continuous_normalization'].n + rates = decay_constant * intensity * f.y + dist_continuous = Tabular(f.x, rates, interpolation) + sources[particle_type].append(dist_continuous) + + # Combine discrete distributions + merged_sources = {} + for particle_type, dist_list in sources.items(): + merged_sources[particle_type] = combine_distributions( + dist_list, [1.0]*len(dist_list)) + + self._sources = merged_sources + return self._sources diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py index f4c914d90..a036da694 100644 --- a/openmc/data/kalbach_mann.py +++ b/openmc/data/kalbach_mann.py @@ -5,6 +5,7 @@ from warnings import warn import numpy as np import openmc.checkvalue as cv +from openmc.mixin import EqualityMixin from openmc.stats import Tabular, Univariate, Discrete, Mixture from .function import Tabulated1D, INTERPOLATION_SCHEME from .angle_energy import AngleEnergy @@ -12,6 +13,240 @@ from .data import EV_PER_MEV from .endf import get_list_record, get_tab2_record +class _AtomicRepresentation(EqualityMixin): + """Atomic representation of an isotope or a particle. + + Parameters + ---------- + z : int + Number of protons (atomic number) + a : int + Number of nucleons (mass number) + + Raises + ------ + ValueError + When the number of protons (z) declared is higher than the number + of nucleons (a) + + Attributes + ---------- + z : int + Number of protons (atomic number) + a : int + Number of nucleons (mass number) + n : int + Number of neutrons + za : int + ZA identifier, 1000*Z + A, where Z is the atomic number and A the mass + number + + """ + def __init__(self, z, a): + # Sanity checks on values + cv.check_type('z', z, Integral) + cv.check_greater_than('z', z, 0, equality=True) + cv.check_type('a', a, Integral) + cv.check_greater_than('a', a, 0, equality=True) + if z > a: + raise ValueError(f"Number of protons ({z}) must be less than or " + f"equal to number of nucleons ({a}).") + + self._z = z + self._a = a + + def __add__(self, other): + """Add two _AtomicRepresentations""" + z = self.z + other.z + a = self.a + other.a + return _AtomicRepresentation(z=z, a=a) + + def __sub__(self, other): + """Substract two _AtomicRepresentations""" + z = self.z - other.z + a = self.a - other.a + return _AtomicRepresentation(z=z, a=a) + + @property + def a(self): + return self._a + + @property + def z(self): + return self._z + + @property + def n(self): + return self.a - self.z + + @property + def za(self): + return self.z * 1000 + self.a + + @classmethod + def from_za(cls, za): + """Instantiate an _AtomicRepresentation from a ZA identifier. + + Parameters + ---------- + za : int + ZA identifier, 1000*Z + A, where Z is the atomic number and A the + mass number + + Returns + ------- + _AtomicRepresentation + Atomic representation of the isotope/particle + + """ + z, a = divmod(za, 1000) + return cls(z, a) + + +def _separation_energy(compound, nucleus, particle): + """Calculates the separation energy as defined in ENDF-6 manual + BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1 + and LANG=2. This function can be used for the incident or emitted + particle of the following reaction: A + a -> C -> B + b + + Parameters + ---------- + compound : _AtomicRepresentation + Atomic representation of the compound (C) + nucleus : _AtomicRepresentation + Atomic representation of the nucleus (A or B) + particle : _AtomicRepresentation + Atomic representation of the particle (a or b) + + Returns + ------- + separation_energy : float + Separation energy in MeV + + """ + # Determine A, Z, and N for compound and nucleus + A_c = compound.a + Z_c = compound.z + N_c = compound.n + A_a = nucleus.a + Z_a = nucleus.z + N_a = nucleus.n + + # Determine breakup energy of incident particle (ENDF-6 Formats Manual, + # Appendix H, Table 3) in MeV + za_to_breaking_energy = { + 1: 0.0, + 1001: 0.0, + 1002: 2.224566, + 1003: 8.481798, + 2003: 7.718043, + 2004: 28.29566 + } + I_a = za_to_breaking_energy[particle.za] + + # Eq. 4 in in doi:10.1103/PhysRevC.37.2350 or ENDF-6 Formats Manual section + # 6.2.3.2 + return ( + 15.68 * (A_c - A_a) - + 28.07 * ((N_c - Z_c)**2 / A_c - (N_a - Z_a)**2 / A_a) - + 18.56 * (A_c**(2./3.) - A_a**(2./3.)) + + 33.22 * ((N_c - Z_c)**2 / A_c**(4./3.) - (N_a - Z_a)**2 / A_a**(4./3.)) - + 0.717 * (Z_c**2 / A_c**(1./3.) - Z_a**2 / A_a**(1./3.)) + + 1.211 * (Z_c**2 / A_c - Z_a**2 / A_a) - + I_a + ) + + +def kalbach_slope(energy_projectile, energy_emitted, za_projectile, + za_emitted, za_target): + """Returns Kalbach-Mann slope from calculations. + + The associated reaction is defined as: + A + a -> C -> B + b + + Where: + + - A is the targeted nucleus, + - a is the projectile, + - C is the compound, + - B is the residual nucleus, + - b is the emitted particle. + + The Kalbach-Mann slope calculation is done as defined in ENDF-6 manual + BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1 and + LANG=2. One exception to this, is that the entrance and emission channel + energies are not calculated with the AWR number, but approximated with + the number of mass instead. + + Parameters + ---------- + energy_projectile : float + Energy of the projectile in the laboratory system in eV + energy_emitted : float + Energy of the emitted particle in the center of mass system in eV + za_projectile : int + ZA identifier of the projectile + za_emitted : int + ZA identifier of the emitted particle + za_target : int + ZA identifier of the targeted nucleus + + Raises + ------ + NotImplementedError + When the projectile is not a neutron + + Returns + ------- + slope : float + Kalbach-Mann slope given with the same format as ACE file. + + """ + # TODO: develop for photons as projectile + # TODO: test for other particles than neutron + if za_projectile != 1: + raise NotImplementedError( + "Developed and tested for neutron projectile only." + ) + + # Special handling of elemental carbon + if za_emitted == 6000: + za_emitted = 6012 + if za_target == 6000: + za_target = 6012 + + projectile = _AtomicRepresentation.from_za(za_projectile) + emitted = _AtomicRepresentation.from_za(za_emitted) + target = _AtomicRepresentation.from_za(za_target) + compound = projectile + target + residual = compound - emitted + + # Calculate entrance and emission channel energy in MeV, defined in section + # 6.2.3.2 in the ENDF-6 Formats Manual + epsilon_a = energy_projectile * target.a / (target.a + projectile.a) / EV_PER_MEV + epsilon_b = energy_emitted * (residual.a + emitted.a) \ + / (residual.a * EV_PER_MEV) + + # Calculate separation energies using Eq. 4 in doi:10.1103/PhysRevC.37.2350 + # or ENDF-6 Formats Manual section 6.2.3.2 + s_a = _separation_energy(compound, target, projectile) + s_b = _separation_energy(compound, residual, emitted) + + # See Eq. 10 in doi:10.1103/PhysRevC.37.2350 or section 6.2.3.2 in the + # ENDF-6 Formats Manual + za_to_M = {1: 1.0, 1001: 1.0, 1002: 1.0, 2004: 0.0} + za_to_m = {1: 0.5, 1001: 1.0, 1002: 1.0, 1003: 1.0, 2003: 1.0, 2004: 2.0} + M = za_to_M[projectile.za] + m = za_to_m[emitted.za] + e_a = epsilon_a + s_a + e_b = epsilon_b + s_b + r_1 = min(e_a, 130.) + r_3 = min(e_a, 41.) + x_1 = r_1 * e_b / e_a + x_3 = r_3 * e_b / e_a + return 0.04 * x_1 + 1.8e-6 * x_1**3 + 6.7e-7 * M * m * x_3**4 + + class KalbachMann(AngleEnergy): """Kalbach-Mann distribution @@ -319,7 +554,7 @@ class KalbachMann(AngleEnergy): n_energy_out = int(ace.xss[idx + 1]) data = ace.xss[idx + 2:idx + 2 + 5*n_energy_out].copy() data.shape = (5, n_energy_out) - data[0,:] *= EV_PER_MEV + data[0, :] *= EV_PER_MEV # Create continuous distribution eout_continuous = Tabular(data[0][n_discrete_lines:], @@ -352,13 +587,28 @@ class KalbachMann(AngleEnergy): return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a) @classmethod - def from_endf(cls, file_obj): - """Generate Kalbach-Mann distribution from an ENDF evaluation + def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass): + """Generate Kalbach-Mann distribution from an ENDF evaluation. + + If the projectile is a neutron, the slope is calculated when it is + not given explicitly. Parameters ---------- file_obj : file-like object ENDF file positioned at the start of the Kalbach-Mann distribution + za_emitted : int + ZA identifier of the emitted particle + za_target : int + ZA identifier of the target + projectile_mass : float + Mass of the projectile + + Warns + ----- + UserWarning + If the mass of the projectile is not equal to 1 (other than + a neutron), the slope is not calculated and set to 0 if missing. Returns ------- @@ -374,6 +624,7 @@ class KalbachMann(AngleEnergy): energy_out = [] precompound = [] slope = [] + calculated_slope = [] for i in range(ne): items, values = get_list_record(file_obj) energy[i] = items[1] @@ -385,19 +636,46 @@ class KalbachMann(AngleEnergy): values.shape = (n_energy_out, n_angle + 2) # Outgoing energy distribution at the i-th incoming energy - eout_i = values[:,0] - eout_p_i = values[:,1] + eout_i = values[:, 0] + eout_p_i = values[:, 1] energy_out_i = Tabular(eout_i, eout_p_i, INTERPOLATION_SCHEME[lep]) energy_out.append(energy_out_i) - # Precompound and slope factors for Kalbach-Mann - r_i = values[:,2] + # Precompound factors for Kalbach-Mann + r_i = values[:, 2] + + # Slope factors for Kalbach-Mann if n_angle == 2: - a_i = values[:,3] + a_i = values[:, 3] + calculated_slope.append(False) else: - a_i = np.zeros_like(r_i) + # Check if the projectile is not a neutron + if not np.isclose(projectile_mass, 1.0, atol=1.0e-12, rtol=0.): + warn( + "Kalbach-Mann slope calculation is only available with " + "neutrons as projectile. Slope coefficients are set to 0." + ) + a_i = np.zeros_like(r_i) + calculated_slope.append(False) + + else: + # TODO: retrieve ZA of the projectile + za_projectile = 1 + a_i = [kalbach_slope(energy_projectile=energy[i], + energy_emitted=e, + za_projectile=za_projectile, + za_emitted=za_emitted, + za_target=za_target) + for e in eout_i] + calculated_slope.append(True) + precompound.append(Tabulated1D(eout_i, r_i)) slope.append(Tabulated1D(eout_i, a_i)) - return cls(tab2.breakpoints, tab2.interpolation, energy, - energy_out, precompound, slope) + km_distribution = cls(tab2.breakpoints, tab2.interpolation, energy, + energy_out, precompound, slope) + + # List of bool to indicate slope calculation by OpenMC + km_distribution._calculated_slope = calculated_slope + + return km_distribution diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py index 305edaf4f..d8ecaae18 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -127,7 +127,7 @@ acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%% 1 0 1 .{ext} / '{library}: {zsymam} at {temperature}'/ {mat} {temperature} -1 1/ +1 1 {ismooth}/ / """ @@ -248,7 +248,8 @@ def make_pendf(filename, pendf='pendf', error=0.001, stdout=False): def make_ace(filename, temperatures=None, acer=True, xsdir=None, output_dir=None, pendf=False, error=0.001, broadr=True, - heatr=True, gaspr=True, purr=True, evaluation=None, **kwargs): + heatr=True, gaspr=True, purr=True, evaluation=None, + smoothing=True, **kwargs): """Generate incident neutron ACE file from an ENDF file File names can be passed to @@ -298,6 +299,8 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None, evaluation : openmc.data.endf.Evaluation, optional If the ENDF file contains multiple material evaluations, this argument indicates which evaluation should be used. + smoothing : bool, optional + If the smoothing option (ACER card 6) is on (True) or off (False). **kwargs Keyword arguments passed to :func:`openmc.data.njoy.run` @@ -380,6 +383,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None, # acer if acer: + ismooth = int(smoothing) nacer_in = nlast for i, temperature in enumerate(temperatures): # Extend input with an ACER run for each temperature diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 5e4287f16..a2431cf1c 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -80,6 +80,14 @@ def _get_products(ev, mt): mt : int The MT value of the reaction to get products for + Raises + ------ + IOError + When the Kalbach-Mann systematics is used, but the product + is not defined in the 'center-of-mass' system. The breakup logic + is not implemented which can lead to this error being raised while + the definition of the product is correct. + Returns ------- products : list of openmc.data.Product @@ -141,7 +149,26 @@ def _get_products(ev, mt): if lang == 1: p.distribution = [CorrelatedAngleEnergy.from_endf(file_obj)] elif lang == 2: - p.distribution = [KalbachMann.from_endf(file_obj)] + # Products need to be described in the center-of-mass system + product_center_of_mass = False + if reference_frame == 'center-of-mass': + product_center_of_mass = True + elif reference_frame == 'light-heavy': + product_center_of_mass = (awr <= 4.0) + # TODO: 'breakup' logic not implemented + + if product_center_of_mass is False: + raise IOError( + "Kalbach-Mann representation must be defined in the " + "'center-of-mass' system" + ) + + zat = ev.target["atomic_number"] * 1000 + ev.target["mass_number"] + projectile_mass = ev.projectile["mass"] + p.distribution = [KalbachMann.from_endf(file_obj, + za, + zat, + projectile_mass)] elif law == 2: # Discrete two-body scattering diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py index 5891555a2..329ee8b52 100644 --- a/openmc/deplete/__init__.py +++ b/openmc/deplete/__init__.py @@ -7,7 +7,10 @@ A depletion front-end tool. from .nuclide import * from .chain import * -from .operator import * +from .openmc_operator import * +from .coupled_operator import * +from .independent_operator import * +from .microxs import * from .reaction_rates import * from .atom_number import * from .stepresult import * diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 8bbe126df..e3587f052 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -1,7 +1,6 @@ -"""function module. +"""abc module. -This module contains the Operator class, which is then passed to an integrator -to run a full depletion simulation. +This module contains Abstract Base Classes for implementing operator, integrator, depletion system solver, and operator helper classes """ from abc import ABC, abstractmethod @@ -29,8 +28,8 @@ from .pool import deplete __all__ = [ - "OperatorResult", "TransportOperator", "ReactionRateHelper", - "NormalizationHelper", "FissionYieldHelper", + "OperatorResult", "TransportOperator", + "ReactionRateHelper", "NormalizationHelper", "FissionYieldHelper", "Integrator", "SIIntegrator", "DepSystemSolver", "add_params"] @@ -83,7 +82,8 @@ class TransportOperator(ABC): operator that takes a vector of material compositions and returns an eigenvalue and reaction rates. This abstract class sets the requirements for such a transport operator. Users should instantiate - :class:`openmc.deplete.Operator` rather than this class. + :class:`openmc.deplete.CoupledOperator` or + :class:`openmc.deplete.IndependentOperator` rather than this class. Parameters ---------- @@ -221,9 +221,11 @@ class ReactionRateHelper(ABC): Parameters ---------- n_nucs : int - Number of burnable nuclides tracked by :class:`openmc.deplete.Operator` + Number of burnable nuclides tracked by + :class:`openmc.deplete.abc.TransportOperator` n_react : int - Number of reactions tracked by :class:`openmc.deplete.Operator` + Number of reactions tracked by + :class:`openmc.deplete.abc.TransportOperator` Attributes ---------- @@ -292,9 +294,9 @@ class NormalizationHelper(ABC): """Abstract class for obtaining normalization factor on tallies This helper class determines how reaction rates calculated by an instance of - :class:`openmc.deplete.Operator` should be normalized for the purpose of - constructing a burnup matrix. Based on the method chosen, the power or - source rate provided by the user, and reaction rates from a + :class:`openmc.deplete.abc.TransportOperator` should be normalized for the + purpose of constructing a burnup matrix. Based on the method chosen, the + power or source rate provided by the user, and reaction rates from a :class:`ReactionRateHelper`, this class will scale reaction rates to the correct values. @@ -302,7 +304,7 @@ class NormalizationHelper(ABC): ---------- nuclides : list of str All nuclides with desired reaction rates. Ordered to be - consistent with :class:`openmc.deplete.Operator` + consistent with :class:`openmc.deplete.abc.TransportOperator` """ @@ -316,9 +318,9 @@ class NormalizationHelper(ABC): def prepare(self, chain_nucs, rate_index): """Perform work needed to obtain energy produced - This method is called prior to the transport simulations - in :meth:`openmc.deplete.Operator.initial_condition`. Only used for - energy-based normalization. + This method is called prior to calculating the reaction rates + in :meth:`openmc.deplete.abc.TransportOperator.initial_condition`. Only + used for energy-based normalization. Parameters ---------- @@ -329,7 +331,7 @@ class NormalizationHelper(ABC): `fission_rates` for :meth:`update`. """ - def update(self, fission_rates): + def update(self, fission_rates, mat_index=None): """Update the normalization based on fission rates (only used for energy-based normalization) @@ -339,6 +341,8 @@ class NormalizationHelper(ABC): fission reaction rate for each isotope in the specified material. Should be ordered corresponding to initial ``rate_index`` used in :meth:`prepare` + mat_index : int + Material index """ @property @@ -429,7 +433,7 @@ class FissionYieldHelper(ABC): def unpack(): """Unpack tally data prior to compute fission yields. - Called after a :meth:`openmc.deplete.Operator.__call__` + Called after a :meth:`openmc.deplete.abc.TransportOperator.__call__` routine during the normalization of reaction rates. Not necessary for all subclasses to implement, unless tallies @@ -450,7 +454,7 @@ class FissionYieldHelper(ABC): mat_indexes : iterable of int Indices of tallied materials that will have their fission yields computed by this helper. Necessary as the - :class:`openmc.deplete.Operator` that uses this helper + :class:`openmc.deplete.CoupledOperator` that uses this helper may only burn a subset of all materials when running in parallel mode. """ @@ -462,14 +466,15 @@ class FissionYieldHelper(ABC): ---------- nuclides : iterable of str Nuclides with non-zero densities from the - :class:`openmc.deplete.Operator` + :class:`openmc.deplete.abc.TransportOperator` Returns ------- nuclides : list of str - Union of nuclides that the :class:`openmc.deplete.Operator` - says have non-zero densities at this stage and those that - have yield data. Sorted by nuclide name + Union of nuclides that the + :class:`openmc.deplete.abc.TransportOperator` says have non-zero + densities at this stage and those that have yield data. Sorted by + nuclide name """ return sorted(self._chain_set & set(nuclides)) @@ -483,7 +488,7 @@ class FissionYieldHelper(ABC): Parameters ---------- - operator : openmc.deplete.TransportOperator + operator : openmc.deplete.abc.TransportOperator Operator with a depletion chain kwargs: optional Additional keyword arguments to be used in constuction @@ -504,7 +509,7 @@ class Integrator(ABC): _params = r""" Parameters ---------- - operator : openmc.deplete.TransportOperator + operator : openmc.deplete.abc.TransportOperator Operator to perform transport simulations timesteps : iterable of float or iterable of tuple Array of timesteps. Note that values are not cumulative. The units are @@ -522,9 +527,10 @@ class Integrator(ABC): power_density : float or iterable of float, optional Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` - is not speficied. + is not specified. source_rates : float or iterable of float, optional - Source rate in [neutron/sec] for each interval in :attr:`timesteps` + Source rate in [neutron/sec] or neutron flux in [neutron/s-cm^2] for + each interval in :attr:`timesteps` .. versionadded:: 0.12.1 timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} @@ -546,7 +552,7 @@ class Integrator(ABC): Attributes ---------- - operator : openmc.deplete.TransportOperator + operator : openmc.deplete.abc.TransportOperator Operator to perform transport simulations chain : openmc.deplete.Chain Depletion chain @@ -841,8 +847,8 @@ class SIIntegrator(Integrator): _params = r""" Parameters ---------- - operator : openmc.deplete.TransportOperator - The operator object to simulate on. + operator : openmc.deplete.abc.TransportOperator + Operator to perform transport simulations timesteps : iterable of float or iterable of tuple Array of timesteps. Note that values are not cumulative. The units are specified by the `timestep_units` argument when `timesteps` is an @@ -859,9 +865,10 @@ class SIIntegrator(Integrator): power_density : float or iterable of float, optional Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` - is not speficied. + is not specified. source_rates : float or iterable of float, optional - Source rate in [neutron/sec] for each interval in :attr:`timesteps` + Source rate in [neutron/sec] or neutron flux in [neutron/s-cm^2] for + each interval in :attr:`timesteps` .. versionadded:: 0.12.1 timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} @@ -886,7 +893,7 @@ class SIIntegrator(Integrator): Attributes ---------- - operator : openmc.deplete.TransportOperator + operator : openmc.deplete.abc.TransportOperator Operator to perform transport simulations chain : openmc.deplete.Chain Depletion chain @@ -1006,7 +1013,7 @@ class DepSystemSolver(ABC): Parameters ---------- A : scipy.sparse.csr_matrix - Sparse transmutation matrix ``A[j, i]`` desribing rates at + Sparse transmutation matrix ``A[j, i]`` describing rates at which isotope ``i`` transmutes to isotope ``j`` n0 : numpy.ndarray Initial compositions, typically given in number of atoms in some diff --git a/openmc/deplete/atom_number.py b/openmc/deplete/atom_number.py index 5430ddc35..78ceecca6 100644 --- a/openmc/deplete/atom_number.py +++ b/openmc/deplete/atom_number.py @@ -126,6 +126,23 @@ class AtomNumber: return [nuc for nuc, ind in self.index_nuc.items() if ind < self.n_nuc_burn] + def get_mat_volume(self, mat): + """Return material volume + + Parameters + ---------- + mat : str, int, openmc.Material, or slice + Material index. + + Returns + ------- + float + Material volume in [cm^3] + + """ + mat = self._get_mat_index(mat) + return self.volume[mat] + def get_atom_density(self, mat, nuc): """Return atom density of given material and nuclide diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 0c1c34cf7..d63bc6bfb 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -264,7 +264,8 @@ class Chain: requires a list of ENDF incident neutron, decay, and neutron fission product yield sublibrary files. The depletion chain used during a depletion simulation is indicated by either an argument to - :class:`openmc.deplete.Operator` or through the + :class:`openmc.deplete.CoupledOperator` or + :class:`openmc.deplete.IndependentOperator`, or through the ``depletion_chain`` item in the :envvar:`OPENMC_CROSS_SECTIONS` environment variable. diff --git a/openmc/deplete/coupled_operator.py b/openmc/deplete/coupled_operator.py new file mode 100644 index 000000000..eba90f60d --- /dev/null +++ b/openmc/deplete/coupled_operator.py @@ -0,0 +1,546 @@ +"""Transport-coupled transport operator for depletion. + +This module implements a transport operator coupled to OpenMC's transport solver +so that it can be used by depletion integrators. The implementation makes use of +the Python bindings to OpenMC's C API so that reading tally results and updating +material number densities is all done in-memory instead of through the +filesystem. + +""" + +import copy +import os +from warnings import warn + +import numpy as np +from uncertainties import ufloat + +import openmc +from openmc.checkvalue import check_value +from openmc.data import DataLibrary +from openmc.exceptions import DataError +import openmc.lib +from openmc.mpi import comm +from .abc import OperatorResult +from .chain import _find_chain_file +from .openmc_operator import OpenMCOperator, _distribute +from .results import Results +from .helpers import ( + DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, + FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper, + SourceRateHelper, FluxCollapseHelper) + + +__all__ = ["CoupledOperator", "Operator", "OperatorResult"] + + +def _find_cross_sections(model): + """Determine cross sections to use for depletion + + Parameters + ---------- + model : openmc.model.Model + Reactor model + + """ + if model.materials and model.materials.cross_sections is not None: + # Prefer info from Model class if available + return model.materials.cross_sections + + # otherwise fallback to environment variable + cross_sections = os.environ.get("OPENMC_CROSS_SECTIONS") + if cross_sections is None: + raise DataError( + "Cross sections were not specified in Model.materials and " + "the OPENMC_CROSS_SECTIONS environment variable is not set." + ) + return cross_sections + +def _get_nuclides_with_data(cross_sections): + """Loads cross_sections.xml file to find nuclides with neutron data + + Parameters + ---------- + cross_sections : str + Path to cross_sections.xml file + + Returns + ------- + nuclides : set of str + Set of nuclide names that have cross secton data + + """ + nuclides = set() + data_lib = DataLibrary.from_xml(cross_sections) + for library in data_lib.libraries: + if library['type'] != 'neutron': + continue + for name in library['materials']: + if name not in nuclides: + nuclides.add(name) + + return nuclides + +class CoupledOperator(OpenMCOperator): + """Transport-coupled transport operator. + + Instances of this class can be used to perform transport-coupled depletion + using OpenMC's transport solver. Normally, a user needn't call methods of + this class directly. Instead, an instance of this class is passed to an + integrator class, such as :class:`openmc.deplete.CECMIntegrator`. + + .. versionchanged:: 0.13.0 + The geometry and settings parameters have been replaced with a + model parameter that takes a :class:`~openmc.model.Model` object + + Parameters + ---------- + model : openmc.model.Model + OpenMC model object + chain_file : str, optional + Path to the depletion chain XML file. Defaults to the file + listed under ``depletion_chain`` in + :envvar:`OPENMC_CROSS_SECTIONS` environment variable. + prev_results : Results, optional + Results from a previous depletion calculation. If this argument is + specified, the depletion calculation will start from the latest state + in the previous results. + diff_burnable_mats : bool, optional + Whether to differentiate burnable materials with multiple instances. + Volumes are divided equally from the original material volume. + normalization_mode : {"energy-deposition", "fission-q", "source-rate"} + Indicate how tally results should be normalized. ``"energy-deposition"`` + computes the total energy deposited in the system and uses the ratio of + the power to the energy produced as a normalization factor. + ``"fission-q"`` uses the fission Q values from the depletion chain to + compute the total energy deposited. ``"source-rate"`` normalizes + tallies based on the source rate (for fixed source calculations). + fission_q : dict, optional + Dictionary of nuclides and their fission Q values [eV]. If not given, + values will be pulled from the ``chain_file``. Only applicable + if ``"normalization_mode" == "fission-q"`` + dilute_initial : float, optional + Initial atom density [atoms/cm^3] to add for nuclides that are zero + in initial condition to ensure they exist in the decay chain. + Only done for nuclides with reaction rates. + fission_yield_mode : {"constant", "cutoff", "average"} + Key indicating what fission product yield scheme to use. The + key determines what fission energy helper is used: + + * "constant": :class:`~openmc.deplete.helpers.ConstantFissionYieldHelper` + * "cutoff": :class:`~openmc.deplete.helpers.FissionYieldCutoffHelper` + * "average": :class:`~openmc.deplete.helpers.AveragedFissionYieldHelper` + + The documentation on these classes describe their methodology + and differences. Default: ``"constant"`` + fission_yield_opts : dict of str to option, optional + Optional arguments to pass to the helper determined by + ``fission_yield_mode``. Will be passed directly on to the + helper. Passing a value of None will use the defaults for + the associated helper. + reaction_rate_mode : {"direct", "flux"}, optional + Indicate how one-group reaction rates should be calculated. The "direct" + method tallies transmutation reaction rates directly. The "flux" method + tallies a multigroup flux spectrum and then collapses one-group reaction + rates after a transport solve (with an option to tally some reaction + rates directly). + + .. versionadded:: 0.12.1 + reaction_rate_opts : dict, optional + Keyword arguments that are passed to the reaction rate helper class. + When ``reaction_rate_mode`` is set to "flux", energy group boundaries + can be set using the "energies" key. See the + :class:`~openmc.deplete.helpers.FluxCollapseHelper` class for all + options. + + .. versionadded:: 0.12.1 + reduce_chain : bool, optional + If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the + depletion chain up to ``reduce_chain_level``. + + .. versionadded:: 0.12 + reduce_chain_level : int, optional + Depth of the search when reducing the depletion chain. Only used + if ``reduce_chain`` evaluates to true. The default value of + ``None`` implies no limit on the depth. + + .. versionadded:: 0.12 + + Attributes + ---------- + model : openmc.model.Model + OpenMC model object + geometry : openmc.Geometry + OpenMC geometry object + settings : openmc.Settings + OpenMC settings object + dilute_initial : float + Initial atom density [atoms/cm^3] to add for nuclides that + are zero in initial condition to ensure they exist in the decay + chain. Only done for nuclides with reaction rates. + output_dir : pathlib.Path + Path to output directory to save results. + round_number : bool + Whether or not to round output to OpenMC to 8 digits. + Useful in testing, as OpenMC is incredibly sensitive to exact values. + number : openmc.deplete.AtomNumber + Total number of atoms in simulation. + nuclides_with_data : set of str + A set listing all unique nuclides available from cross_sections.xml. + chain : openmc.deplete.Chain + The depletion chain information necessary to form matrices and tallies. + reaction_rates : openmc.deplete.ReactionRates + Reaction rates from the last operator step. + burnable_mats : list of str + All burnable material IDs + heavy_metal : float + Initial heavy metal inventory [g] + local_mats : list of str + All burnable material IDs being managed by a single process + prev_res : Results or None + Results from a previous depletion calculation. ``None`` if no + results are to be used. + cleanup_when_done : bool + Whether to finalize and clear the shared library memory when the + depletion operation is complete. Defaults to clearing the library. + """ + _fission_helpers = { + "average": AveragedFissionYieldHelper, + "constant": ConstantFissionYieldHelper, + "cutoff": FissionYieldCutoffHelper, + } + + def __init__(self, model, chain_file=None, prev_results=None, + diff_burnable_mats=False, normalization_mode="fission-q", + fission_q=None, dilute_initial=1.0e3, + fission_yield_mode="constant", fission_yield_opts=None, + reaction_rate_mode="direct", reaction_rate_opts=None, + reduce_chain=False, reduce_chain_level=None): + + # check for old call to constructor + if isinstance(model, openmc.Geometry): + msg = "As of version 0.13.0 openmc.deplete.CoupledOperator " \ + "requires an openmc.Model object rather than the " \ + "openmc.Geometry and openmc.Settings parameters. Please use " \ + "the geometry and settings objects passed here to create a " \ + " model with which to generate the transport Operator." + raise TypeError(msg) + + # Determine cross sections / depletion chain + cross_sections = _find_cross_sections(model) + if chain_file is None: + chain_file = _find_chain_file(cross_sections) + + check_value('fission yield mode', fission_yield_mode, + self._fission_helpers.keys()) + check_value('normalization mode', normalization_mode, + ('energy-deposition', 'fission-q', 'source-rate')) + if normalization_mode != "fission-q": + if fission_q is not None: + warn("Fission Q dictionary will not be used") + fission_q = None + self.model = model + self.settings = model.settings + self.geometry = model.geometry + + # determine set of materials in the model + if not model.materials: + model.materials = openmc.Materials( + model.geometry.get_all_materials().values() + ) + + self.cleanup_when_done = True + + if reaction_rate_opts is None: + reaction_rate_opts = {} + if fission_yield_opts is None: + fission_yield_opts = {} + helper_kwargs = { + 'reaction_rate_mode': reaction_rate_mode, + 'normalization_mode': normalization_mode, + 'fission_yield_mode': fission_yield_mode, + 'reaction_rate_opts': reaction_rate_opts, + 'fission_yield_opts': fission_yield_opts + } + + super().__init__( + model.materials, + cross_sections, + chain_file, + prev_results, + diff_burnable_mats, + fission_q, + dilute_initial, + helper_kwargs, + reduce_chain, + reduce_chain_level) + + def _differentiate_burnable_mats(self): + """Assign distribmats for each burnable material""" + + # Count the number of instances for each cell and material + self.geometry.determine_paths(instances_only=True) + + # Extract all burnable materials which have multiple instances + distribmats = set( + [mat for mat in self.materials + if mat.depletable and mat.num_instances > 1]) + + for mat in distribmats: + if mat.volume is None: + raise RuntimeError("Volume not specified for depletable " + "material with ID={}.".format(mat.id)) + mat.volume /= mat.num_instances + + if distribmats: + # Assign distribmats to cells + for cell in self.geometry.get_all_material_cells().values(): + if cell.fill in distribmats: + mat = cell.fill + cell.fill = [mat.clone() + for i in range(cell.num_instances)] + + self.materials = openmc.Materials( + self.model.geometry.get_all_materials().values() + ) + + def _load_previous_results(self): + """Load results from a previous depletion simulation""" + # Reload volumes into geometry + self.prev_res[-1].transfer_volumes(self.model) + + # Store previous results in operator + # Distribute reaction rates according to those tracked + # on this process + if comm.size != 1: + prev_results = self.prev_res + self.prev_res = Results() + mat_indexes = _distribute(range(len(self.burnable_mats))) + for res_obj in prev_results: + new_res = res_obj.distribute(self.local_mats, mat_indexes) + self.prev_res.append(new_res) + + def _get_nuclides_with_data(self, cross_sections): + """Loads cross_sections.xml file to find nuclides with neutron data + + Parameters + ---------- + cross_sections : str + Path to cross_sections.xml file + + Returns + ------- + nuclides : set of str + Set of nuclide names that have cross secton data + + """ + return _get_nuclides_with_data(cross_sections) + + def _get_helper_classes(self, helper_kwargs): + """Create the ``_rate_helper``, ``_normalization_helper``, and + ``_yield_helper`` objects. + + Parameters + ---------- + helper_kwargs : dict + Keyword arguments for helper classes + + """ + reaction_rate_mode = helper_kwargs['reaction_rate_mode'] + normalization_mode = helper_kwargs['normalization_mode'] + fission_yield_mode = helper_kwargs['fission_yield_mode'] + reaction_rate_opts = helper_kwargs['reaction_rate_opts'] + fission_yield_opts = helper_kwargs['fission_yield_opts'] + + # Get classes to assist working with tallies + if reaction_rate_mode == "direct": + self._rate_helper = DirectReactionRateHelper( + self.reaction_rates.n_nuc, self.reaction_rates.n_react) + elif reaction_rate_mode == "flux": + # Ensure energy group boundaries were specified + if 'energies' not in reaction_rate_opts: + raise ValueError( + "Energy group boundaries must be specified in the " + "reaction_rate_opts argument when reaction_rate_mode is" + "set to 'flux'.") + + self._rate_helper = FluxCollapseHelper( + self.reaction_rates.n_nuc, + self.reaction_rates.n_react, + **reaction_rate_opts + ) + else: + raise ValueError("Invalid reaction rate mode.") + + if normalization_mode == "fission-q": + self._normalization_helper = ChainFissionHelper() + elif normalization_mode == "energy-deposition": + score = "heating" if self.settings.photon_transport else "heating-local" + self._normalization_helper = EnergyScoreHelper(score) + else: + self._normalization_helper = SourceRateHelper() + + # Select and create fission yield helper + fission_helper = self._fission_helpers[fission_yield_mode] + self._yield_helper = fission_helper.from_operator( + self, **fission_yield_opts) + + def initial_condition(self): + """Performs final setup and returns initial condition. + + Returns + ------- + list of numpy.ndarray + Total density for initial conditions. + + """ + + # Create XML files + if comm.rank == 0: + self.geometry.export_to_xml() + self.settings.export_to_xml() + self._generate_materials_xml() + + # Initialize OpenMC library + comm.barrier() + if not openmc.lib.is_initialized: + openmc.lib.init(intracomm=comm) + + # Generate tallies in memory + materials = [openmc.lib.materials[int(i)] for i in self.burnable_mats] + + return super().initial_condition(materials) + + def _generate_materials_xml(self): + """Creates materials.xml from self.number. + + Due to uncertainty with how MPI interacts with OpenMC API, this + constructs the XML manually. The long term goal is to do this + through direct memory writing. + + """ + # Sort nuclides according to order in AtomNumber object + nuclides = list(self.number.nuclides) + for mat in self.materials: + mat._nuclides.sort(key=lambda x: nuclides.index(x[0])) + + self.materials.export_to_xml() + + def __call__(self, vec, source_rate): + """Runs a simulation. + + Simulation will abort under the following circumstances: + + 1) No energy is computed using OpenMC tallies. + + Parameters + ---------- + vec : list of numpy.ndarray + Total atoms to be used in function. + source_rate : float + Power in [W] or source rate in [neutron/sec] + + Returns + ------- + openmc.deplete.OperatorResult + Eigenvalue and reaction rates resulting from transport operator + + """ + # Reset results in OpenMC + openmc.lib.reset() + + self._update_materials_and_nuclides(vec) + + # If the source rate is zero, return zero reaction rates without running + # a transport solve + if source_rate == 0.0: + rates = self.reaction_rates.copy() + rates.fill(0.0) + return OperatorResult(ufloat(0.0, 0.0), rates) + + # Run OpenMC + openmc.lib.run() + openmc.lib.reset_timers() + + # Extract results + rates = self._calculate_reaction_rates(source_rate) + + # Get k and uncertainty + keff = ufloat(*openmc.lib.keff()) + + op_result = OperatorResult(keff, rates) + + return copy.deepcopy(op_result) + + def _update_materials(self): + """Updates material compositions in OpenMC on all processes.""" + + for rank in range(comm.size): + number_i = comm.bcast(self.number, root=rank) + + for mat in number_i.materials: + nuclides = [] + densities = [] + for nuc in number_i.nuclides: + if nuc in self.nuclides_with_data: + val = 1.0e-24 * number_i.get_atom_density(mat, nuc) + + # If nuclide is zero, do not add to the problem. + if val > 0.0: + if self.round_number: + val_magnitude = np.floor(np.log10(val)) + val_scaled = val / 10**val_magnitude + val_round = round(val_scaled, 8) + + val = val_round * 10**val_magnitude + + nuclides.append(nuc) + densities.append(val) + else: + # Only output warnings if values are significantly + # negative. CRAM does not guarantee positive + # values. + if val < -1.0e-21: + print(f'WARNING: nuclide {nuc} in material' + f'{mat} is negative (density = {val}' + + ' atom/b-cm)') + + number_i[mat, nuc] = 0.0 + + # Update densities on C API side + mat_internal = openmc.lib.materials[int(mat)] + mat_internal.set_densities(nuclides, densities) + + # TODO Update densities on the Python side, otherwise the + # summary.h5 file contains densities at the first time step + + @staticmethod + def write_bos_data(step): + """Write a state-point file with beginning of step data + + Parameters + ---------- + step : int + Current depletion step including restarts + + """ + openmc.lib.statepoint_write( + "openmc_simulation_n{}.h5".format(step), + write_source=False) + + def finalize(self): + """Finalize a depletion simulation and release resources.""" + if self.cleanup_when_done: + openmc.lib.finalize() + +# Retain deprecated name for the time being +def Operator(*args, **kwargs): + # warn of name change + warn( + "The Operator(...) class has been renamed and will " + "be removed in a future version of OpenMC. Use " + "CoupledOperator(...) instead.", + FutureWarning + ) + return CoupledOperator(*args, **kwargs) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 33bd182d7..1e382576e 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -1,5 +1,5 @@ """ -Class for normalizing fission energy deposition +Classes for collecting and calculating quantities for reaction rate operators """ import bisect from abc import abstractmethod @@ -68,7 +68,7 @@ class TalliedFissionYieldHelper(FissionYieldHelper): mat_indexes : iterable of int Indices of tallied materials that will have their fission yields computed by this helper. Necessary as the - :class:`openmc.deplete.Operator` that uses this helper + :class:`openmc.deplete.CoupledOperator` that uses this helper may only burn a subset of all materials when running in parallel mode. """ @@ -133,9 +133,11 @@ class DirectReactionRateHelper(ReactionRateHelper): Parameters ---------- n_nucs : int - Number of burnable nuclides tracked by :class:`openmc.deplete.Operator` + Number of burnable nuclides tracked by + :class:`openmc.deplete.CoupledOperator` n_react : int - Number of reactions tracked by :class:`openmc.deplete.Operator` + Number of reactions tracked by an instance of + :class:`openmc.deplete.CoupledOperator` Attributes ---------- @@ -217,9 +219,10 @@ class FluxCollapseHelper(ReactionRateHelper): Parameters ---------- n_nucs : int - Number of burnable nuclides tracked by :class:`openmc.deplete.Operator` + Number of burnable nuclides tracked by + :class:`openmc.deplete.CoupledOperator` n_react : int - Number of reactions tracked by :class:`openmc.deplete.Operator` + Number of reactions tracked by :class:`openmc.deplete.CoupledOperator` energies : iterable of float Energy group boundaries for flux spectrum in [eV] reactions : iterable of str @@ -385,7 +388,7 @@ class ChainFissionHelper(EnergyNormalizationHelper): ---------- nuclides : list of str All nuclides with desired reaction rates. Ordered to be - consistent with :class:`openmc.deplete.Operator` + consistent with :class:`openmc.deplete.CoupledOperator` energy : float Total energy [J/s/source neutron] produced in a transport simulation. Updated in the material iteration with :meth:`update`. @@ -423,7 +426,7 @@ class ChainFissionHelper(EnergyNormalizationHelper): self._fission_q_vector = fission_qs - def update(self, fission_rates): + def update(self, fission_rates, mat_index=None): """Update energy produced with fission rates in a material Parameters @@ -432,6 +435,8 @@ class ChainFissionHelper(EnergyNormalizationHelper): fission reaction rate for each isotope in the specified material. Should be ordered corresponding to initial ``rate_index`` used in :meth:`prepare` + mat_index : int + Unused """ self._energy += dot(fission_rates, self._fission_q_vector) @@ -552,7 +557,7 @@ class ConstantFissionYieldHelper(FissionYieldHelper): Parameters ---------- - operator : openmc.deplete.TransportOperator + operator : openmc.deplete.abc.TransportOperator operator with a depletion chain kwargs: Additional keyword arguments to be used in construction @@ -610,7 +615,6 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): Default: 0.0253 [eV] fast_energy : float, optional Energy of yield data corresponding to fast yields. - Default: 500 [kev] Attributes ---------- @@ -632,10 +636,10 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): Array of fission rate fractions with shape ``(n_mats, 2, n_nucs)``. ``results[:, 0]`` corresponds to the fraction of all fissions - that occured below ``cutoff``. The number + that occurred below ``cutoff``. The number of materials in the first axis corresponds to the number of materials burned by the - :class:`openmc.deplete.Operator` + :class:`openmc.deplete.CoupledOperator` """ def __init__(self, chain_nuclides, n_bmats, cutoff=112.0, @@ -692,7 +696,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): Parameters ---------- - operator : openmc.deplete.Operator + operator : openmc.deplete.CoupledOperator Operator with a chain and burnable materials kwargs: Additional keyword arguments to be used in construction @@ -718,7 +722,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): mat_indexes : iterable of int Indices of tallied materials that will have their fission yields computed by this helper. Necessary as the - :class:`openmc.deplete.Operator` that uses this helper + :class:`openmc.deplete.CoupledOperator` that uses this helper may only burn a subset of all materials when running in parallel mode. """ @@ -788,7 +792,7 @@ class FissionYieldCutoffHelper(TalliedFissionYieldHelper): class AveragedFissionYieldHelper(TalliedFissionYieldHelper): r"""Class that computes fission yields based on average fission energy - Computes average energy at which fission events occured with + Computes average energy at which fission events occurred with .. math:: @@ -843,7 +847,7 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): mat_indexes : iterable of int Indices of tallied materials that will have their fission yields computed by this helper. Necessary as the - :class:`openmc.deplete.Operator` that uses this helper + :class:`openmc.deplete.CoupledOperator` that uses this helper may only burn a subset of all materials when running in parallel mode. """ @@ -906,7 +910,7 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): Use the computed average energy of fission events to determine fission yields. If average energy is between two sets of yields, linearly - interpolate bewteen the two. + interpolate between the two. Otherwise take the closet set of yields. Parameters @@ -956,7 +960,7 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper): Parameters ---------- - operator : openmc.deplete.TransportOperator + operator : openmc.deplete.CoupledOperator Operator with a depletion chain kwargs : Additional keyword arguments to be used in construction diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py new file mode 100644 index 000000000..b7956205a --- /dev/null +++ b/openmc/deplete/independent_operator.py @@ -0,0 +1,452 @@ +"""Transport-independent transport operator for depletion. + +This module implements a transport operator that runs independently of any +transport solver by using user-provided one-group cross sections. + +""" + +import copy +from collections import OrderedDict +from warnings import warn +from itertools import product + +import numpy as np +from uncertainties import ufloat + +import openmc +from openmc.checkvalue import check_type +from openmc.mpi import comm +from .abc import ReactionRateHelper, OperatorResult +from .openmc_operator import OpenMCOperator, _distribute +from .microxs import MicroXS +from .results import Results +from .helpers import ChainFissionHelper, ConstantFissionYieldHelper, SourceRateHelper + +class IndependentOperator(OpenMCOperator): + """Transport-independent transport operator that uses one-group cross + sections to calculate reaction rates. + + Instances of this class can be used to perform depletion using one-group + cross sections and constant flux or constant power. Normally, a user needn't + call methods of this class directly. Instead, an instance of this class is + passed to an integrator class, such as + :class:`openmc.deplete.CECMIntegrator`. + + Parameters + ---------- + materials : openmc.Materials + Materials to deplete. + micro_xs : MicroXS + One-group microscopic cross sections in [b] . + chain_file : str + Path to the depletion chain XML file. + keff : 2-tuple of float, optional + keff eigenvalue and uncertainty from transport calculation. + Default is None. + prev_results : Results, optional + Results from a previous depletion calculation. + normalization_mode : {"fission-q", "source-rate"} + Indicate how reaction rates should be calculated. + ``"fission-q"`` uses the fission Q values from the depletion chain to + compute the flux based on the power. ``"source-rate"`` uses a the + source rate (assumed to be neutron flux) to calculate the + reaction rates. + fission_q : dict, optional + Dictionary of nuclides and their fission Q values [eV]. If not given, + values will be pulled from the ``chain_file``. Only applicable + if ``"normalization_mode" == "fission-q"``. + reduce_chain : bool, optional + If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the + depletion chain up to ``reduce_chain_level``. + reduce_chain_level : int, optional + Depth of the search when reducing the depletion chain. Only used + if ``reduce_chain`` evaluates to true. The default value of + ``None`` implies no limit on the depth. + fission_yield_opts : dict of str to option, optional + Optional arguments to pass to the + :class:`openmc.deplete.helpers.FissionYieldHelper` object. Will be + passed directly on to the helper. Passing a value of None will use + the defaults for the associated helper. + + Attributes + ---------- + materials : openmc.Materials + All materials present in the model + cross_sections : MicroXS + Object containing one-group cross-sections in [cm^2]. + dilute_initial : float + Initial atom density [atoms/cm^3] to add for nuclides that + are zero in initial condition to ensure they exist in the decay + chain. Only done for nuclides with reaction rates. + output_dir : pathlib.Path + Path to output directory to save results. + round_number : bool + Whether or not to round output to OpenMC to 8 digits. + Useful in testing, as OpenMC is incredibly sensitive to exact values. + number : openmc.deplete.AtomNumber + Total number of atoms in simulation. + nuclides_with_data : set of str + A set listing all unique nuclides available from cross_sections.xml. + chain : openmc.deplete.Chain + The depletion chain information necessary to form matrices and tallies. + reaction_rates : openmc.deplete.ReactionRates + Reaction rates from the last operator step. + burnable_mats : list of str + All burnable material IDs + heavy_metal : float + Initial heavy metal inventory [g] + local_mats : list of str + All burnable material IDs being managed by a single process + prev_res : Results or None + Results from a previous depletion calculation. ``None`` if no + results are to be used. + + """ + + def __init__(self, + materials, + micro_xs, + chain_file, + keff=None, + normalization_mode='source-rate', + fission_q=None, + prev_results=None, + reduce_chain=False, + reduce_chain_level=None, + fission_yield_opts=None): + # Validate micro-xs parameters + check_type('materials', materials, openmc.Materials) + check_type('micro_xs', micro_xs, MicroXS) + if keff is not None: + check_type('keff', keff, tuple, float) + keff = ufloat(*keff) + + self._keff = keff + + if fission_yield_opts is None: + fission_yield_opts = {} + helper_kwargs = {'normalization_mode': normalization_mode, + 'fission_yield_opts': fission_yield_opts} + + cross_sections = micro_xs * 1e-24 + super().__init__( + materials, + cross_sections, + chain_file, + prev_results, + fission_q=fission_q, + helper_kwargs=helper_kwargs, + reduce_chain=reduce_chain, + reduce_chain_level=reduce_chain_level) + + @classmethod + def from_nuclides(cls, volume, nuclides, + micro_xs, + chain_file, + nuc_units='atom/b-cm', + keff=None, + normalization_mode='source-rate', + fission_q=None, + prev_results=None, + reduce_chain=False, + reduce_chain_level=None, + fission_yield_opts=None): + """ + Alternate constructor from a dictionary of nuclide concentrations + + volume : float + Volume of the material being depleted in [cm^3] + nuclides : dict of str to float + Dictionary with nuclide names as keys and nuclide concentrations as + values. + micro_xs : MicroXS + One-group microscopic cross sections. + chain_file : str + Path to the depletion chain XML file. + nuc_units : {'atom/cm3', 'atom/b-cm'} + Units for nuclide concentration. + keff : 2-tuple of float, optional + keff eigenvalue and uncertainty from transport calculation. + Default is None. + normalization_mode : {"fission-q", "source-rate"} + Indicate how reaction rates should be calculated. + ``"fission-q"`` uses the fission Q values from the depletion + chain to compute the flux based on the power. ``"source-rate"`` uses + the source rate (assumed to be neutron flux) to calculate the + reaction rates. + fission_q : dict, optional + Dictionary of nuclides and their fission Q values [eV]. If not + given, values will be pulled from the ``chain_file``. Only + applicable if ``"normalization_mode" == "fission-q"``. + prev_results : Results, optional + Results from a previous depletion calculation. + reduce_chain : bool, optional + If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the + depletion chain up to ``reduce_chain_level``. Default is False. + reduce_chain_level : int, optional + Depth of the search when reducing the depletion chain. Only used + if ``reduce_chain`` evaluates to true. The default value of + ``None`` implies no limit on the depth. + fission_yield_opts : dict of str to option, optional + Optional arguments to pass to the + :class:`openmc.deplete.helpers.FissionYieldHelper` class. Will be + passed directly on to the helper. Passing a value of None will use + the defaults for the associated helper. + + """ + check_type('nuclides', nuclides, dict, str) + materials = cls._consolidate_nuclides_to_material(nuclides, nuc_units, volume) + return cls(materials, + micro_xs, + chain_file, + keff=keff, + normalization_mode=normalization_mode, + fission_q=fission_q, + prev_results=prev_results, + reduce_chain=reduce_chain, + reduce_chain_level=reduce_chain_level, + fission_yield_opts=fission_yield_opts) + + + @staticmethod + def _consolidate_nuclides_to_material(nuclides, nuc_units, volume): + """Puts nuclide list into an openmc.Materials object. + + """ + openmc.reset_auto_ids() + mat = openmc.Material() + if nuc_units == 'atom/b-cm': + for nuc, conc in nuclides.items(): + mat.add_nuclide(nuc, conc) + elif nuc_units == 'atom/cm3': + for nuc, conc in nuclides.items(): + mat.add_nuclide(nuc, conc * 1e-24) # convert to at/b-cm + else: + raise ValueError(f"Unit '{nuc_units}' is invalid.") + + mat.volume = volume + mat.depletable = True + + return openmc.Materials([mat]) + + def _load_previous_results(self): + """Load results from a previous depletion simulation""" + # Reload volumes into geometry + model = openmc.Model(materials=self.materials) + self.prev_res[-1].transfer_volumes(model) + self.materials = model.materials + + + # Store previous results in operator + # Distribute reaction rates according to those tracked + # on this process + if comm.size != 1: + prev_results = self.prev_res + self.prev_res = Results() + mat_indexes = _distribute(range(len(self.burnable_mats))) + for res_obj in prev_results: + new_res = res_obj.distribute(self.local_mats, mat_indexes) + self.prev_res.append(new_res) + + + def _get_nuclides_with_data(self, cross_sections): + """Finds nuclides with cross section data""" + return set(cross_sections.index) + + class _IndependentNormalizationHelper(ChainFissionHelper): + """Class for calculating one-group flux based on a power. + + flux = Power / X, where X = volume * sum_i(Q_i * fission_micro_xs_i * density_i) + + Parameters + ---------- + op : openmc.deplete.IndependentOperator + Reference to the object encapsulating _IndependentNormalizationHelper. + We pass this so we don't have to duplicate :attr:`IndependentOperator.number`. + + """ + + def __init__(self, op): + rates = op.reaction_rates + self.nuc_ind_map = {ind: nuc for nuc, ind in rates.index_nuc.items()} + self._op = op + super().__init__() + + def update(self, fission_rates, mat_index=None): + """Update 'energy' produced with fission rates in a material. What + this actually calculates is the quantity X. + + Parameters + ---------- + fission_rates : numpy.ndarray + fission reaction rate for each isotope in the specified + material. Should be ordered corresponding to initial + ``rate_index`` used in :meth:`prepare` + mat_index : int + Material index + + """ + volume = self._op.number.get_mat_volume(mat_index) + densities = np.empty_like(fission_rates) + for i_nuc, nuc in self.nuc_ind_map.items(): + densities[i_nuc] = self._op.number.get_atom_density(mat_index, nuc) + fission_rates *= volume * densities + + super().update(fission_rates) + + class _IndependentRateHelper(ReactionRateHelper): + """Class for generating one-group reaction rates with flux and + one-group cross sections. + + This class does not generate tallies, and instead stores cross sections + for each nuclide and transmutation reaction relevant for a depletion + calculation. The reaction rate is calculated by multiplying the flux by the + cross sections. + + Parameters + ---------- + op : openmc.deplete.IndependentOperator + Reference to the object encapsulate _IndependentRateHelper. + We pass this so we don't have to duplicate the :attr:`IndependentOperator.number` object. + + + Attributes + ---------- + nuc_ind_map : dict of int to str + Dictionary mapping the nuclide index to nuclide name + rxn_ind_map : dict of int to str + Dictionary mapping reaction index to reaction name + + """ + + def __init__(self, op): + rates = op.reaction_rates + super().__init__(rates.n_nuc, rates.n_react) + + self.nuc_ind_map = {ind: nuc for nuc, ind in rates.index_nuc.items()} + self.rxn_ind_map = {ind: rxn for rxn, ind in rates.index_rx.items()} + self._op = op + + def generate_tallies(self, materials, scores): + """Unused in this case""" + pass + + def get_material_rates(self, mat_id, nuc_index, react_index): + """Return 2D array of [nuclide, reaction] reaction rates + + Parameters + ---------- + mat_id : int + Unique ID for the requested material + nuc_index : list of str + Ordering of desired nuclides + react_index : list of str + Ordering of reactions + """ + self._results_cache.fill(0.0) + + volume = self._op.number.get_mat_volume(mat_id) + for i_nuc, i_react in product(nuc_index, react_index): + nuc = self.nuc_ind_map[i_nuc] + rxn = self.rxn_ind_map[i_react] + density = self._op.number.get_atom_density(mat_id, nuc) + self._results_cache[i_nuc, + i_react] = self._op.cross_sections[rxn][nuc] * density * volume + + return self._results_cache + + def _get_helper_classes(self, helper_kwargs): + """Get helper classes for calculating reaction rates and fission yields + + Parameters + ---------- + helper_kwargs : dict + Keyword arguments for helper classes + + """ + + normalization_mode = helper_kwargs['normalization_mode'] + fission_yield_opts = helper_kwargs['fission_yield_opts'] + + self._rate_helper = self._IndependentRateHelper(self) + if normalization_mode == "fission-q": + self._normalization_helper = self._IndependentNormalizationHelper(self) + else: + self._normalization_helper = SourceRateHelper() + + # Select and create fission yield helper + fission_helper = ConstantFissionYieldHelper + self._yield_helper = fission_helper.from_operator( + self, **fission_yield_opts) + + def initial_condition(self): + """Performs final setup and returns initial condition. + + Returns + ------- + list of numpy.ndarray + Total density for initial conditions. + """ + + # Return number density vector + return super().initial_condition(self.materials) + + def __call__(self, vec, source_rate): + """Obtain the reaction rates + + Parameters + ---------- + vec : list of numpy.ndarray + Total atoms to be used in function. + source_rate : float + Power in [W] or flux in [neutron/cm^2-s] + + Returns + ------- + openmc.deplete.OperatorResult + Eigenvalue and reaction rates resulting from transport operator + + """ + + self._update_materials_and_nuclides(vec) + + rates = self._calculate_reaction_rates(source_rate) + keff = self._keff + + op_result = OperatorResult(keff, rates) + return copy.deepcopy(op_result) + + def _update_materials(self): + """Updates material compositions in OpenMC on all processes.""" + + for rank in range(comm.size): + number_i = comm.bcast(self.number, root=rank) + + for mat in number_i.materials: + nuclides = [] + densities = [] + for nuc in number_i.nuclides: + if nuc in self.nuclides_with_data: + val = 1.0e-24 * number_i.get_atom_density(mat, nuc) + + # If nuclide is zero, do not add to the problem. + if val > 0.0: + if self.round_number: + val_magnitude = np.floor(np.log10(val)) + val_scaled = val / 10**val_magnitude + val_round = round(val_scaled, 8) + + val = val_round * 10**val_magnitude + + nuclides.append(nuc) + densities.append(val) + else: + # Only output warnings if values are significantly + # negative. CRAM does not guarantee positive + # values. + if val < -1.0e-21: + print(f'WARNING: nuclide {nuc} in material' + f'{mat} is negative (density = {val}' + + ' atom/b-cm)') + number_i[mat, nuc] = 0.0 diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py index c19ef076a..74a3cebdb 100644 --- a/openmc/deplete/integrators.py +++ b/openmc/deplete/integrators.py @@ -103,7 +103,7 @@ class CECMIntegrator(Integrator): op_results : list of openmc.deplete.OperatorResult Eigenvalue and reaction rates from transport simulations """ - # deplete across first half of inteval + # deplete across first half of interval time0, x_middle = self._timed_deplete(conc, rates, dt / 2) res_middle = self.operator(x_middle, source_rate) diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py new file mode 100644 index 000000000..c4094970f --- /dev/null +++ b/openmc/deplete/microxs.py @@ -0,0 +1,230 @@ +"""MicroXS module + +A pandas.DataFrame storing microscopic cross section data with +nuclides names as row indices and reaction names as column indices. +""" + +import tempfile +from pathlib import Path +from copy import deepcopy + +from pandas import DataFrame, read_csv, concat +import numpy as np + +from openmc.checkvalue import check_type, check_value, check_iterable_type +from openmc.mgxs import EnergyGroups, ArbitraryXS, FissionXS +from openmc.data import DataLibrary +from openmc import Tallies, StatePoint, Materials, Material + + +from .chain import Chain, REACTIONS +from .coupled_operator import _find_cross_sections, _get_nuclides_with_data + +_valid_rxns = list(REACTIONS) +_valid_rxns.append('fission') + + +class MicroXS(DataFrame): + """Stores microscopic cross section data for use in + independent depletion. + """ + + @classmethod + def from_model(cls, + model, + reaction_domain, + chain_file, + dilute_initial=1.0e3, + energy_bounds=(0, 20e6)): + """Generate a one-group cross-section dataframe using + OpenMC. Note that the ``openmc`` executable must be compiled. + + Parameters + ---------- + model : openmc.Model + OpenMC model object. Must contain geometry, materials, and settings. + reaction_domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.RegularMesh + Domain in which to tally reaction rates. + chain_file : str + Path to the depletion chain XML file that will be used in depletion + simulation. Used to determine cross sections for materials not + present in the inital composition. + dilute_initial : float + Initial atom density [atoms/cm^3] to add for nuclides that + are zero in initial condition to ensure they exist in the cross + section data. Only done for nuclides with reaction rates. + reactions : list of str, optional + Reaction names to tally + energy_bound : 2-tuple of float, optional + Bounds for the energy group. + + Returns + ------- + MicroXS + Cross section data in [b] + + """ + groups = EnergyGroups(energy_bounds) + + # Set up the reaction tallies + original_tallies = model.tallies + original_materials = deepcopy(model.materials) + tallies = Tallies() + xs = {} + reactions, diluted_materials = cls._add_dilute_nuclides(chain_file, + model, + dilute_initial) + model.materials = diluted_materials + + for rxn in reactions: + if rxn == 'fission': + xs[rxn] = FissionXS(domain=reaction_domain, groups=groups, by_nuclide=True) + else: + xs[rxn] = ArbitraryXS(rxn, domain=reaction_domain, groups=groups, by_nuclide=True) + tallies += xs[rxn].tallies.values() + + model.tallies = tallies + + # create temporary run + with tempfile.TemporaryDirectory() as temp_dir: + statepoint_path = model.run(cwd=temp_dir) + + with StatePoint(statepoint_path) as sp: + for rxn in xs: + xs[rxn].load_from_statepoint(sp) + + # Build the DataFrame + series = {} + for rx in xs: + df = xs[rx].get_pandas_dataframe(xs_type='micro') + series[rx] = df.set_index('nuclide')['mean'] + + # Revert to the original tallies and materials + model.tallies = original_tallies + model.materials = original_materials + + return cls(series) + + @classmethod + def _add_dilute_nuclides(cls, chain_file, model, dilute_initial): + """ + Add nuclides not present in burnable materials that have neutron data + and are present in the depletion chain to those materials. This allows + us to tally those specific nuclides for reactions to create one-group + cross sections. + + Parameters + ---------- + chain_file : str + Path to the depletion chain XML file that will be used in depletion + simulation. Used to determine cross sections for materials not + present in the inital composition. + model : openmc.Model + Model object + dilute_initial : float + Initial atom density [atoms/cm^3] to add for nuclides that + are zero in initial condition to ensure they exist in the cross + section data. Only done for nuclides with reaction rates. + + Returns + ------- + reactions : list of str + List of reaction names + diluted_materials : openmc.Materials + :class:`openmc.Materials` object with nuclides added to burnable + materials. + """ + chain = Chain.from_xml(chain_file) + reactions = chain.reactions + cross_sections = _find_cross_sections(model) + nuclides_with_data = _get_nuclides_with_data(cross_sections) + burnable_nucs = [nuc.name for nuc in chain.nuclides + if nuc.name in nuclides_with_data] + diluted_materials = Materials() + for material in model.materials: + if material.depletable: + nuc_densities = material.get_nuclide_atom_densities() + dilute_density = 1.e-24 * dilute_initial + material.set_density('sum') + for nuc, density in nuc_densities.items(): + material.remove_nuclide(nuc) + material.add_nuclide(nuc, density) + for burn_nuc in burnable_nucs: + if burn_nuc not in nuc_densities: + material.add_nuclide(burn_nuc, + dilute_density) + diluted_materials.append(material) + + return reactions, diluted_materials + + @classmethod + def from_array(cls, nuclides, reactions, data): + """ + Creates a ``MicroXS`` object from arrays. + + Parameters + ---------- + nuclides : list of str + List of nuclide symbols for that have data for at least one + reaction. + reactions : list of str + List of reactions. All reactions must match those in + :data:`openmc.deplete.chain.REACTIONS` + data : ndarray of floats + Array containing one-group microscopic cross section values for + each nuclide and reaction. Cross section values are assumed to be + in [b]. + + Returns + ------- + MicroXS + """ + + # Validate inputs + if data.shape != (len(nuclides), len(reactions)): + raise ValueError( + f'Nuclides list of length {len(nuclides)} and ' + f'reactions array of length {len(reactions)} do not ' + f'match dimensions of data array of shape {data.shape}') + + cls._validate_micro_xs_inputs( + nuclides, reactions, data) + micro_xs = cls(index=nuclides, columns=reactions, data=data) + + return micro_xs + + @classmethod + def from_csv(cls, csv_file, **kwargs): + """ + Load a ``MicroXS`` object from a ``.csv`` file. + + Parameters + ---------- + csv_file : str + Relative path to csv-file containing microscopic cross section + data. Cross section values are assumed to be in [b] + **kwargs : dict + Keyword arguments to pass to :func:`pandas.read_csv()`. + + Returns + ------- + MicroXS + + """ + if 'float_precision' not in kwargs: + kwargs['float_precision'] = 'round_trip' + + micro_xs = cls(read_csv(csv_file, index_col=0, **kwargs)) + + cls._validate_micro_xs_inputs(list(micro_xs.index), + list(micro_xs.columns), + micro_xs.to_numpy()) + return micro_xs + + @staticmethod + def _validate_micro_xs_inputs(nuclides, reactions, data): + check_iterable_type('nuclides', nuclides, str) + check_iterable_type('reactions', reactions, str) + check_type('data', data, np.ndarray, expected_iter_type=float) + for reaction in reactions: + check_value('reactions', reaction, _valid_rxns) diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py new file mode 100644 index 000000000..6e6f8dd8d --- /dev/null +++ b/openmc/deplete/openmc_operator.py @@ -0,0 +1,571 @@ +"""OpenMC transport operator + +This module implements functions shared by both OpenMC transport-coupled and +transport-independent transport operators. + +""" + +from abc import abstractmethod +from collections import OrderedDict + +import numpy as np + +import openmc +from openmc.mpi import comm +from .abc import TransportOperator, OperatorResult +from .atom_number import AtomNumber +from .reaction_rates import ReactionRates + +__all__ = ["OpenMCOperator", "OperatorResult"] + + +def _distribute(items): + """Distribute items across MPI communicator + + Parameters + ---------- + items : list + List of items of distribute + + Returns + ------- + list + Items assigned to process that called + + """ + min_size, extra = divmod(len(items), comm.size) + j = 0 + for i in range(comm.size): + chunk_size = min_size + int(i < extra) + if comm.rank == i: + return items[j:j + chunk_size] + j += chunk_size + + +class OpenMCOperator(TransportOperator): + """Abstract class holding OpenMC-specific functions for running + depletion calculations. + + Specific classes for running transport-coupled or transport-independent + depletion calculations are implemented as subclasses of OpenMCOperator. + + Parameters + ---------- + materials : openmc.Materials + List of all materials in the model + cross_sections : str or pandas.DataFrame + Path to continuous energy cross section library, or object containing + one-group cross-sections. + chain_file : str, optional + Path to the depletion chain XML file. Defaults to the file + listed under ``depletion_chain`` in + :envvar:`OPENMC_CROSS_SECTIONS` environment variable. + prev_results : Results, optional + Results from a previous depletion calculation. If this argument is + specified, the depletion calculation will start from the latest state + in the previous results. + diff_burnable_mats : bool, optional + Whether to differentiate burnable materials with multiple instances. + Volumes are divided equally from the original material volume. + fission_q : dict, optional + Dictionary of nuclides and their fission Q values [eV]. + dilute_initial : float, optional + Initial atom density [atoms/cm^3] to add for nuclides that are zero + in initial condition to ensure they exist in the decay chain. + Only done for nuclides with reaction rates. + helper_kwargs : dict + Keyword arguments for helper classes + reduce_chain : bool, optional + If True, use :meth:`openmc.deplete.Chain.reduce()` to reduce the + depletion chain up to ``reduce_chain_level``. + reduce_chain_level : int, optional + Depth of the search when reducing the depletion chain. Only used + if ``reduce_chain`` evaluates to true. The default value of + ``None`` implies no limit on the depth. + + + Attributes + ---------- + materials : openmc.Materials + All materials present in the model + cross_sections : str or MicroXS + Path to continuous energy cross section library, or object + containing one-group cross-sections. + dilute_initial : float + Initial atom density [atoms/cm^3] to add for nuclides that + are zero in initial condition to ensure they exist in the decay + chain. Only done for nuclides with reaction rates. + output_dir : pathlib.Path + Path to output directory to save results. + round_number : bool + Whether or not to round output to OpenMC to 8 digits. + Useful in testing, as OpenMC is incredibly sensitive to exact values. + number : openmc.deplete.AtomNumber + Total number of atoms in simulation. + nuclides_with_data : set of str + A set listing all unique nuclides available from cross_sections.xml. + chain : openmc.deplete.Chain + The depletion chain information necessary to form matrices and tallies. + reaction_rates : openmc.deplete.ReactionRates + Reaction rates from the last operator step. + burnable_mats : list of str + All burnable material IDs + heavy_metal : float + Initial heavy metal inventory [g] + local_mats : list of str + All burnable material IDs being managed by a single process + prev_res : Results or None + Results from a previous depletion calculation. ``None`` if no + results are to be used. + + """ + + def __init__( + self, + materials=None, + cross_sections=None, + chain_file=None, + prev_results=None, + diff_burnable_mats=False, + fission_q=None, + dilute_initial=0.0, + helper_kwargs=None, + reduce_chain=False, + reduce_chain_level=None): + + super().__init__(chain_file, fission_q, dilute_initial, prev_results) + self.round_number = False + self.materials = materials + self.cross_sections = cross_sections + + # Reduce the chain to only those nuclides present + if reduce_chain: + init_nuclides = set() + for material in self.materials: + if not material.depletable: + continue + for name, _dens_percent, _dens_type in material.nuclides: + init_nuclides.add(name) + + self.chain = self.chain.reduce(init_nuclides, reduce_chain_level) + + if diff_burnable_mats: + self._differentiate_burnable_mats() + + # Determine which nuclides have cross section data + # This nuclides variables contains every nuclides + # for which there is an entry in the micro_xs parameter + openmc.reset_auto_ids() + self.burnable_mats, volumes, all_nuclides = self._get_burnable_mats() + self.local_mats = _distribute(self.burnable_mats) + + self._mat_index_map = { + lm: self.burnable_mats.index(lm) for lm in self.local_mats} + + if self.prev_res is not None: + self._load_previous_results() + + self.nuclides_with_data = self._get_nuclides_with_data( + self.cross_sections) + + # Select nuclides with data that are also in the chain + self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides + if nuc.name in self.nuclides_with_data] + + # Extract number densities from the geometry / previous depletion run + self._extract_number(self.local_mats, + volumes, + all_nuclides, + self.prev_res) + + # Create reaction rates array + self.reaction_rates = ReactionRates( + self.local_mats, self._burnable_nucs, self.chain.reactions) + + self._get_helper_classes(helper_kwargs) + + def _differentiate_burnable_mats(self): + """Assign distribmats for each burnable material""" + pass + + def _get_burnable_mats(self): + """Determine depletable materials, volumes, and nuclides + + Returns + ------- + burnable_mats : list of str + List of burnable material IDs + volume : OrderedDict of str to float + Volume of each material in [cm^3] + nuclides : list of str + Nuclides in order of how they'll appear in the simulation. + + """ + + burnable_mats = set() + model_nuclides = set() + volume = OrderedDict() + + self.heavy_metal = 0.0 + + # Iterate once through the geometry to get dictionaries + for mat in self.materials: + for nuclide in mat.get_nuclides(): + model_nuclides.add(nuclide) + if mat.depletable: + burnable_mats.add(str(mat.id)) + if mat.volume is None: + raise RuntimeError("Volume not specified for depletable " + "material with ID={}.".format(mat.id)) + volume[str(mat.id)] = mat.volume + self.heavy_metal += mat.fissionable_mass + + # Make sure there are burnable materials + if not burnable_mats: + raise RuntimeError( + "No depletable materials were found in the model.") + + # Sort the sets + burnable_mats = sorted(burnable_mats, key=int) + model_nuclides = sorted(model_nuclides) + + # Construct a global nuclide dictionary, burned first + nuclides = list(self.chain.nuclide_dict) + for nuc in model_nuclides: + if nuc not in nuclides: + nuclides.append(nuc) + + return burnable_mats, volume, nuclides + + def _load_previous_results(self): + """Load results from a previous depletion simulation""" + pass + + @abstractmethod + def _get_nuclides_with_data(self, cross_sections): + """Find nuclides with cross section data + + Parameters + ---------- + cross_sections : str or pandas.DataFrame + Path to continuous energy cross section library, or object + containing one-group cross-sections. + + Returns + ------- + nuclides : set of str + Set of nuclide names that have cross secton data + + """ + + def _extract_number(self, local_mats, volume, all_nuclides, prev_res=None): + """Construct AtomNumber using geometry + + Parameters + ---------- + local_mats : list of str + Material IDs to be managed by this process + volume : OrderedDict of str to float + Volumes for the above materials in [cm^3] + all_nuclides : list of str + Nuclides to be used in the simulation. + prev_res : Results, optional + Results from a previous depletion calculation + + """ + self.number = AtomNumber(local_mats, all_nuclides, volume, len(self.chain)) + + if self.dilute_initial != 0.0: + for nuc in self._burnable_nucs: + self.number.set_atom_density( + np.s_[:], nuc, self.dilute_initial) + + # Now extract and store the number densities + # From the geometry if no previous depletion results + if prev_res is None: + for mat in self.materials: + if str(mat.id) in local_mats: + self._set_number_from_mat(mat) + + # Else from previous depletion results + else: + for mat in self.materials: + if str(mat.id) in local_mats: + self._set_number_from_results(mat, prev_res) + + def _set_number_from_mat(self, mat): + """Extracts material and number densities from openmc.Material + + Parameters + ---------- + mat : openmc.Material + The material to read from + + """ + mat_id = str(mat.id) + + for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items(): + atom_per_cc = atom_per_bcm * 1.0e24 + self.number.set_atom_density(mat_id, nuclide, atom_per_cc) + + def _set_number_from_results(self, mat, prev_res): + """Extracts material nuclides and number densities. + + If the nuclide concentration's evolution is tracked, the densities come + from depletion results. Else, densities are extracted from the geometry + in the summary. + + Parameters + ---------- + mat : openmc.Material + The material to read from + prev_res : Results + Results from a previous depletion calculation + + """ + mat_id = str(mat.id) + + # Get nuclide lists from geometry and depletion results + depl_nuc = prev_res[-1].nuc_to_ind + geom_nuc_densities = mat.get_nuclide_atom_densities() + + # Merge lists of nuclides, with the same order for every calculation + geom_nuc_densities.update(depl_nuc) + + for nuclide, atom_per_bcm in geom_nuc_densities.items(): + if nuclide in depl_nuc: + concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1] + volume = prev_res[-1].volume[mat_id] + atom_per_cc = concentration / volume + else: + atom_per_cc = atom_per_bcm * 1.0e24 + + self.number.set_atom_density(mat_id, nuclide, atom_per_cc) + + @abstractmethod + def _get_helper_classes(self, helper_kwargs): + """Create the ``_rate_helper``, ``_normalization_helper``, and + ``_yield_helper`` objects. + + Parameters + ---------- + helper_kwargs : dict + Keyword arguments for helper classes + + """ + + def initial_condition(self, materials): + """Performs final setup and returns initial condition. + + Parameters + ---------- + materials : list of str + list of material IDs + + Returns + ------- + list of numpy.ndarray + Total density for initial conditions. + + """ + + self._rate_helper.generate_tallies(materials, self.chain.reactions) + self._normalization_helper.prepare( + self.chain.nuclides, self.reaction_rates.index_nuc) + # Tell fission yield helper what materials this process is + # responsible for + self._yield_helper.generate_tallies( + materials, tuple(sorted(self._mat_index_map.values()))) + + # Return number density vector + return list(self.number.get_mat_slice(np.s_[:])) + + def _update_materials_and_nuclides(self, vec): + """Update the number density, material compositions, and nuclide + lists in helper objects + + Parameters + ---------- + vec : list of numpy.ndarray + Total atoms. + + """ + + # Update the number densities regardless of the source rate + self.number.set_density(vec) + self._update_materials() + + # Prevent OpenMC from complaining about re-creating tallies + openmc.reset_auto_ids() + + # Update tally nuclides data in preparation for transport solve + nuclides = self._get_reaction_nuclides() + self._rate_helper.nuclides = nuclides + self._normalization_helper.nuclides = nuclides + self._yield_helper.update_tally_nuclides(nuclides) + + @abstractmethod + def _update_materials(self): + """Updates material compositions in OpenMC on all processes.""" + + def write_bos_data(self, step): + """Document beginning of step data for a given step + + Called at the beginning of a depletion step and at + the final point in the simulation. + + Parameters + ---------- + step : int + Current depletion step including restarts + """ + # Since we aren't running a transport simulation, we simply pass + pass + + def _get_reaction_nuclides(self): + """Determine nuclides that should be tallied for reaction rates. + + This method returns a list of all nuclides that have cross section data + and are listed in the depletion chain. Technically, we should count + nuclides that may not appear in the depletion chain because we still + need to get the fission reaction rate for these nuclides in order to + normalize power, but that is left as a future exercise. + + Returns + ------- + list of str + Nuclides with reaction rates + + """ + nuc_set = set() + + # Create the set of all nuclides in the decay chain in materials marked + # for burning in which the number density is greater than zero. + for nuc in self.number.nuclides: + if nuc in self.nuclides_with_data: + if np.sum(self.number[:, nuc]) > 0.0: + nuc_set.add(nuc) + + # Communicate which nuclides have nonzeros to rank 0 + if comm.rank == 0: + for i in range(1, comm.size): + nuc_newset = comm.recv(source=i, tag=i) + nuc_set |= nuc_newset + else: + comm.send(nuc_set, dest=0, tag=comm.rank) + + if comm.rank == 0: + # Sort nuclides in the same order as self.number + nuc_list = [nuc for nuc in self.number.nuclides + if nuc in nuc_set] + else: + nuc_list = None + + # Store list of nuclides on each process + nuc_list = comm.bcast(nuc_list) + return [nuc for nuc in nuc_list if nuc in self.chain] + + def _calculate_reaction_rates(self, source_rate): + """Unpack tallies from OpenMC and return an operator result + + This method uses OpenMC's C API bindings to determine the k-effective + value and reaction rates from the simulation. The reaction rates are + normalized by a helper class depending on the method being used. + + Parameters + ---------- + source_rate : float + Power in [W] or source rate in [neutron/sec] + + Returns + ------- + rates : openmc.deplete.ReactionRates + Reaction rates for nuclides + + """ + rates = self.reaction_rates + rates.fill(0.0) + + # Extract reaction nuclides + rxn_nuclides = self._rate_helper.nuclides + + # Form fast map + nuc_ind = [rates.index_nuc[nuc] for nuc in rxn_nuclides] + react_ind = [rates.index_rx[react] for react in self.chain.reactions] + + # Keep track of energy produced from all reactions in eV per source + # particle + self._normalization_helper.reset() + self._yield_helper.unpack() + + # Store fission yield dictionaries + fission_yields = [] + + # Create arrays to store fission Q values, reaction rates, and nuclide + # numbers, zeroed out in material iteration + number = np.empty(rates.n_nuc) + + fission_ind = rates.index_rx.get("fission") + + # Extract results + for i, mat in enumerate(self.local_mats): + # Get tally index + mat_index = self._mat_index_map[mat] + + # Zero out reaction rates and nuclide numbers + number.fill(0.0) + + # Get new number densities + for nuc, i_nuc_results in zip(rxn_nuclides, nuc_ind): + number[i_nuc_results] = self.number[mat, nuc] + + tally_rates = self._rate_helper.get_material_rates( + mat_index, nuc_ind, react_ind) + + # Compute fission yields for this material + fission_yields.append(self._yield_helper.weighted_yields(i)) + + # Accumulate energy from fission + if fission_ind is not None: + self._normalization_helper.update( + tally_rates[:, fission_ind], + mat_index=mat_index) + + # Divide by total number and store + rates[i] = self._rate_helper.divide_by_adens(number) + + # Scale reaction rates to obtain units of reactions/sec + rates *= self._normalization_helper.factor(source_rate) + + # Store new fission yields on the chain + self.chain.fission_yields = fission_yields + + return rates + + def get_results_info(self): + """Returns volume list, material lists, and nuc lists. + + Returns + ------- + volume : dict of str float + Volumes corresponding to materials in full_burn_dict + nuc_list : list of str + A list of all nuclide names. Used for sorting the simulation. + burn_list : list of int + A list of all material IDs to be burned. Used for sorting the simulation. + full_burn_list : list + List of all burnable material IDs + + """ + nuc_list = self.number.burnable_nuclides + burn_list = self.local_mats + + volume = {} + for i, mat in enumerate(burn_list): + volume[mat] = self.number.volume[i] + + # Combine volume dictionaries across processes + volume_list = comm.allgather(volume) + volume = {k: v for d in volume_list for k, v in d.items()} + + return volume, nuc_list, burn_list, self.burnable_mats diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py deleted file mode 100644 index 7ef874e87..000000000 --- a/openmc/deplete/operator.py +++ /dev/null @@ -1,822 +0,0 @@ -"""OpenMC transport operator - -This module implements a transport operator for OpenMC so that it can be used by -depletion integrators. The implementation makes use of the Python bindings to -OpenMC's C API so that reading tally results and updating material number -densities is all done in-memory instead of through the filesystem. - -""" - -import copy -from collections import OrderedDict -import os -from warnings import warn - -import numpy as np -from uncertainties import ufloat - -import openmc -from openmc.checkvalue import check_value -from openmc.data import DataLibrary -from openmc.exceptions import DataError -import openmc.lib -from openmc.mpi import comm -from .abc import TransportOperator, OperatorResult -from .atom_number import AtomNumber -from .chain import _find_chain_file -from .reaction_rates import ReactionRates -from .results import Results -from .helpers import ( - DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper, - FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper, - SourceRateHelper, FluxCollapseHelper) - - -__all__ = ["Operator", "OperatorResult"] - - -def _distribute(items): - """Distribute items across MPI communicator - - Parameters - ---------- - items : list - List of items of distribute - - Returns - ------- - list - Items assigned to process that called - - """ - min_size, extra = divmod(len(items), comm.size) - j = 0 - for i in range(comm.size): - chunk_size = min_size + int(i < extra) - if comm.rank == i: - return items[j:j + chunk_size] - j += chunk_size - - -def _find_cross_sections(model): - """Determine cross sections to use for depletion""" - if model.materials and model.materials.cross_sections is not None: - # Prefer info from Model class if available - return model.materials.cross_sections - - # otherwise fallback to environment variable - cross_sections = os.environ.get("OPENMC_CROSS_SECTIONS") - if cross_sections is None: - raise DataError( - "Cross sections were not specified in Model.materials and " - "the OPENMC_CROSS_SECTIONS environment variable is not set." - ) - return cross_sections - - -class Operator(TransportOperator): - """OpenMC transport operator for depletion. - - Instances of this class can be used to perform depletion using OpenMC as the - transport operator. Normally, a user needn't call methods of this class - directly. Instead, an instance of this class is passed to an integrator - class, such as :class:`openmc.deplete.CECMIntegrator`. - - .. versionchanged:: 0.13.0 - The geometry and settings parameters have been replaced with a - model parameter that takes a :class:`~openmc.model.Model` object - - Parameters - ---------- - model : openmc.model.Model - OpenMC model object - chain_file : str, optional - Path to the depletion chain XML file. Defaults to the file - listed under ``depletion_chain`` in - :envvar:`OPENMC_CROSS_SECTIONS` environment variable. - prev_results : Results, optional - Results from a previous depletion calculation. If this argument is - specified, the depletion calculation will start from the latest state - in the previous results. - diff_burnable_mats : bool, optional - Whether to differentiate burnable materials with multiple instances. - Volumes are divided equally from the original material volume. - Default: False. - normalization_mode : {"energy-deposition", "fission-q", "source-rate"} - Indicate how tally results should be normalized. ``"energy-deposition"`` - computes the total energy deposited in the system and uses the ratio of - the power to the energy produced as a normalization factor. - ``"fission-q"`` uses the fission Q values from the depletion chain to - compute the total energy deposited. ``"source-rate"`` normalizes - tallies based on the source rate (for fixed source calculations). - fission_q : dict, optional - Dictionary of nuclides and their fission Q values [eV]. If not given, - values will be pulled from the ``chain_file``. Only applicable - if ``"normalization_mode" == "fission-q"`` - dilute_initial : float, optional - Initial atom density [atoms/cm^3] to add for nuclides that are zero - in initial condition to ensure they exist in the decay chain. - Only done for nuclides with reaction rates. - Defaults to 1.0e3. - fission_yield_mode : {"constant", "cutoff", "average"} - Key indicating what fission product yield scheme to use. The - key determines what fission energy helper is used: - - * "constant": :class:`~openmc.deplete.helpers.ConstantFissionYieldHelper` - * "cutoff": :class:`~openmc.deplete.helpers.FissionYieldCutoffHelper` - * "average": :class:`~openmc.deplete.helpers.AveragedFissionYieldHelper` - - The documentation on these classes describe their methodology - and differences. Default: ``"constant"`` - fission_yield_opts : dict of str to option, optional - Optional arguments to pass to the helper determined by - ``fission_yield_mode``. Will be passed directly on to the - helper. Passing a value of None will use the defaults for - the associated helper. - reaction_rate_mode : {"direct", "flux"}, optional - Indicate how one-group reaction rates should be calculated. The "direct" - method tallies transmutation reaction rates directly. The "flux" method - tallies a multigroup flux spectrum and then collapses one-group reaction - rates after a transport solve (with an option to tally some reaction - rates directly). - - .. versionadded:: 0.12.1 - reaction_rate_opts : dict, optional - Keyword arguments that are passed to the reaction rate helper class. - When ``reaction_rate_mode`` is set to "flux", energy group boundaries - can be set using the "energies" key. See the - :class:`~openmc.deplete.helpers.FluxCollapseHelper` class for all - options. - - .. versionadded:: 0.12.1 - reduce_chain : bool, optional - If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the - depletion chain up to ``reduce_chain_level``. Default is False. - - .. versionadded:: 0.12 - reduce_chain_level : int, optional - Depth of the search when reducing the depletion chain. Only used - if ``reduce_chain`` evaluates to true. The default value of - ``None`` implies no limit on the depth. - - .. versionadded:: 0.12 - - Attributes - ---------- - model : openmc.model.Model - OpenMC model object - geometry : openmc.Geometry - OpenMC geometry object - settings : openmc.Settings - OpenMC settings object - dilute_initial : float - Initial atom density [atoms/cm^3] to add for nuclides that - are zero in initial condition to ensure they exist in the decay - chain. Only done for nuclides with reaction rates. - output_dir : pathlib.Path - Path to output directory to save results. - round_number : bool - Whether or not to round output to OpenMC to 8 digits. - Useful in testing, as OpenMC is incredibly sensitive to exact values. - number : openmc.deplete.AtomNumber - Total number of atoms in simulation. - nuclides_with_data : set of str - A set listing all unique nuclides available from cross_sections.xml. - chain : openmc.deplete.Chain - The depletion chain information necessary to form matrices and tallies. - reaction_rates : openmc.deplete.ReactionRates - Reaction rates from the last operator step. - burnable_mats : list of str - All burnable material IDs - heavy_metal : float - Initial heavy metal inventory [g] - local_mats : list of str - All burnable material IDs being managed by a single process - prev_res : Results or None - Results from a previous depletion calculation. ``None`` if no - results are to be used. - cleanup_when_done : bool - Whether to finalize and clear the shared library memory when the - depletion operation is complete. Defaults to clearing the library. - """ - _fission_helpers = { - "average": AveragedFissionYieldHelper, - "constant": ConstantFissionYieldHelper, - "cutoff": FissionYieldCutoffHelper, - } - - def __init__(self, model, chain_file=None, prev_results=None, - diff_burnable_mats=False, normalization_mode="fission-q", - fission_q=None, dilute_initial=1.0e3, - fission_yield_mode="constant", fission_yield_opts=None, - reaction_rate_mode="direct", reaction_rate_opts=None, - reduce_chain=False, reduce_chain_level=None): - # check for old call to constructor - if isinstance(model, openmc.Geometry): - msg = "As of version 0.13.0 openmc.deplete.Operator requires an " \ - "openmc.Model object rather than the openmc.Geometry and " \ - "openmc.Settings parameters. Please use the geometry and " \ - "settings objects passed here to create a model with which " \ - "to generate the depletion Operator." - raise TypeError(msg) - - # Determine cross sections / depletion chain - cross_sections = _find_cross_sections(model) - if chain_file is None: - chain_file = _find_chain_file(cross_sections) - - check_value('fission yield mode', fission_yield_mode, - self._fission_helpers.keys()) - check_value('normalization mode', normalization_mode, - ('energy-deposition', 'fission-q', 'source-rate')) - if normalization_mode != "fission-q": - if fission_q is not None: - warn("Fission Q dictionary will not be used") - fission_q = None - super().__init__(chain_file, fission_q, dilute_initial, prev_results) - self.round_number = False - self.model = model - self.settings = model.settings - self.geometry = model.geometry - - # determine set of materials in the model - if not model.materials: - model.materials = openmc.Materials( - model.geometry.get_all_materials().values() - ) - self.materials = model.materials - - self.cleanup_when_done = True - - # Reduce the chain before we create more materials - if reduce_chain: - all_isotopes = set() - for material in self.materials: - if not material.depletable: - continue - for name, _dens_percent, _dens_type in material.nuclides: - all_isotopes.add(name) - self.chain = self.chain.reduce(all_isotopes, reduce_chain_level) - - # Differentiate burnable materials with multiple instances - if diff_burnable_mats: - self._differentiate_burnable_mats() - self.materials = openmc.Materials( - model.geometry.get_all_materials().values() - ) - - # Clear out OpenMC, create task lists, distribute - openmc.reset_auto_ids() - self.burnable_mats, volume, nuclides = self._get_burnable_mats() - self.local_mats = _distribute(self.burnable_mats) - - # Generate map from local materials => material index - self._mat_index_map = { - lm: self.burnable_mats.index(lm) for lm in self.local_mats} - - if self.prev_res is not None: - # Reload volumes into geometry - prev_results[-1].transfer_volumes(self.model) - - # Store previous results in operator - # Distribute reaction rates according to those tracked - # on this process - if comm.size == 1: - self.prev_res = prev_results - else: - self.prev_res = Results() - mat_indexes = _distribute(range(len(self.burnable_mats))) - for res_obj in prev_results: - new_res = res_obj.distribute(self.local_mats, mat_indexes) - self.prev_res.append(new_res) - - # Determine which nuclides have incident neutron data - self.nuclides_with_data = self._get_nuclides_with_data(cross_sections) - - # Select nuclides with data that are also in the chain - self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides - if nuc.name in self.nuclides_with_data] - - # Extract number densities from the geometry / previous depletion run - self._extract_number(self.local_mats, volume, nuclides, self.prev_res) - - # Create reaction rates array - self.reaction_rates = ReactionRates( - self.local_mats, self._burnable_nucs, self.chain.reactions) - - # Get classes to assist working with tallies - if reaction_rate_mode == "direct": - self._rate_helper = DirectReactionRateHelper( - self.reaction_rates.n_nuc, self.reaction_rates.n_react) - elif reaction_rate_mode == "flux": - if reaction_rate_opts is None: - reaction_rate_opts = {} - - # Ensure energy group boundaries were specified - if 'energies' not in reaction_rate_opts: - raise ValueError( - "Energy group boundaries must be specified in the " - "reaction_rate_opts argument when reaction_rate_mode is" - "set to 'flux'.") - - self._rate_helper = FluxCollapseHelper( - self.reaction_rates.n_nuc, - self.reaction_rates.n_react, - **reaction_rate_opts - ) - else: - raise ValueError("Invalid reaction rate mode.") - - if normalization_mode == "fission-q": - self._normalization_helper = ChainFissionHelper() - elif normalization_mode == "energy-deposition": - score = "heating" if self.settings.photon_transport else "heating-local" - self._normalization_helper = EnergyScoreHelper(score) - else: - self._normalization_helper = SourceRateHelper() - - # Select and create fission yield helper - fission_helper = self._fission_helpers[fission_yield_mode] - fission_yield_opts = ( - {} if fission_yield_opts is None else fission_yield_opts) - self._yield_helper = fission_helper.from_operator( - self, **fission_yield_opts) - - def __call__(self, vec, source_rate): - """Runs a simulation. - - Simulation will abort under the following circumstances: - - 1) No energy is computed using OpenMC tallies. - - Parameters - ---------- - vec : list of numpy.ndarray - Total atoms to be used in function. - source_rate : float - Power in [W] or source rate in [neutron/sec] - - Returns - ------- - openmc.deplete.OperatorResult - Eigenvalue and reaction rates resulting from transport operator - - """ - # Reset results in OpenMC - openmc.lib.reset() - - # Update the number densities regardless of the source rate - self.number.set_density(vec) - self._update_materials() - - # If the source rate is zero, return zero reaction rates without running - # a transport solve - if source_rate == 0.0: - rates = self.reaction_rates.copy() - rates.fill(0.0) - return OperatorResult(ufloat(0.0, 0.0), rates) - - # Prevent OpenMC from complaining about re-creating tallies - openmc.reset_auto_ids() - - # Update tally nuclides data in preparation for transport solve - nuclides = self._get_tally_nuclides() - self._rate_helper.nuclides = nuclides - self._normalization_helper.nuclides = nuclides - self._yield_helper.update_tally_nuclides(nuclides) - - # Run OpenMC - openmc.lib.run() - openmc.lib.reset_timers() - - # Extract results - op_result = self._unpack_tallies_and_normalize(source_rate) - - return copy.deepcopy(op_result) - - @staticmethod - def write_bos_data(step): - """Write a state-point file with beginning of step data - - Parameters - ---------- - step : int - Current depletion step including restarts - """ - openmc.lib.statepoint_write( - "openmc_simulation_n{}.h5".format(step), - write_source=False) - - def _differentiate_burnable_mats(self): - """Assign distribmats for each burnable material - - """ - - # Count the number of instances for each cell and material - self.geometry.determine_paths(instances_only=True) - - # Extract all burnable materials which have multiple instances - distribmats = set( - [mat for mat in self.materials - if mat.depletable and mat.num_instances > 1]) - - for mat in distribmats: - if mat.volume is None: - raise RuntimeError("Volume not specified for depletable " - "material with ID={}.".format(mat.id)) - mat.volume /= mat.num_instances - - if distribmats: - # Assign distribmats to cells - for cell in self.geometry.get_all_material_cells().values(): - if cell.fill in distribmats: - mat = cell.fill - cell.fill = [mat.clone() - for i in range(cell.num_instances)] - - def _get_burnable_mats(self): - """Determine depletable materials, volumes, and nuclides - - Returns - ------- - burnable_mats : list of str - List of burnable material IDs - volume : OrderedDict of str to float - Volume of each material in [cm^3] - nuclides : list of str - Nuclides in order of how they'll appear in the simulation. - - """ - - burnable_mats = set() - model_nuclides = set() - volume = OrderedDict() - - self.heavy_metal = 0.0 - - # Iterate once through the geometry to get dictionaries - for mat in self.materials: - for nuclide in mat.get_nuclides(): - model_nuclides.add(nuclide) - if mat.depletable: - burnable_mats.add(str(mat.id)) - if mat.volume is None: - raise RuntimeError("Volume not specified for depletable " - "material with ID={}.".format(mat.id)) - volume[str(mat.id)] = mat.volume - self.heavy_metal += mat.fissionable_mass - - # Make sure there are burnable materials - if not burnable_mats: - raise RuntimeError( - "No depletable materials were found in the model.") - - # Sort the sets - burnable_mats = sorted(burnable_mats, key=int) - model_nuclides = sorted(model_nuclides) - - # Construct a global nuclide dictionary, burned first - nuclides = list(self.chain.nuclide_dict) - for nuc in model_nuclides: - if nuc not in nuclides: - nuclides.append(nuc) - - return burnable_mats, volume, nuclides - - def _extract_number(self, local_mats, volume, nuclides, prev_res=None): - """Construct AtomNumber using geometry - - Parameters - ---------- - local_mats : list of str - Material IDs to be managed by this process - volume : OrderedDict of str to float - Volumes for the above materials in [cm^3] - nuclides : list of str - Nuclides to be used in the simulation. - prev_res : Results, optional - Results from a previous depletion calculation - - """ - self.number = AtomNumber(local_mats, nuclides, volume, len(self.chain)) - - if self.dilute_initial != 0.0: - for nuc in self._burnable_nucs: - self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial) - - # Now extract and store the number densities - # From the geometry if no previous depletion results - if prev_res is None: - for mat in self.materials: - if str(mat.id) in local_mats: - self._set_number_from_mat(mat) - - # Else from previous depletion results - else: - for mat in self.materials: - if str(mat.id) in local_mats: - self._set_number_from_results(mat, prev_res) - - def _set_number_from_mat(self, mat): - """Extracts material and number densities from openmc.Material - - Parameters - ---------- - mat : openmc.Material - The material to read from - - """ - mat_id = str(mat.id) - - for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items(): - atom_per_cc = atom_per_bcm * 1.0e24 - self.number.set_atom_density(mat_id, nuclide, atom_per_cc) - - def _set_number_from_results(self, mat, prev_res): - """Extracts material nuclides and number densities. - - If the nuclide concentration's evolution is tracked, the densities come - from depletion results. Else, densities are extracted from the geometry - in the summary. - - Parameters - ---------- - mat : openmc.Material - The material to read from - prev_res : Results - Results from a previous depletion calculation - - """ - mat_id = str(mat.id) - - # Get nuclide lists from geometry and depletion results - depl_nuc = prev_res[-1].nuc_to_ind - geom_nuc_densities = mat.get_nuclide_atom_densities() - - # Merge lists of nuclides, with the same order for every calculation - geom_nuc_densities.update(depl_nuc) - - for nuclide, atom_per_bcm in geom_nuc_densities.items(): - if nuclide in depl_nuc: - concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1] - volume = prev_res[-1].volume[mat_id] - atom_per_cc = concentration / volume - else: - atom_per_cc = atom_per_bcm * 1.0e24 - - self.number.set_atom_density(mat_id, nuclide, atom_per_cc) - - def initial_condition(self): - """Performs final setup and returns initial condition. - - Returns - ------- - list of numpy.ndarray - Total density for initial conditions. - """ - - # Create XML files - if comm.rank == 0: - self.geometry.export_to_xml() - self.settings.export_to_xml() - self._generate_materials_xml() - - # Initialize OpenMC library - comm.barrier() - if not openmc.lib.is_initialized: - openmc.lib.init(intracomm=comm) - - # Generate tallies in memory - materials = [openmc.lib.materials[int(i)] - for i in self.burnable_mats] - self._rate_helper.generate_tallies(materials, self.chain.reactions) - self._normalization_helper.prepare( - self.chain.nuclides, self.reaction_rates.index_nuc) - # Tell fission yield helper what materials this process is - # responsible for - self._yield_helper.generate_tallies( - materials, tuple(sorted(self._mat_index_map.values()))) - - # Return number density vector - return list(self.number.get_mat_slice(np.s_[:])) - - def finalize(self): - """Finalize a depletion simulation and release resources.""" - if self.cleanup_when_done: - openmc.lib.finalize() - - def _update_materials(self): - """Updates material compositions in OpenMC on all processes.""" - - for rank in range(comm.size): - number_i = comm.bcast(self.number, root=rank) - - for mat in number_i.materials: - nuclides = [] - densities = [] - for nuc in number_i.nuclides: - if nuc in self.nuclides_with_data: - val = 1.0e-24 * number_i.get_atom_density(mat, nuc) - - # If nuclide is zero, do not add to the problem. - if val > 0.0: - if self.round_number: - val_magnitude = np.floor(np.log10(val)) - val_scaled = val / 10**val_magnitude - val_round = round(val_scaled, 8) - - val = val_round * 10**val_magnitude - - nuclides.append(nuc) - densities.append(val) - else: - # Only output warnings if values are significantly - # negative. CRAM does not guarantee positive values. - if val < -1.0e-21: - print("WARNING: nuclide ", nuc, " in material ", mat, - " is negative (density = ", val, " at/barn-cm)") - number_i[mat, nuc] = 0.0 - - # Update densities on C API side - mat_internal = openmc.lib.materials[int(mat)] - mat_internal.set_densities(nuclides, densities) - - #TODO Update densities on the Python side, otherwise the - # summary.h5 file contains densities at the first time step - - def _generate_materials_xml(self): - """Creates materials.xml from self.number. - - Due to uncertainty with how MPI interacts with OpenMC API, this - constructs the XML manually. The long term goal is to do this - through direct memory writing. - - """ - # Sort nuclides according to order in AtomNumber object - nuclides = list(self.number.nuclides) - for mat in self.materials: - mat._nuclides.sort(key=lambda x: nuclides.index(x[0])) - - self.materials.export_to_xml() - - def _get_tally_nuclides(self): - """Determine nuclides that should be tallied for reaction rates. - - This method returns a list of all nuclides that have neutron data and - are listed in the depletion chain. Technically, we should tally nuclides - that may not appear in the depletion chain because we still need to get - the fission reaction rate for these nuclides in order to normalize - power, but that is left as a future exercise. - - Returns - ------- - list of str - Tally nuclides - - """ - nuc_set = set() - - # Create the set of all nuclides in the decay chain in materials marked - # for burning in which the number density is greater than zero. - for nuc in self.number.nuclides: - if nuc in self.nuclides_with_data: - if np.sum(self.number[:, nuc]) > 0.0: - nuc_set.add(nuc) - - # Communicate which nuclides have nonzeros to rank 0 - if comm.rank == 0: - for i in range(1, comm.size): - nuc_newset = comm.recv(source=i, tag=i) - nuc_set |= nuc_newset - else: - comm.send(nuc_set, dest=0, tag=comm.rank) - - if comm.rank == 0: - # Sort nuclides in the same order as self.number - nuc_list = [nuc for nuc in self.number.nuclides - if nuc in nuc_set] - else: - nuc_list = None - - # Store list of tally nuclides on each process - nuc_list = comm.bcast(nuc_list) - return [nuc for nuc in nuc_list if nuc in self.chain] - - def _unpack_tallies_and_normalize(self, source_rate): - """Unpack tallies from OpenMC and return an operator result - - This method uses OpenMC's C API bindings to determine the k-effective - value and reaction rates from the simulation. The reaction rates are - normalized by a helper class depending on the method being used. - - Parameters - ---------- - source_rate : float - Power in [W] or source rate in [neutron/sec] - - Returns - ------- - openmc.deplete.OperatorResult - Eigenvalue and reaction rates resulting from transport operator - - """ - rates = self.reaction_rates - rates.fill(0.0) - - # Get k and uncertainty - keff = ufloat(*openmc.lib.keff()) - - # Extract tally bins - nuclides = self._rate_helper.nuclides - - # Form fast map - nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides] - react_ind = [rates.index_rx[react] for react in self.chain.reactions] - - # Keep track of energy produced from all reactions in eV per source - # particle - self._normalization_helper.reset() - self._yield_helper.unpack() - - # Store fission yield dictionaries - fission_yields = [] - - # Create arrays to store fission Q values, reaction rates, and nuclide - # numbers, zeroed out in material iteration - number = np.empty(rates.n_nuc) - - fission_ind = rates.index_rx.get("fission") - - # Extract results - for i, mat in enumerate(self.local_mats): - # Get tally index - mat_index = self._mat_index_map[mat] - - # Zero out reaction rates and nuclide numbers - number.fill(0.0) - - # Get new number densities - for nuc, i_nuc_results in zip(nuclides, nuc_ind): - number[i_nuc_results] = self.number[mat, nuc] - - tally_rates = self._rate_helper.get_material_rates( - mat_index, nuc_ind, react_ind) - - # Compute fission yields for this material - fission_yields.append(self._yield_helper.weighted_yields(i)) - - # Accumulate energy from fission - if fission_ind is not None: - self._normalization_helper.update(tally_rates[:, fission_ind]) - - # Divide by total number and store - rates[i] = self._rate_helper.divide_by_adens(number) - - # Scale reaction rates to obtain units of reactions/sec - rates *= self._normalization_helper.factor(source_rate) - - # Store new fission yields on the chain - self.chain.fission_yields = fission_yields - - return OperatorResult(keff, rates) - - def _get_nuclides_with_data(self, cross_sections): - """Loads cross_sections.xml file to find nuclides with neutron data""" - nuclides = set() - data_lib = DataLibrary.from_xml(cross_sections) - for library in data_lib.libraries: - if library['type'] != 'neutron': - continue - for name in library['materials']: - if name not in nuclides: - nuclides.add(name) - - return nuclides - - def get_results_info(self): - """Returns volume list, material lists, and nuc lists. - - Returns - ------- - volume : dict of str float - Volumes corresponding to materials in full_burn_dict - nuc_list : list of str - A list of all nuclide names. Used for sorting the simulation. - burn_list : list of int - A list of all material IDs to be burned. Used for sorting the simulation. - full_burn_list : list - List of all burnable material IDs - - """ - nuc_list = self.number.burnable_nuclides - burn_list = self.local_mats - - volume = {} - for i, mat in enumerate(burn_list): - volume[mat] = self.number.volume[i] - - # Combine volume dictionaries across processes - volume_list = comm.allgather(volume) - volume = {k: v for d in volume_list for k, v in d.items()} - - return volume, nuc_list, burn_list, self.burnable_mats diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index c238002fe..b89bae088 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -16,6 +16,19 @@ __all__ = ["Results", "ResultsList"] def _get_time_as(seconds, units): + """Converts the time in seconds to time in different units + + Parameters + ---------- + seconds : float + The time to convert expressed in seconds + units : {"s", "min", "h", "d", "a"} + The units to convert time into. Available options are seconds ``"s"``, + minutes ``"min"``, hours ``"h"`` days ``"d"``, Julian years ``"a"`` + + """ + if units == "a": + return seconds / (60 * 60 * 24 * 365.25) # 365.25 due to the leap year if units == "d": return seconds / (60 * 60 * 24) elif units == "h": @@ -81,9 +94,9 @@ class Results(list): .. note:: Initial values for some isotopes that do not appear in initial concentrations may be non-zero, depending on the - value of the :attr:`openmc.deplete.Operator.dilute_initial` - attribute. The :class:`openmc.deplete.Operator` class adds isotopes - according to this setting, which can be set to zero. + value of the :attr:`openmc.deplete.CoupledOperator.dilute_initial` + attribute. The :class:`openmc.deplete.CoupledOperator` class adds + isotopes according to this setting, which can be set to zero. Parameters ---------- @@ -95,9 +108,10 @@ class Results(list): Units for the returned concentration. Default is ``"atoms"`` .. versionadded:: 0.12 - time_units : {"s", "min", "h", "d"}, optional + time_units : {"s", "min", "h", "d", "a"}, optional Units for the returned time array. Default is ``"s"`` to - return the value in seconds. + return the value in seconds. Other options are minutes ``"min"``, + hours ``"h"``, days ``"d"``, and Julian years ``"a"``. .. versionadded:: 0.12 @@ -109,7 +123,7 @@ class Results(list): Concentration of specified nuclide in units of ``nuc_units`` """ - cv.check_value("time_units", time_units, {"s", "d", "min", "h"}) + cv.check_value("time_units", time_units, {"s", "d", "min", "h", "a"}) cv.check_value("nuc_units", nuc_units, {"atoms", "atom/b-cm", "atom/cm3"}) @@ -145,8 +159,8 @@ class Results(list): Initial values for some isotopes that do not appear in initial concentrations may be non-zero, depending on the - value of :class:`openmc.deplete.Operator` ``dilute_initial`` - The :class:`openmc.deplete.Operator` adds isotopes according + value of :class:`openmc.deplete.CoupledOperator` ``dilute_initial`` + The :class:`openmc.deplete.CoupledOperator` adds isotopes according to this setting, which can be set to zero. Parameters @@ -190,8 +204,10 @@ class Results(list): Parameters ---------- - time_units : {"s", "d", "h", "min"}, optional - Desired units for the times array + time_units : {"s", "d", "min", "h", "a"}, optional + Desired units for the times array. Options are seconds ``"s"``, + minutes ``"min"``, hours ``"h"``, days ``"d"``, and Julian years + ``"a"``. Returns ------- @@ -203,7 +219,7 @@ class Results(list): 1 contains the associated uncertainty """ - cv.check_value("time_units", time_units, {"s", "d", "min", "h"}) + cv.check_value("time_units", time_units, {"s", "d", "min", "h", "a"}) times = np.empty_like(self, dtype=float) eigenvalues = np.empty((len(self), 2), dtype=float) @@ -257,9 +273,10 @@ class Results(list): Parameters ---------- - time_units : {"s", "d", "h", "min"}, optional + time_units : {"s", "d", "min", "h", "a"}, optional Return the vector in these units. Default is to - convert to days + convert to days ``"d"``. Other options are seconds ``"s"``, minutes + ``"min"``, hours ``"h"``, days ``"d"``, and Julian years ``"a"``. Returns ------- @@ -267,7 +284,7 @@ class Results(list): 1-D vector of time points """ - cv.check_value("time_units", time_units, {"s", "d", "min", "h"}) + cv.check_value("time_units", time_units, {"s", "d", "min", "h", "a"}) times = np.fromiter( (r.time[0] for r in self), @@ -296,8 +313,9 @@ class Results(list): ---------- time : float Desired point in time - time_units : {"s", "d", "min", "h"}, optional - Units on ``time``. Default: days + time_units : {"s", "d", "min", "h", "a"}, optional + Units on ``time``. Default: days ``"d"``. Other options are seconds + ``"s"``, minutes ``"min"``, hours ``"h"`` and Julian years ``"a"``. atol : float, optional Absolute tolerance (in ``time_units``) if ``time`` is not found. diff --git a/openmc/deplete/stepresult.py b/openmc/deplete/stepresult.py index 6e2373c45..c54f9edcf 100644 --- a/openmc/deplete/stepresult.py +++ b/openmc/deplete/stepresult.py @@ -462,7 +462,7 @@ class StepResult: Parameters ---------- - op : openmc.deplete.TransportOperator + op : openmc.deplete.abc.TransportOperator The operator used to generate these results. x : list of list of numpy.array The prior x vectors. Indexed [i][cell] using the above equation. @@ -495,7 +495,13 @@ class StepResult: for mat_i in range(n_mat): results[i, mat_i, :] = x[i][mat_i] - results.k = [(r.k.nominal_value, r.k.std_dev) for r in op_results] + ks = [] + for r in op_results: + if isinstance(r.k, type(None)): + ks += [(None, None)] + else: + ks += [(r.k.nominal_value, r.k.std_dev)] + results.k = ks results.rates = [r.rates for r in op_results] results.time = t results.source_rate = source_rate diff --git a/openmc/material.py b/openmc/material.py index 6c2641680..913d53d9d 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -861,7 +861,7 @@ class Material(IDManagerMixin): elif self.density_units == 'atom/b-cm': density = self.density elif self.density_units == 'atom/cm3' or self.density_units == 'atom/cc': - density = 1.E-24 * self.density + density = 1.e-24 * self.density # For ease of processing split out nuc, nuc_density, # and nuc_density_type into separate arrays @@ -897,7 +897,7 @@ class Material(IDManagerMixin): # Convert the mass density to an atom density if not density_in_atom: - density = -density / self.average_molar_mass * 1.E-24 \ + density = -density / self.average_molar_mass * 1.e-24 \ * openmc.data.AVOGADRO nuc_densities = density * nuc_densities diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 6271154fb..51019b0d8 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -451,10 +451,10 @@ class RectangularParallelepiped(CompositeSurface): self.zmax = openmc.ZPlane(z0=zmax, **kwargs) def __neg__(self): - return +self.xmin & -self.xmax & +self.ymin & -self.ymax & +self.zmin & -self.zmax + return -self.xmax & +self.xmin & -self.ymax & +self.ymin & -self.zmax & +self.zmin def __pos__(self): - return -self.xmin | +self.xmax | -self.ymin | +self.ymax | -self.zmin | +self.zmax + return +self.xmax | -self.xmin | +self.ymax | -self.ymin | +self.zmax | -self.zmin class XConeOneSided(CompositeSurface): diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 5528bdb30..58b7c6172 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -1,6 +1,7 @@ from abc import ABC, abstractmethod from collections import defaultdict from collections.abc import Iterable +from copy import deepcopy from numbers import Real from xml.etree import ElementTree as ET @@ -95,9 +96,9 @@ class Discrete(Univariate): Attributes ---------- - x : Iterable of float + x : numpy.ndarray Values of the random variable - p : Iterable of float + p : numpy.ndarray Discrete probability for each value """ @@ -122,7 +123,7 @@ class Discrete(Univariate): if isinstance(x, Real): x = [x] cv.check_type('discrete values', x, Iterable, Real) - self._x = x + self._x = np.array(x, dtype=float) @p.setter def p(self, p): @@ -131,14 +132,15 @@ class Discrete(Univariate): cv.check_type('discrete probabilities', p, Iterable, Real) for pk in p: cv.check_greater_than('discrete probability', pk, 0.0, True) - self._p = p + self._p = np.array(p, dtype=float) def cdf(self): return np.insert(np.cumsum(self.p), 0, 0.0) def sample(self, n_samples=1, seed=None): np.random.seed(seed) - return np.random.choice(self.x, n_samples, p=self.p) + p = self.p / self.p.sum() + return np.random.choice(self.x, n_samples, p=p) def normalize(self): """Normalize the probabilities stored on the distribution""" @@ -220,6 +222,18 @@ class Discrete(Univariate): p_arr = np.array([p_merged[x] for x in x_arr]) return cls(x_arr, p_arr) + def integral(self): + """Return integral of distribution + + .. versionadded:: 0.13.1 + + Returns + ------- + float + Integral of discrete distribution + """ + return np.sum(self.p) + class Uniform(Univariate): """Distribution with constant probability over a finite interval [a,b] @@ -825,9 +839,9 @@ class Tabular(Univariate): Attributes ---------- - x : Iterable of float + x : numpy.ndarray Tabulated values of the random variable - p : Iterable of float + p : numpy.ndarray Tabulated probabilities interpolation : {'histogram', 'linear-linear', 'linear-log', 'log-linear', 'log-log'}, optional Indicate whether the density function is constant between tabulated @@ -860,7 +874,7 @@ class Tabular(Univariate): @x.setter def x(self, x): cv.check_type('tabulated values', x, Iterable, Real) - self._x = x + self._x = np.array(x, dtype=float) @p.setter def p(self, p): @@ -868,7 +882,7 @@ class Tabular(Univariate): if not self._ignore_negative: for pk in p: cv.check_greater_than('tabulated probability', pk, 0.0, True) - self._p = p + self._p = np.array(p, dtype=float) @interpolation.setter def interpolation(self, interpolation): @@ -881,8 +895,8 @@ class Tabular(Univariate): 'distributions using histogram or ' 'linear-linear interpolation') c = np.zeros_like(self.x) - x = np.asarray(self.x) - p = np.asarray(self.p) + x = self.x + p = self.p if self.interpolation == 'histogram': c[1:] = p[:-1] * np.diff(x) @@ -922,7 +936,7 @@ class Tabular(Univariate): def normalize(self): """Normalize the probabilities stored on the distribution""" - self.p = np.asarray(self.p) / self.cdf().max() + self.p /= self.cdf().max() def sample(self, n_samples=1, seed=None): if not self.interpolation in ('histogram', 'linear-linear'): @@ -931,10 +945,11 @@ class Tabular(Univariate): 'linear-linear interpolation') np.random.seed(seed) xi = np.random.rand(n_samples) - cdf = self.cdf() - cdf /= cdf.max() + # always use normalized probabilities when sampling + cdf = self.cdf() p = self.p / cdf.max() + cdf /= cdf.max() # get CDF bins that are above the # sampled values @@ -949,7 +964,7 @@ class Tabular(Univariate): # the random number is less than the next cdf # entry x_i = self.x[cdf_idx] - p_i = self.p[cdf_idx] + p_i = p[cdf_idx] if self.interpolation == 'histogram': # mask where probability is greater than zero @@ -967,7 +982,7 @@ class Tabular(Univariate): # get variable and probability values for the # next entry x_i1 = self.x[cdf_idx + 1] - p_i1 = self.p[cdf_idx + 1] + p_i1 = p[cdf_idx + 1] # compute slope between entries m = (p_i1 - p_i) / (x_i1 - x_i) # set values for zero slope @@ -1027,6 +1042,24 @@ class Tabular(Univariate): p = params[len(params)//2:] return cls(x, p, interpolation) + def integral(self): + """Return integral of distribution + + .. versionadded: 0.13.1 + + Returns + ------- + float + Integral of tabular distrbution + """ + if self.interpolation == 'histogram': + return np.sum(np.diff(self.x) * self.p[:-1]) + elif self.interpolation == 'linear-linear': + return np.trapz(self.p, self.x) + else: + raise NotImplementedError( + f'integral() not supported for {self.inteprolation} interpolation') + class Legendre(Univariate): r"""Probability density given by a Legendre polynomial expansion @@ -1135,9 +1168,10 @@ class Mixture(Univariate): def sample(self, n_samples=1, seed=None): np.random.seed(seed) - idx = np.random.choice(self.distribution, n_samples, p=self.probability) + idx = np.random.choice(range(len(self.distribution)), + n_samples, p=self.probability) - out = np.zeros_like(idx) + out = np.empty_like(idx, dtype=float) for i in np.unique(idx): n_dist_samples = np.count_nonzero(idx == i) samples = self.distribution[i].sample(n_dist_samples) @@ -1199,3 +1233,68 @@ class Mixture(Univariate): distribution.append(Univariate.from_xml_element(pair.find("dist"))) return cls(probability, distribution) + + def integral(self): + """Return integral of the distribution + + .. versionadded:: 0.13.1 + + Returns + ------- + float + Integral of the distribution + """ + return sum([ + p*dist.integral() + for p, dist in zip(self.probability, self.distribution) + ]) + + +def combine_distributions(dists, probs): + """Combine distributions with specified probabilities + + This function can be used to combine multiple instances of + :class:`~openmc.stats.Discrete` and `~openmc.stats.Tabular`. Multiple + discrete distributions are merged into a single distribution and the + remainder of the distributions are put into a :class:`~openmc.stats.Mixture` + distribution. + + .. versionadded:: 0.13.1 + + Parameters + ---------- + dists : iterable of openmc.stats.Univariate + Distributions to combine + probs : iterable of float + Probability (or intensity) of each distribution + + """ + # Get copy of distribution list so as not to modify the argument + dist_list = deepcopy(dists) + + # Get list of discrete/continuous distribution indices + discrete_index = [i for i, d in enumerate(dist_list) if isinstance(d, Discrete)] + cont_index = [i for i, d in enumerate(dist_list) if isinstance(d, Tabular)] + + # Apply probabilites to continuous distributions + for i in cont_index: + dist = dist_list[i] + dist.p *= probs[i] + + if discrete_index: + # Create combined discrete distribution + dist_discrete = [dist_list[i] for i in discrete_index] + discrete_probs = [probs[i] for i in discrete_index] + combined_dist = Discrete.merge(dist_discrete, discrete_probs) + + # Replace multiple discrete distributions with merged + for idx in reversed(discrete_index): + dist_list.pop(idx) + dist_list.append(combined_dist) + + # Combine discrete and continuous if present + if len(dist_list) > 1: + probs = [d.integral() for d in dist_list] + dist_list[:] = [Mixture(probs, dist_list.copy())] + + return dist_list[0] diff --git a/openmc/universe.py b/openmc/universe.py index 78f75bbfc..3ef9e6650 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -334,6 +334,13 @@ class Universe(UniverseBase): if seed is not None: model.settings.seed = seed + # Determine whether any materials contains macroscopic data and if + # so, set energy mode accordingly + for mat in self.get_all_materials().values(): + if mat._macroscopic is not None: + model.settings.energy_mode = 'multi-group' + break + # Create plot object matching passed arguments plot = openmc.Plot() plot.origin = origin @@ -785,7 +792,7 @@ class DAGMCUniverse(UniverseBase): Universe instance """ - bounding_cell = openmc.Cell(fill=self, cell_id =bounding_cell_id, region=self.bounding_region(**kwargs)) + bounding_cell = openmc.Cell(fill=self, cell_id=bounding_cell_id, region=self.bounding_region(**kwargs)) return openmc.Universe(cells=[bounding_cell]) @classmethod diff --git a/tests/micro_xs_simple.csv b/tests/micro_xs_simple.csv new file mode 100644 index 000000000..787ce74c3 --- /dev/null +++ b/tests/micro_xs_simple.csv @@ -0,0 +1,13 @@ +nuclide,"(n,gamma)",fission +U234,22.231989822002465,0.49620744663749855 +U235,10.479008971197121,48.41787337164604 +U238,0.8673334105437324,0.10467880588762352 +U236,8.65171044607122,0.31948392400019293 +O16,7.497851000107524e-05,0.0 +O17,0.0004079227797153371,0.0 +I135,6.842395323713927,0.0 +Xe135,227463.86426990604,0.0 +Xe136,0.02317896034753588,0.0 +Cs135,2.1721665580713623,0.0 +Gd157,12786.09939237018,0.0 +Gd156,3.4006085445846983,0.0 diff --git a/tests/regression_tests/cpp_driver/test.py b/tests/regression_tests/cpp_driver/test.py index 95912b3f7..e8b1c62ac 100644 --- a/tests/regression_tests/cpp_driver/test.py +++ b/tests/regression_tests/cpp_driver/test.py @@ -48,8 +48,8 @@ def cpp_driver(request): finally: # Remove local build directory when test is complete - shutil.rmtree('build') - os.remove('CMakeLists.txt') + shutil.rmtree(request.node.path.parent / 'build') + os.remove(request.node.path.parent / 'CMakeLists.txt') @pytest.fixture diff --git a/tests/regression_tests/deplete_no_transport/__init__.py b/tests/regression_tests/deplete_no_transport/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/deplete_no_transport/test.py b/tests/regression_tests/deplete_no_transport/test.py new file mode 100644 index 000000000..95b842537 --- /dev/null +++ b/tests/regression_tests/deplete_no_transport/test.py @@ -0,0 +1,106 @@ +""" Transport-free depletion test suite """ + +from pathlib import Path + +import numpy as np +import pytest +import openmc +import openmc.deplete +from openmc.deplete import IndependentOperator, MicroXS + +@pytest.fixture(scope="module") +def fuel(): + fuel = openmc.Material(name="uo2") + fuel.add_element("U", 1, percent_type="ao", enrichment=4.25) + fuel.add_element("O", 2) + fuel.set_density("g/cc", 10.4) + fuel.depletable = True + + fuel.volume = np.pi * 0.42 ** 2 + + return fuel + +@pytest.mark.parametrize("multiproc, from_nuclides, normalization_mode, power, flux", [ + (True, True,'source-rate', None, 1164719970082145.0), + (False, True, 'source-rate', None, 1164719970082145.0), + (True, True, 'fission-q', 174, None), + (False, True, 'fission-q', 174, None), + (True, False,'source-rate', None, 1164719970082145.0), + (False, False, 'source-rate', None, 1164719970082145.0), + (True, False, 'fission-q', 174, None), + (False, False, 'fission-q', 174, None)]) +def test_no_transport_from_nuclides(run_in_tmpdir, fuel, multiproc, from_nuclides, normalization_mode, power, flux): + """Transport free system test suite. + + Runs an OpenMC transport-free depletion calculation and verifies + that the outputs match a reference file. + + """ + # Create operator + micro_xs_file = Path(__file__).parents[2] / 'micro_xs_simple.csv' + micro_xs = MicroXS.from_csv(micro_xs_file) + chain_file = Path(__file__).parents[2] / 'chain_simple.xml' + + if from_nuclides: + nuclides = {} + for nuc, dens in fuel.get_nuclide_atom_densities().items(): + nuclides[nuc] = dens + + op = IndependentOperator.from_nuclides( + fuel.volume, nuclides, micro_xs, chain_file, normalization_mode=normalization_mode) + + else: + op = IndependentOperator(openmc.Materials([fuel]), micro_xs, chain_file, normalization_mode=normalization_mode) + + # Power and timesteps + dt = [30] # single step + + # Perform simulation using the predictor algorithm + openmc.deplete.pool.USE_MULTIPROCESSING = multiproc + openmc.deplete.PredictorIntegrator( + op, dt, power=power, source_rates=flux, timestep_units='d').integrate() + + # Get path to test and reference results + path_test = op.output_dir / 'depletion_results.h5' + if flux is not None: + ref_path = 'test_reference_source_rate.h5' + else: + ref_path = 'test_reference_fission_q.h5' + path_reference = Path(__file__).with_name(ref_path) + + # Load the reference/test results + res_test = openmc.deplete.Results(path_test) + res_ref = openmc.deplete.Results(path_reference) + + # Assert same mats + for mat in res_ref[0].mat_to_ind: + assert mat in res_test[0].mat_to_ind, \ + "Material {} not in new results.".format(mat) + for nuc in res_ref[0].nuc_to_ind: + assert nuc in res_test[0].nuc_to_ind, \ + "Nuclide {} not in new results.".format(nuc) + + for mat in res_test[0].mat_to_ind: + assert mat in res_ref[0].mat_to_ind, \ + "Material {} not in old results.".format(mat) + for nuc in res_test[0].nuc_to_ind: + assert nuc in res_ref[0].nuc_to_ind, \ + "Nuclide {} not in old results.".format(nuc) + + tol = 1.0e-6 + for mat in res_test[0].mat_to_ind: + for nuc in res_test[0].nuc_to_ind: + _, y_test = res_test.get_atoms(mat, nuc) + _, y_old = res_ref.get_atoms(mat, nuc) + + # Test each point + correct = True + for i, ref in enumerate(y_old): + if ref != y_test[i]: + if ref != 0.0: + correct = np.abs(y_test[i] - ref) / ref <= tol + else: + correct = False + + assert correct, "Discrepancy in mat {} and nuc {}\n{}\n{}".format( + mat, nuc, y_old, y_test) diff --git a/tests/regression_tests/deplete_no_transport/test_reference_fission_q.h5 b/tests/regression_tests/deplete_no_transport/test_reference_fission_q.h5 new file mode 100644 index 000000000..82bf11567 Binary files /dev/null and b/tests/regression_tests/deplete_no_transport/test_reference_fission_q.h5 differ diff --git a/tests/regression_tests/deplete_no_transport/test_reference_source_rate.h5 b/tests/regression_tests/deplete_no_transport/test_reference_source_rate.h5 new file mode 100644 index 000000000..88ed73b17 Binary files /dev/null and b/tests/regression_tests/deplete_no_transport/test_reference_source_rate.h5 differ diff --git a/tests/regression_tests/deplete_with_transport/__init__.py b/tests/regression_tests/deplete_with_transport/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/deplete/example_geometry.py b/tests/regression_tests/deplete_with_transport/example_geometry.py similarity index 100% rename from tests/regression_tests/deplete/example_geometry.py rename to tests/regression_tests/deplete_with_transport/example_geometry.py diff --git a/tests/regression_tests/deplete/last_step_reference_materials.xml b/tests/regression_tests/deplete_with_transport/last_step_reference_materials.xml similarity index 100% rename from tests/regression_tests/deplete/last_step_reference_materials.xml rename to tests/regression_tests/deplete_with_transport/last_step_reference_materials.xml diff --git a/tests/regression_tests/deplete/test.py b/tests/regression_tests/deplete_with_transport/test.py similarity index 99% rename from tests/regression_tests/deplete/test.py rename to tests/regression_tests/deplete_with_transport/test.py index 6c431fb33..0c2a4cd97 100644 --- a/tests/regression_tests/deplete/test.py +++ b/tests/regression_tests/deplete_with_transport/test.py @@ -12,7 +12,7 @@ from openmc.data import JOULE_PER_EV import openmc.deplete from tests.regression_tests import config -from example_geometry import generate_problem +from .example_geometry import generate_problem @pytest.fixture(scope="module") diff --git a/tests/regression_tests/deplete/test_reference.h5 b/tests/regression_tests/deplete_with_transport/test_reference.h5 similarity index 100% rename from tests/regression_tests/deplete/test_reference.h5 rename to tests/regression_tests/deplete_with_transport/test_reference.h5 diff --git a/tests/regression_tests/microxs/__init.py__ b/tests/regression_tests/microxs/__init.py__ new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/microxs/test.py b/tests/regression_tests/microxs/test.py new file mode 100644 index 000000000..dbfe3357b --- /dev/null +++ b/tests/regression_tests/microxs/test.py @@ -0,0 +1,52 @@ +"""Test one-group cross section generation""" +from pathlib import Path + +import numpy as np +import pytest +import openmc + +from openmc.deplete import MicroXS + +CHAIN_FILE = Path(__file__).parents[2] / "chain_simple.xml" + +@pytest.fixture(scope="module") +def model(): + fuel = openmc.Material(name="uo2") + fuel.add_element("U", 1, percent_type="ao", enrichment=4.25) + fuel.add_element("O", 2) + fuel.set_density("g/cc", 10.4) + + clad = openmc.Material(name="clad") + clad.add_element("Zr", 1) + clad.set_density("g/cc", 6) + + water = openmc.Material(name="water") + water.add_element("O", 1) + water.add_element("H", 2) + water.set_density("g/cc", 1.0) + water.add_s_alpha_beta("c_H_in_H2O") + + radii = [0.42, 0.45] + fuel.volume = np.pi * radii[0] ** 2 + + materials = openmc.Materials([fuel, clad, water]) + + pin_surfaces = [openmc.ZCylinder(r=r) for r in radii] + pin_univ = openmc.model.pin(pin_surfaces, materials) + bound_box = openmc.rectangular_prism(1.24, 1.24, boundary_type="reflective") + root_cell = openmc.Cell(fill=pin_univ, region=bound_box) + geometry = openmc.Geometry([root_cell]) + + settings = openmc.Settings() + settings.particles = 1000 + settings.inactive = 10 + settings.batches = 50 + + return openmc.Model(geometry, materials, settings) + + +def test_from_model(model): + ref_xs = MicroXS.from_csv('test_reference.csv') + test_xs = MicroXS.from_model(model, model.materials[0], CHAIN_FILE) + + np.testing.assert_allclose(test_xs, ref_xs, rtol=1e-11) diff --git a/tests/regression_tests/microxs/test_reference.csv b/tests/regression_tests/microxs/test_reference.csv new file mode 100644 index 000000000..787ce74c3 --- /dev/null +++ b/tests/regression_tests/microxs/test_reference.csv @@ -0,0 +1,13 @@ +nuclide,"(n,gamma)",fission +U234,22.231989822002465,0.49620744663749855 +U235,10.479008971197121,48.41787337164604 +U238,0.8673334105437324,0.10467880588762352 +U236,8.65171044607122,0.31948392400019293 +O16,7.497851000107524e-05,0.0 +O17,0.0004079227797153371,0.0 +I135,6.842395323713927,0.0 +Xe135,227463.86426990604,0.0 +Xe136,0.02317896034753588,0.0 +Cs135,2.1721665580713623,0.0 +Gd157,12786.09939237018,0.0 +Gd156,3.4006085445846983,0.0 diff --git a/tests/regression_tests/torus/inputs_true.dat b/tests/regression_tests/torus/inputs_true.dat index 0045d7e46..0ccbda400 100644 --- a/tests/regression_tests/torus/inputs_true.dat +++ b/tests/regression_tests/torus/inputs_true.dat @@ -3,7 +3,7 @@ - + diff --git a/tests/unit_tests/test_data_decay.py b/tests/unit_tests/test_data_decay.py index 06b8e6bed..f0bda1bd9 100644 --- a/tests/unit_tests/test_data_decay.py +++ b/tests/unit_tests/test_data_decay.py @@ -23,6 +23,14 @@ def nb90(): return openmc.data.Decay.from_endf(filename) +@pytest.fixture(scope='module') +def ba137m(): + """Ba137_m1 decay data.""" + endf_data = os.environ['OPENMC_ENDF_DATA'] + filename = os.path.join(endf_data, 'decay', 'dec-056_Ba_137m1.endf') + return openmc.data.Decay.from_endf(filename) + + @pytest.fixture(scope='module') def u235_yields(): """U235 fission product yield data.""" @@ -48,6 +56,7 @@ def test_nb90_halflife(nb90): ufloat_close(nb90.decay_constant, log(2.)/nb90.half_life) ufloat_close(nb90.decay_energy, ufloat(2265527.5, 25159.400474401213)) + def test_nb90_nuclide(nb90): assert nb90.nuclide['atomic_number'] == 41 assert nb90.nuclide['mass_number'] == 90 @@ -91,3 +100,29 @@ def test_fpy(u235_yields): assert len(u235_yields.independent) == 3 thermal = u235_yields.independent[0] ufloat_close(thermal['I135'], ufloat(0.0292737, 0.000819663)) + + +def test_sources(ba137m, nb90): + # Running .sources twice should give same objects + sources = ba137m.sources + sources2 = ba137m.sources + for key in sources: + assert sources[key] is sources2[key] + + # Each source should be a univariate distribution + for dist in sources.values(): + assert isinstance(dist, openmc.stats.Univariate) + + # Check for presence of 662 keV gamma ray in decay of Ba137m + gamma_source = ba137m.sources['photon'] + assert isinstance(gamma_source, openmc.stats.Discrete) + b = np.isclose(gamma_source.x, 661657.) + assert np.count_nonzero(b) == 1 + + # Check value of decay/s/atom + idx = np.flatnonzero(b)[0] + assert gamma_source.p[idx] == pytest.approx(0.004069614) + + # Nb90 decays by β+ and should emit positrons, electrons, and photons + sources = nb90.sources + assert len(set(sources.keys()) ^ {'positron', 'electron', 'photon'}) == 0 diff --git a/tests/unit_tests/test_data_kalbach_mann.py b/tests/unit_tests/test_data_kalbach_mann.py new file mode 100644 index 000000000..3b837af7d --- /dev/null +++ b/tests/unit_tests/test_data_kalbach_mann.py @@ -0,0 +1,176 @@ +"""Test of the Kalbach-Mann slope calculation when data are +retrieved from ENDF files.""" + +import os +from pathlib import Path +import pytest + +import numpy as np + +from openmc.data import IncidentNeutron +from openmc.data.kalbach_mann import _separation_energy, _AtomicRepresentation +from openmc.data import kalbach_slope +from openmc.data import KalbachMann + +from . import needs_njoy + + +@pytest.fixture(scope='module') +def neutron(): + """Neutron AtomicRepresentation.""" + return _AtomicRepresentation(z=0, a=1) + + +@pytest.fixture(scope='module') +def triton(): + """Triton AtomicRepresentation.""" + return _AtomicRepresentation(z=1, a=3) + + +@pytest.fixture(scope='module') +def b10(): + """B10 AtomicRepresentation.""" + return _AtomicRepresentation(z=5, a=10) + + +@pytest.fixture(scope='module') +def c12(): + """C12 AtomicRepresentation.""" + return _AtomicRepresentation(z=6, a=12) + + +@pytest.fixture(scope='module') +def c13(): + """C13 AtomicRepresentation.""" + return _AtomicRepresentation(z=6, a=13) + + +@pytest.fixture(scope='module') +def na23(): + """Na23 AtomicRepresentation.""" + return _AtomicRepresentation(z=11, a=23) + + +def test_atomic_representation(neutron, triton, b10, c12, c13, na23): + """Test the _AtomicRepresentation class.""" + # Test instantiation from_za + assert b10 == _AtomicRepresentation.from_za(5010) + + # Test addition + assert c13 + b10 == na23 + + # Test substraction + assert c13 - c12 == neutron + assert c13 - b10 == triton + + # Test properties when no information for Kalbach-Mann are given + assert c13.a == 13 + assert c13.z == 6 + assert c13.n == 7 + assert c13.za == 6013 + + # Test properties when information for Kalbach-Mann are given + assert triton.a == 3 + assert triton.z == 1 + assert triton.n == 2 + assert triton.za == 1003 + + # Test instantiation errors + with pytest.raises(ValueError): + _AtomicRepresentation(z=5, a=1) + with pytest.raises(ValueError): + _AtomicRepresentation(z=-1, a=1) + with pytest.raises(ValueError): + _AtomicRepresentation(z=5, a=0) + with pytest.raises(ValueError): + _AtomicRepresentation(z=5, a=-2) + with pytest.raises(ValueError): + neutron - triton + + +def test_separation_energy(triton, b10, c13): + """Comparison to hand-calculations on a simple example.""" + assert _separation_energy( + compound=c13, + nucleus=b10, + particle=triton + ) == pytest.approx(18.6880713) + + +def test_kalbach_slope(): + """Comparison to hand-calculations for n + c12 -> c13 -> triton + b10.""" + energy_projectile = 10.2 # [eV] + energy_emitted = 5.4 # [eV] + + # Check that NotImplementedError is raised if the projectile is not + # a neutron + with pytest.raises(NotImplementedError): + kalbach_slope( + energy_projectile=energy_projectile, + energy_emitted=energy_emitted, + za_projectile=1000, + za_emitted=1, + za_target=6012 + ) + + assert kalbach_slope( + energy_projectile=energy_projectile, + energy_emitted=energy_emitted, + za_projectile=1, + za_emitted=1003, + za_target=6012 + ) == pytest.approx(0.8409921475) + + +@pytest.mark.parametrize( + "hdf5_filename, endf_filename", [ + ('O16.h5', 'n-008_O_016.endf'), + ('Ca46.h5', 'n-020_Ca_046.endf'), + ('Hg204.h5', 'n-080_Hg_204.endf') + ] +) +def test_comparison_slope_hdf5(hdf5_filename, endf_filename): + """Test the calculation of the Kalbach-Mann slope done by OpenMC + by comparing it to HDF5 data. The test is based on the first product + of MT=5 (neutron). The isotopes tested have been selected because the + corresponding products in ENDF/B-VII.1 are described using MF=6, LAW=1, + LANG=2 (i.e., Kalbach-Mann systematics) and the slope is not given + explicitly. + + If an error occurs during the "validity check", this means that + the nuclear data evaluation has evolved and the distribution might + no longer be described using Kalbach-Mann systematics. Another + isotope needs to be identified and tested. + + Warning: This test is valid as long as ENDF files are not directly + used to generate the HDF5 files used in the tests. + + """ + # HDF5 data + hdf5_directory = Path(os.environ['OPENMC_CROSS_SECTIONS']).parent + hdf5_data = IncidentNeutron.from_hdf5(hdf5_directory / hdf5_filename) + hdf5_product = hdf5_data[5].products[0] + hdf5_distribution = hdf5_product.distribution[0] + + # ENDF data + endf_directory = Path(os.environ['OPENMC_ENDF_DATA']) + endf_path = endf_directory / 'neutrons' / endf_filename + endf_data = IncidentNeutron.from_endf(endf_path) + endf_product = endf_data[5].products[0] + endf_distribution = endf_product.distribution[0] + + # Validity check + assert isinstance(endf_distribution, KalbachMann) + assert isinstance(hdf5_distribution, KalbachMann) + assert endf_product.particle == hdf5_product.particle + assert len(endf_distribution.slope) == len(hdf5_distribution.slope) + + # Results check + for i, hdf5_slope in enumerate(hdf5_distribution.slope): + assert endf_distribution._calculated_slope[i] + + np.testing.assert_array_almost_equal( + endf_distribution.slope[i].y, + hdf5_slope.y, + decimal=5 + ) diff --git a/tests/unit_tests/test_deplete_coupled_operator.py b/tests/unit_tests/test_deplete_coupled_operator.py new file mode 100644 index 000000000..6de0342b0 --- /dev/null +++ b/tests/unit_tests/test_deplete_coupled_operator.py @@ -0,0 +1,55 @@ +"""Basic unit tests for openmc.deplete.CoupledOperator instantiation + +""" + +from pathlib import Path + +import pytest +from openmc.deplete import CoupledOperator +import openmc +import numpy as np + +CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml" + +@pytest.fixture(scope="module") +def model(): + fuel = openmc.Material(name="uo2") + fuel.add_element("U", 1, percent_type="ao", enrichment=4.25) + fuel.add_element("O", 2) + fuel.set_density("g/cc", 10.4) + + clad = openmc.Material(name="clad") + clad.add_element("Zr", 1) + clad.set_density("g/cc", 6) + + water = openmc.Material(name="water") + water.add_element("O", 1) + water.add_element("H", 2) + water.set_density("g/cc", 1.0) + water.add_s_alpha_beta("c_H_in_H2O") + + radii = [0.42, 0.45] + fuel.volume = np.pi * radii[0] ** 2 + clad.volume = np.pi * (radii[1]**2 - radii[0]**2) + water.volume = 1.24**2 - (np.pi * radii[1]**2) + + materials = openmc.Materials([fuel, clad, water]) + + pin_surfaces = [openmc.ZCylinder(r=r) for r in radii] + pin_univ = openmc.model.pin(pin_surfaces, materials) + bound_box = openmc.rectangular_prism(1.24, 1.24, boundary_type="reflective") + root_cell = openmc.Cell(fill=pin_univ, region=bound_box) + geometry = openmc.Geometry([root_cell]) + + settings = openmc.Settings() + settings.particles = 1000 + settings.inactive = 10 + settings.batches = 50 + + return openmc.Model(geometry, materials, settings) + +def test_operator_init(model): + """The test uses a temporary dummy chain. This file will be removed + at the end of the test, and only contains a depletion_chain node.""" + + CoupledOperator(model, CHAIN_PATH) diff --git a/tests/unit_tests/test_deplete_independent_operator.py b/tests/unit_tests/test_deplete_independent_operator.py new file mode 100644 index 000000000..f6cc7e200 --- /dev/null +++ b/tests/unit_tests/test_deplete_independent_operator.py @@ -0,0 +1,36 @@ +"""Basic unit tests for openmc.deplete.IndependentOperator instantiation + +""" + +from pathlib import Path + +import pytest +from openmc.deplete import IndependentOperator, MicroXS +from openmc import Material, Materials +import numpy as np + +CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml" +ONE_GROUP_XS = Path(__file__).parents[1] / "micro_xs_simple.csv" + +def test_operator_init(): + """The test uses a temporary dummy chain. This file will be removed + at the end of the test, and only contains a depletion_chain node.""" + volume = 1 + nuclides = {'U234': 8.922411359424315e+18, + 'U235': 9.98240191860822e+20, + 'U238': 2.2192386373095893e+22, + 'U236': 4.5724195495061115e+18, + 'O16': 4.639065406771322e+22, + 'O17': 1.7588724018066158e+19} + micro_xs = MicroXS.from_csv(ONE_GROUP_XS) + IndependentOperator.from_nuclides( + volume, nuclides, micro_xs, CHAIN_PATH, nuc_units='atom/cm3') + + fuel = Material(name="uo2") + fuel.add_element("U", 1, percent_type="ao", enrichment=4.25) + fuel.add_element("O", 2) + fuel.set_density("g/cc", 10.4) + fuel.depletable=True + fuel.volume = 1 + materials = Materials([fuel]) + IndependentOperator(materials, micro_xs, CHAIN_PATH) diff --git a/tests/unit_tests/test_deplete_microxs.py b/tests/unit_tests/test_deplete_microxs.py new file mode 100644 index 000000000..557082eda --- /dev/null +++ b/tests/unit_tests/test_deplete_microxs.py @@ -0,0 +1,59 @@ +"""Basic unit tests for openmc.deplete.IndependentOperator instantiation + +Modifies and resets environment variable OPENMC_CROSS_SECTIONS +to a custom file with new depletion_chain node +""" + +from os import remove +from pathlib import Path + +import pytest +from openmc.deplete import MicroXS +import numpy as np + +ONE_GROUP_XS = Path(__file__).parents[1] / "micro_xs_simple.csv" + + +def test_from_array(): + nuclides = [ + 'U234', + 'U235', + 'U238', + 'U236', + 'O16', + 'O17', + 'I135', + 'Xe135', + 'Xe136', + 'Cs135', + 'Gd157', + 'Gd156'] + reactions = ['fission', '(n,gamma)'] + # These values are placeholders and are not at all + # physically meaningful. + data = np.array([[0.1, 0.], + [0.1, 0.], + [0.9, 0.], + [0.4, 0.], + [0., 0.], + [0., 0.], + [0., 0.1], + [0., 0.9], + [0., 0.], + [0., 0.], + [0., 0.1], + [0., 0.1]]) + + MicroXS.from_array(nuclides, reactions, data) + with pytest.raises(ValueError, match=r'Nuclides list of length \d* and ' + r'reactions array of length \d* do not ' + r'match dimensions of data array of shape \(\d*\,d*\)'): + MicroXS.from_array(nuclides, reactions, data[:, 0]) + +def test_csv(): + ref_xs = MicroXS.from_csv(ONE_GROUP_XS) + ref_xs.to_csv('temp_xs.csv') + temp_xs = MicroXS.from_csv('temp_xs.csv') + assert np.all(ref_xs == temp_xs) + remove('temp_xs.csv') + diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index 6ce15601f..63073cf57 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -11,7 +11,7 @@ import openmc.deplete @pytest.fixture def res(): """Load the reference results""" - filename = (Path(__file__).parents[1] / 'regression_tests' / 'deplete' + filename = (Path(__file__).parents[1] / 'regression_tests' / 'deplete_with_transport' / 'test_reference.h5') return openmc.deplete.Results(filename) @@ -70,14 +70,16 @@ def test_get_keff(res): np.testing.assert_allclose(k[:, 1], u_ref) -@pytest.mark.parametrize("unit", ("s", "d", "min", "h")) +@pytest.mark.parametrize("unit", ("s", "d", "min", "h", "a")) def test_get_steps(unit): # Make a Results full of near-empty Result instances # Just fill out a time schedule results = openmc.deplete.Results() # Time in units of unit times = np.linspace(0, 100, num=5) - if unit == "d": + if unit == "a": + conversion_to_seconds = 60 * 60 * 24 * 365.25 + elif unit == "d": conversion_to_seconds = 60 * 60 * 24 elif unit == "h": conversion_to_seconds = 60 * 60 diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py index b57f9578a..b8bb94f37 100644 --- a/tests/unit_tests/test_stats.py +++ b/tests/unit_tests/test_stats.py @@ -6,6 +6,12 @@ import openmc import openmc.stats +def assert_sample_mean(samples, expected_mean): + std_dev = samples.std() / np.sqrt(samples.size) + assert np.abs(expected_mean - samples.mean()) < 3*std_dev + + + def test_discrete(): x = [0.0, 1.0, 10.0] p = [0.3, 0.2, 0.5] @@ -13,8 +19,8 @@ def test_discrete(): elem = d.to_xml_element('distribution') d = openmc.stats.Discrete.from_xml_element(elem) - assert d.x == x - assert d.p == p + np.testing.assert_array_equal(d.x, x) + np.testing.assert_array_equal(d.p, p) assert len(d) == len(x) d = openmc.stats.Univariate.from_xml_element(elem) @@ -33,13 +39,11 @@ def test_discrete(): d3 = openmc.stats.Discrete(vals, probs) - # sample discrete distribution + # sample discrete distribution and check that the mean of the samples is + # within 3 std. dev. of the expected mean n_samples = 1_000_000 samples = d3.sample(n_samples, seed=100) - # check that the mean of the samples is within 3 std. dev. - # of the expected mean - std_dev = samples.std() / np.sqrt(n_samples) - assert np.abs(exp_mean - samples.mean()) < 3*std_dev + assert_sample_mean(samples, exp_mean) def test_merge_discrete(): @@ -76,17 +80,16 @@ def test_uniform(): assert len(d) == 2 t = d.to_tabular() - assert t.x == [a, b] - assert t.p == [1/(b-a), 1/(b-a)] + np.testing.assert_array_equal(t.x, [a, b]) + np.testing.assert_array_equal(t.p, [1/(b-a), 1/(b-a)]) assert t.interpolation == 'histogram' + # Sample distribution and check that the mean of the samples is within 3 + # std. dev. of the expected mean exp_mean = 0.5 * (a + b) n_samples = 1_000_000 samples = d.sample(n_samples, seed=100) - # check that the mean of the samples is within 3 std. dev. - # of the expected mean - std_dev = samples.std() / np.sqrt(n_samples) - assert np.abs(exp_mean - samples.mean()) < 3*std_dev + assert_sample_mean(samples, exp_mean) def test_powerlaw(): @@ -102,13 +105,11 @@ def test_powerlaw(): exp_mean = 100.0 * (n+1) / (n+2) - # sample power law distribution + # sample power law distribution and check that the mean of the samples is + # within 3 std. dev. of the expected mean n_samples = 1_000_000 samples = d.sample(n_samples, seed=100) - # check that the mean of the samples is within 3 std. dev. - # of the expected mean - std_dev = samples.std() / np.sqrt(n_samples) - assert np.abs(exp_mean - samples.mean()) < 3*std_dev + assert_sample_mean(samples, exp_mean) def test_maxwell(): @@ -122,21 +123,15 @@ def test_maxwell(): exp_mean = 3/2 * theta - # sample maxwell distribution + # sample maxwell distribution and check that the mean of the samples is + # within 3 std. dev. of the expected mean n_samples = 1_000_000 samples = d.sample(n_samples, seed=100) - # check that the mean of the samples is within 3 std. dev. - # of the expected mean - std_dev = samples.std() / np.sqrt(n_samples) - assert np.abs(exp_mean - samples.mean()) < 3*std_dev + assert_sample_mean(samples, exp_mean) # A second sample with a different seed samples_2 = d.sample(n_samples, seed=200) - # check that the mean of the samples is within 3 std. dev. - # of the expected mean - std_dev = samples_2.std() / np.sqrt(n_samples) - assert np.abs(exp_mean - samples_2.mean()) < 3*std_dev - + assert_sample_mean(samples_2, exp_mean) assert samples_2.mean() != samples.mean() @@ -156,13 +151,11 @@ def test_watt(): # https://doi.org/10.1016/j.physletb.2003.09.048 exp_mean = 3/2 * a + a**2 * b / 4 - # sample Watt distribution + # sample Watt distribution and check that the mean of the samples is within + # 3 std. dev. of the expected mean n_samples = 1_000_000 samples = d.sample(n_samples, seed=100) - # check that the mean of the samples is within 3 std. dev. - # of the expected mean - std_dev = samples.std() / np.sqrt(n_samples) - assert np.abs(exp_mean - samples.mean()) < 3*std_dev + assert_sample_mean(samples, exp_mean) def test_tabular(): @@ -228,6 +221,11 @@ def test_mixture(): assert mix.distribution == [d1, d2] assert len(mix) == 4 + # Sample and make sure sample mean is close to expected mean + n_samples = 1_000_000 + samples = mix.sample(n_samples) + assert_sample_mean(samples, (2.5 + 5.0)/2) + elem = mix.to_xml_element('distribution') d = openmc.stats.Mixture.from_xml_element(elem) @@ -396,3 +394,47 @@ def test_muir(): assert within_2_sigma / n_samples >= 0.95 within_3_sigma = np.count_nonzero(samples < 3*d.std_dev) assert within_3_sigma / n_samples >= 0.99 + + +def test_combine_distributions(): + # Combine two discrete (same data as in test_merge_discrete) + x1 = [0.0, 1.0, 10.0] + p1 = [0.3, 0.2, 0.5] + d1 = openmc.stats.Discrete(x1, p1) + x2 = [0.5, 1.0, 5.0] + p2 = [0.4, 0.5, 0.1] + d2 = openmc.stats.Discrete(x2, p2) + + # Merged distribution should have x values sorted and probabilities + # appropriately combined. Duplicate x values should appear once. + merged = openmc.stats.combine_distributions([d1, d2], [0.6, 0.4]) + assert isinstance(merged, openmc.stats.Discrete) + assert merged.x == pytest.approx([0.0, 0.5, 1.0, 5.0, 10.0]) + assert merged.p == pytest.approx( + [0.6*0.3, 0.4*0.4, 0.6*0.2 + 0.4*0.5, 0.4*0.1, 0.6*0.5]) + + # Probabilities add up but are not normalized + d1 = openmc.stats.Discrete([3.0], [1.0]) + triple = openmc.stats.combine_distributions([d1, d1, d1], [1.0, 2.0, 3.0]) + assert triple.x == pytest.approx([3.0]) + assert triple.p == pytest.approx([6.0]) + + # Combine discrete and tabular + t1 = openmc.stats.Tabular(x2, p2) + mixed = openmc.stats.combine_distributions([d1, t1], [0.5, 0.5]) + assert isinstance(mixed, openmc.stats.Mixture) + assert len(mixed.distribution) == 2 + assert len(mixed.probability) == 2 + + # Combine 1 discrete and 2 tabular -- the tabular distributions should + # combine to produce a uniform distribution with mean 0.5. The combined + # distribution should have a mean of 0.25. + t1 = openmc.stats.Tabular([0., 1.], [2.0, 0.0]) + t2 = openmc.stats.Tabular([0., 1.], [0.0, 2.0]) + d1 = openmc.stats.Discrete([0.0], [1.0]) + combined = openmc.stats.combine_distributions([t1, t2, d1], [0.25, 0.25, 0.5]) + + # Sample the combined distribution and make sure the sample mean is within + # uncertainty of the expected value + samples = combined.sample(1000) + assert_sample_mean(samples, 0.25)