diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 34b73fc55c5..73ff862009e 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -1380,7 +1380,7 @@ mesh-based weight windows. *Default*: true :method: - Method used to update weight window values (currently only 'magic' is supported) + Method used to update weight window values (one of 'magic' or 'fw_cadis') *Default*: magic diff --git a/docs/source/methods/index.rst b/docs/source/methods/index.rst index 08aa7d05341..75c421c8773 100644 --- a/docs/source/methods/index.rst +++ b/docs/source/methods/index.rst @@ -20,4 +20,5 @@ Theory and Methodology energy_deposition parallelization cmfd + variance_reduction random_ray \ No newline at end of file diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst index c827c351dd0..eb22b0544da 100644 --- a/docs/source/methods/random_ray.rst +++ b/docs/source/methods/random_ray.rst @@ -1060,6 +1060,8 @@ random ray and Monte Carlo, however. develop the scattering source by way of inactive batches before beginning active batches. +.. _adjoint: + ------------------------ Adjoint Flux Solver Mode ------------------------ diff --git a/docs/source/methods/variance_reduction.rst b/docs/source/methods/variance_reduction.rst new file mode 100644 index 00000000000..353ae5077ea --- /dev/null +++ b/docs/source/methods/variance_reduction.rst @@ -0,0 +1,134 @@ +.. _methods_variance_reduction: + +================== +Variance Reduction +================== + +.. _methods_variance_reduction_intro: + +------------ +Introduction +------------ + +Transport problems can sometimes involve a significant degree of attenuation +between the source and a detector (tally) region, which can result in a flux +differential of ten orders of magnitude (or more) throughout the simulation +domain. As Monte Carlo uncertainties tend to be inversely proportional to the +physical flux density, it can be extremely difficult to accurately resolve +tallies in locations that are optically far from the source. This issue is +particularly common in fixed source simulations, where some tally locations may +not experience a single scoring event, even after billions of analog histories. + +Variance reduction techniques aim to either flatten the global uncertainty +distribution, such that all regions of phase space have a fairly similar +uncertainty, or to reduce the uncertainty in specific locations (such as a +detector). There are two strategies available in OpenMC for variance reduction: +the Monte Carlo MAGIC method and the FW-CADIS method. Both strategies work by +developing a weight window mesh that can be utilized by subsequent Monte Carlo +solves to split particles heading towards areas of lower flux densities while +terminating particles in higher flux regions---all while maintaining a fair +game. + +------------ +MAGIC Method +------------ + +The Method of Automatic Generation of Importances by Calculation, or `MAGIC +method `_, is an iterative +technique that uses spatial flux information :math:`\phi(r)` obtained from a +normal Monte Carlo solve to produce weight windows :math:`w(r)` that can be +utilized by a subsequent iteration of Monte Carlo. While the first generation of +weight windows produced may only help to reduce variance slightly, use of these +weights to generate another set of weight windows results in a progressively +improving iterative scheme. + +Equation :eq:`magic` defines how the lower bound of weight windows +:math:`w_{\ell}(r)` are generated with MAGIC using forward flux information. +Here, we can see that the flux at location :math:`r` is normalized by the +maximum flux in any group at that location. We can also see that the weights are +divided by a factor of two, which accounts for the typical :math:`5\times` +factor separating the lower and upper weight window bounds in OpenMC. + +.. math:: + :label: magic + + w_{\ell}(r) = \frac{\phi(r)}{2\,\text{max}(\phi(r))} + +A major advantage of this technique is that it does not require any special +transport machinery; it simply uses multiple Monte Carlo simulations to +iteratively improve a set of weight windows (which are typically defined on a +mesh covering the simulation domain). The downside to this method is that as the +flux differential increases between areas near and far from the source, it +requires more outer Monte Carlo iterations, each of which can be expensive in +itself. Additionally, computation of weight windows based on regular (forward) +neutron flux tally information does not produce the most numerically effective +set of weight windows. Nonetheless, MAGIC remains a simple and effective +technique for generating weight windows. + +-------- +FW-CADIS +-------- + +As discussed in the previous section, computation of weight windows based on +regular (forward) neutron flux tally information does not produce the most +numerically efficient set of weight windows. It is highly preferable to generate +weight windows based on spatial adjoint flux :math:`\phi^{\dag}(r)` +information. The adjoint flux is essentially the "reverse" simulation problem, +where we sample a random point and assume this is where a particle was absorbed, +and then trace it backwards (upscattering in energy), until we sample the point +where it was born from. + +The Forward-Weighted Consistent Adjoint Driven Importance Sampling method, or +`FW-CADIS method `_, produces weight windows +for global variance reduction given adjoint flux information throughout the +entire domain. The weight window lower bound is defined in Equation +:eq:`fw_cadis`, and also involves a normalization step not shown here. + +.. math:: + :label: fw_cadis + + w_{\ell}(r) = \frac{1}{2\phi^{\dag}(r)} + +While the algorithm itself is quite simple, it requires estimates of the global +adjoint flux distribution, which is difficult to generate directly with Monte +Carlo transport. Thus, FW-CADIS typically uses an alternative solver (often +deterministic) that can be more readily adapted for generating adjoint flux +information, and which is often much cheaper than Monte Carlo given that a rough +solution is often sufficient for weight window generation. + +The FW-CADIS implementation in OpenMC utilizes its own internal random ray +multigroup transport solver to generate the adjoint source distribution. No +coupling to any external transport is solver is necessary. The random ray solver +operates on the same geometry as the Monte Carlo solver, so no redefinition of +the simulation geometry is required. More details on how the adjoint flux is +computed are given in the :ref:`adjoint methods section `. + +More information on the workflow is available in the :ref:`user guide +`, but generally production of weight windows with FW-CADIS +involves several stages (some of which are highly automated). These tasks +include generation of approximate multigroup cross section data for use by the +random ray solver, running of the random ray solver in normal (forward flux) +mode to generate a source for the adjoint solver, running of the random ray +solver in adjoint mode to generate adjoint flux tallies, and finally the +production of weight windows via the FW-CADIS method. As is discussed in the +user guide, most of these steps are automated together, making the additional +burden on the user fairly small. + +The major advantage of this technique is that it typically produces much more +numerically efficient weight windows as compared to those generated with MAGIC, +sometimes with an order-of-magnitude improvement in the figure of merit +(Equation :eq:`variance_fom`), which accounts for both the variance and the +execution time. Another major advantage is that the cost of the random ray +solver is typically negligible compared to the cost of the subsequent Monte +Carlo solve itself, making it a very cheap method to deploy. The downside to +this method is that it introduces a second transport method into the mix (random +ray), such that there are more free input parameters for the user to know about +and adjust, potentially making the method more complex to use. However, as many +of the parameters have natural choices, much of this parameterization can be +handled automatically behind the scenes without the need for the user to be +aware of this. + +.. math:: + :label: variance_fom + + \text{FOM} = \frac{1}{\text{Time} \times \sigma^2} diff --git a/docs/source/usersguide/index.rst b/docs/source/usersguide/index.rst index 03a77e87063..74651c0114b 100644 --- a/docs/source/usersguide/index.rst +++ b/docs/source/usersguide/index.rst @@ -25,6 +25,7 @@ essential aspects of using OpenMC to perform simulations. processing parallel volume + variance_reduction random_ray troubleshoot \ No newline at end of file diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst index 2b9cf67240b..b6852b7ca8f 100644 --- a/docs/source/usersguide/random_ray.rst +++ b/docs/source/usersguide/random_ray.rst @@ -435,10 +435,11 @@ Inputting Multigroup Cross Sections (MGXS) Multigroup cross sections for use with OpenMC's random ray solver are input the same way as with OpenMC's traditional multigroup Monte Carlo mode. There is more information on generating multigroup cross sections via OpenMC in the -:ref:`multigroup materials ` user guide. You may also wish to -use an existing multigroup library. An example of using OpenMC's Python -interface to generate a correctly formatted ``mgxs.h5`` input file is given -in the `OpenMC Jupyter notebook collection +:ref:`multigroup materials ` user guide. You may also wish to use +an existing ``mgxs.h5`` MGXS library file, or define your own given a known set +of cross section data values (e.g., as taken from a benchmark specification). An +example of using OpenMC's Python interface to generate a correctly formatted +``mgxs.h5`` input file is given in the `OpenMC Jupyter notebook collection `_. .. note:: @@ -447,6 +448,184 @@ in the `OpenMC Jupyter notebook collection separate materials can be defined each with a separate multigroup dataset corresponding to a given temperature. +.. _mgxs_gen: + +------------------------------------------- +Generating Multigroup Cross Sections (MGXS) +------------------------------------------- + +OpenMC is capable of generating multigroup cross sections by way of flux +collapsing data based on flux solutions obtained from a continuous energy Monte +Carlo solve. While it is a circular excercise in some respects to use continuous +energy Monte Carlo to generate cross sections to be used by a reduced-fidelity +multigroup transport solver, there are many use cases where this is nonetheless +highly desirable. For instance, generation of a multigroup library may enable +the same set of approximate multigroup cross section data to be used across a +variety of problem types (or through a multidimensional parameter sweep of +design variables) with only modest errors and at greatly reduced cost as +compared to using only continuous energy Monte Carlo. + +We give here a quick summary of how to produce a multigroup cross section data +file (``mgxs.h5``) from a starting point of a typical continuous energy Monte +Carlo input file. Notably, continuous energy input files define materials as a +mixture of nuclides with different densities, whereas multigroup materials are +simply defined by which name they correspond to in a ``mgxs.h5`` library file. + +To generate the cross section data, we begin with a continuous energy Monte +Carlo input deck and add in the required tallies that will be needed to generate +our library. In this example, we will specify material-wise cross sections and a +two group energy decomposition:: + + # Define geometry + ... + ... + geometry = openmc.Geometry() + ... + ... + + # Initialize MGXS library with a finished OpenMC geometry object + mgxs_lib = openmc.mgxs.Library(geometry) + + # Pick energy group structure + groups = mgxs.EnergyGroups(mgxs.GROUP_STRUCTURES['CASMO-2']) + mgxs_lib.energy_groups = groups + + # Disable transport correction + mgxs_lib.correction = None + + # Specify needed cross sections for random ray + mgxs_lib.mgxs_types = ['total', 'absorption', 'nu-fission', 'fission', + 'nu-scatter matrix', 'multiplicity matrix', 'chi'] + + # Specify a "cell" domain type for the cross section tally filters + mgxs_lib.domain_type = "material" + + # Specify the cell domains over which to compute multi-group cross sections + mgxs_lib.domains = geom.get_all_materials().values() + + # Do not compute cross sections on a nuclide-by-nuclide basis + mgxs_lib.by_nuclide = False + + # Check the library - if no errors are raised, then the library is satisfactory. + mgxs_lib.check_library_for_openmc_mgxs() + + # Construct all tallies needed for the multi-group cross section library + mgxs_lib.build_library() + + # Create a "tallies.xml" file for the MGXS Library + tallies = openmc.Tallies() + mgxs_lib.add_to_tallies_file(tallies, merge=True) + + # Export + tallies.export_to_xml() + + ... + +When selecting an energy decomposition, you can manually define group boundaries +or pick out a group structure already known to OpenMC (a list of which can be +found at :class:`openmc.mgxs.GROUP_STRUCTURES`). Once the above input deck has +been run, the resulting statepoint file will contain the needed flux and +reaction rate tally data so that a MGXS library file can be generated. Below is +the postprocessing script needed to generate the ``mgxs.h5`` library file given +a statepoint file (e.g., ``statepoint.100.h5``) file and summary file (e.g., +``summary.h5``) that resulted from running our previous example:: + + import openmc + import openmc.mgxs as mgxs + + summary = openmc.Summary('summary.h5') + geom = summary.geometry + mats = summary.materials + + statepoint_filename = 'statepoint.100.h5' + sp = openmc.StatePoint(statepoint_filename) + + groups = mgxs.EnergyGroups(mgxs.GROUP_STRUCTURES['CASMO-2']) + mgxs_lib = openmc.mgxs.Library(geom) + mgxs_lib.energy_groups = groups + mgxs_lib.correction = None + mgxs_lib.mgxs_types = ['total', 'absorption', 'nu-fission', 'fission', + 'nu-scatter matrix', 'multiplicity matrix', 'chi'] + + # Specify a "cell" domain type for the cross section tally filters + mgxs_lib.domain_type = "material" + + # Specify the cell domains over which to compute multi-group cross sections + mgxs_lib.domains = geom.get_all_materials().values() + + # Do not compute cross sections on a nuclide-by-nuclide basis + mgxs_lib.by_nuclide = False + + # Check the library - if no errors are raised, then the library is satisfactory. + mgxs_lib.check_library_for_openmc_mgxs() + + # Construct all tallies needed for the multi-group cross section library + mgxs_lib.build_library() + + mgxs_lib.load_from_statepoint(sp) + + names = [] + for mat in mgxs_lib.domains: names.append(mat.name) + + # Create a MGXS File which can then be written to disk + mgxs_file = mgxs_lib.create_mg_library(xs_type='macro', xsdata_names=names) + + # Write the file to disk using the default filename of "mgxs.h5" + mgxs_file.export_to_hdf5("mgxs.h5") + +Notably, the postprocessing script needs to match the same +:class:`openmc.mgxs.Library` settings that were used to generate the tallies, +but otherwise is able to discern the rest of the simulation details from the +statepoint and summary files. Once the postprocessing script is successfully +run, the ``mgxs.h5`` file can be loaded by subsequent runs of OpenMC. + +If you want to convert continuous energy material objects in an OpenMC input +deck to multigroup ones from a ``mgxs.h5`` library, you can follow the below +example. Here we begin with the original continuous energy materials we used to +generate our MGXS library:: + + fuel = openmc.Material(name='UO2 (2.4%)') + fuel.set_density('g/cm3', 10.29769) + fuel.add_nuclide('U234', 4.4843e-6) + fuel.add_nuclide('U235', 5.5815e-4) + fuel.add_nuclide('U238', 2.2408e-2) + fuel.add_nuclide('O16', 4.5829e-2) + + water = openmc.Material(name='Hot borated water') + water.set_density('g/cm3', 0.740582) + water.add_nuclide('H1', 4.9457e-2) + water.add_nuclide('O16', 2.4672e-2) + water.add_nuclide('B10', 8.0042e-6) + water.add_nuclide('B11', 3.2218e-5) + water.add_s_alpha_beta('c_H_in_H2O') + + materials = openmc.Materials([fuel, water]) + +Once the ``mgxs.h5`` library file has been generated, we can then manually make +the necessary edits to the material definitions so that they load from the +multigroup library instead of defining their isotopic contents, as:: + + # Instantiate some Macroscopic Data + fuel_data = openmc.Macroscopic('UO2 (2.4%)') + water_data = openmc.Macroscopic('Hot borated water') + + # Instantiate some Materials and register the appropriate Macroscopic objects + fuel= openmc.Material(name='UO2 (2.4%)') + fuel.set_density('macro', 1.0) + fuel.add_macroscopic(fuel_data) + + water= openmc.Material(name='Hot borated water') + water.set_density('macro', 1.0) + water.add_macroscopic(water_data) + + # Instantiate a Materials collection and export to XML + materials = openmc.Materials([fuel, water]) + materials.cross_sections = "mgxs.h5" + +In the above example, our ``fuel`` and ``water`` materials will now load MGXS +data from the ``mgxs.h5`` file instead of loading continuous energy isotopic +cross section data. + -------------- Linear Sources -------------- @@ -597,9 +776,7 @@ estimator, the following code would be used: Adjoint Flux Mode ----------------- -The adjoint flux random ray solver mode can be enabled as: -entire -:: +The adjoint flux random ray solver mode can be enabled as:: settings.random_ray['adjoint'] = True diff --git a/docs/source/usersguide/variance_reduction.rst b/docs/source/usersguide/variance_reduction.rst new file mode 100644 index 00000000000..124f342034a --- /dev/null +++ b/docs/source/usersguide/variance_reduction.rst @@ -0,0 +1,162 @@ +.. _variance_reduction: + +================== +Variance Reduction +================== + +Global variance reduction in OpenMC is accomplished by weight windowing +techniques. OpenMC is capable of generating weight windows using either the +MAGIC or FW-CADIS methods. Both techniques will produce a ``weight_windows.h5`` +file that can be loaded and used later on. In this section, we break down the +steps required to both generate and then apply weight windows. + +.. _ww_generator: + +------------------------------------ +Generating Weight Windows with MAGIC +------------------------------------ + +As discussed in the :ref:`methods section `, MAGIC +is an iterative method that uses flux tally information from a Monte Carlo +simulation to produce weight windows for a user-defined mesh. While generating +the weight windows, OpenMC is capable of applying the weight windows generated +from a previous batch while processing the next batch, allowing for progressive +improvement in the weight window quality across iterations. + +The typical way of generating weight windows is to define a mesh and then add an +:class:`openmc.WeightWindowGenerator` object to an :attr:`openmc.Settings` +instance, as follows:: + + # Define weight window spatial mesh + ww_mesh = openmc.RegularMesh() + ww_mesh.dimension = (10, 10, 10) + ww_mesh.lower_left = (0.0, 0.0, 0.0) + ww_mesh.upper_right = (100.0, 100.0, 100.0) + + # Create weight window object and adjust parameters + wwg = openmc.WeightWindowGenerator( + method='magic', + mesh=ww_mesh, + max_realizations=settings.batches + ) + + # Add generator to Settings instance + settings.weight_window_generators = wwg + +Notably, the :attr:`max_realizations` attribute is adjusted to the number of +batches, such that all iterations are used to refine the weight window +parameters. + +With the :class:`~openmc.WeightWindowGenerator` instance added to the +:attr:`~openmc.Settings`, the rest of the problem can be defined as normal. When +running, note that the second iteration and beyond may be several orders of +magnitude slower than the first. As the weight windows are applied in each +iteration, particles may be agressively split, resulting in a large number of +secondary (split) particles being generated per initial source particle. This is +not necessarily a bad thing, as the split particles are much more efficient at +exploring low flux regions of phase space as compared to initial particles. +Thus, even though the reported "particles/second" metric of OpenMC may be much +lower when generating (or just applying) weight windows as compared to analog +MC, it typically leads to an overall improvement in the figure of merit +accounting for the reduction in the variance. + +.. warning:: + The number of particles per batch may need to be adjusted downward + significantly to result in reasonable runtimes when weight windows are being + generated or used. + +At the end of the simulation, a ``weight_windows.h5`` file will be saved to disk +for later use. Loading it in another subsequent simulation will be discussed in +the "Using Weight Windows" section below. + +------------------------------------------------------ +Generating Weight Windows with FW-CADIS and Random Ray +------------------------------------------------------ + +Weight window generation with FW-CADIS and random ray in OpenMC uses the same +exact strategy as with MAGIC. An :class:`openmc.WeightWindowGenerator` object is +added to the :attr:`openmc.Settings` object, and a ``weight_windows.h5`` will be +generated at the end of the simulation. The only difference is that the code +must be run in random ray mode. A full description of how to enable and setup +random ray mode can be found in the :ref:`Random Ray User Guide `. + +.. note:: + It is a long term goal for OpenMC to be able to generate FW-CADIS weight + windows with only a few tweaks to an existing continuous energy Monte Carlo + input deck. However, at the present time, the workflow requires several + steps to generate multigroup cross section data and to configure the random + ray solver. A high level overview of the current workflow for generation of + weight windows with FW-CADIS using random ray is given below. + +1. Produce approximate multigroup cross section data (stored in a ``mgxs.h5`` + library). There is more information on generating multigroup cross sections + via OpenMC in the :ref:`multigroup materials ` user guide, and a + specific example of generating cross section data for use with random ray in + the :ref:`random ray MGXS guide `. + +2. Make a copy of your continuous energy Python input file. You'll edit the new + file to work in multigroup mode with random ray for producing weight windows. + +3. Adjust the material definitions in your new multigroup Python file to utilize + the multigroup cross sections instead of nuclide-wise continuous energy data. + There is a specific example of making this conversion in the :ref:`random ray + MGXS guide `. + +4. Configure OpenMC to run in random ray mode (by adding several standard random + ray input flags and settings to the :attr:`openmc.Settings.random_ray` + dictionary). More information can be found in the :ref:`Random Ray User + Guide `. + +5. Add in a :class:`~openmc.WeightWindowGenerator` in a similar manner as for + MAGIC generation with Monte Carlo and set the :attr:`method` attribute set to + ``"fw_cadis"``:: + + # Define weight window spatial mesh + ww_mesh = openmc.RegularMesh() + ww_mesh.dimension = (10, 10, 10) + ww_mesh.lower_left = (0.0, 0.0, 0.0) + ww_mesh.upper_right = (100.0, 100.0, 100.0) + + # Create weight window object and adjust parameters + wwg = openmc.WeightWindowGenerator( + method='fw_cadis', + mesh=ww_mesh, + max_realizations=settings.batches + ) + + # Add generator to openmc.settings object + settings.weight_window_generators = wwg + + +.. warning:: + If using FW-CADIS weight window generation, ensure that the selected weight + window mesh does not subdivide any cells in the problem. In the future, this + restriction is intended to be relaxed, but for now subdivision of cells by a + mesh tally will result in undefined behavior. + +6. When running your multigroup random ray input deck, OpenMC will automatically + run a forward solve followed by an adjoint solve, with a + ``weight_windows.h5`` file generated at the end. The ``weight_windows.h5`` + file will contain FW-CADIS generated weight windows. This file can be used in + identical manner as one generated with MAGIC, as described below. + +-------------------- +Using Weight Windows +-------------------- + +To use a ``weight_windows.h5`` weight window file with OpenMC's Monte Carlo +solver, the Python input just needs to load the h5 file:: + + settings.weight_window_checkpoints = {'collision': True, 'surface': True} + settings.survival_biasing = False + settings.weight_windows = openmc.hdf5_to_wws('weight_windows.h5') + settings.weight_windows_on = True + +The :class:`~openmc.WeightWindowGenerator` instance is not needed to load an +existing ``weight_windows.h5`` file. Inclusion of a +:class:`~openmc.WeightWindowGenerator` instance will cause OpenMC to generate +*new* weight windows and thus overwrite the existing ``weight_windows.h5`` file. +Weight window mesh information is embedded into the weight window file, so the +mesh does not need to be redefined. Monte Carlo solves that load a weight window +file as above will utilize weight windows to reduce the variance of the +simulation. diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index 9852fe679d0..2d4d556948b 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -17,9 +17,7 @@ namespace openmc { -enum class WeightWindowUpdateMethod { - MAGIC, -}; +enum class WeightWindowUpdateMethod { MAGIC, FW_CADIS }; //============================================================================== // Constants @@ -127,8 +125,9 @@ class WeightWindows { //! \param[in] threshold Relative error threshold. Results over this //! threshold will be ignored \param[in] ratio Ratio of upper to lower //! weight window bounds - void update_magic(const Tally* tally, const std::string& value = "mean", - double threshold = 1.0, double ratio = 5.0); + void update_weights(const Tally* tally, const std::string& value = "mean", + double threshold = 1.0, double ratio = 5.0, + WeightWindowUpdateMethod method = WeightWindowUpdateMethod::MAGIC); // NOTE: This is unused for now but may be used in the future //! Write weight window settings to an HDF5 file @@ -221,12 +220,11 @@ class WeightWindowsGenerator { void create_tally(); // Data members - int32_t tally_idx_; //!< Index of the tally used to update the weight windows - int32_t ww_idx_; //!< Index of the weight windows object being generated - std::string method_; //!< Method used to update weight window. Only "magic" - //!< is valid for now. - int32_t max_realizations_; //!< Maximum number of tally realizations - int32_t update_interval_; //!< Determines how often updates occur + int32_t tally_idx_; //!< Index of the tally used to update the weight windows + int32_t ww_idx_; //!< Index of the weight windows object being generated + WeightWindowUpdateMethod method_; //!< Method used to update weight window. + int32_t max_realizations_; //!< Maximum number of tally realizations + int32_t update_interval_; //!< Determines how often updates occur bool on_the_fly_; //!< Whether or not to keep tally results between batches or //!< realizations diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index 1ddbfe10fe1..0502d309f5e 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -660,9 +660,8 @@ class WeightWindowGenerator: maximum and minimum energy for the data available at runtime. particle_type : {'neutron', 'photon'} Particle type the weight windows apply to - method : {'magic'} - The weight window generation methodology applied during an update. Only - 'magic' is currently supported. + method : {'magic', 'fw_cadis'} + The weight window generation methodology applied during an update. max_realizations : int The upper limit for number of tally realizations when generating weight windows. @@ -680,9 +679,8 @@ class WeightWindowGenerator: energies in [eV] for a single bin particle_type : {'neutron', 'photon'} Particle type the weight windows apply to - method : {'magic'} - The weight window generation methodology applied during an update. Only - 'magic' is currently supported. + method : {'magic', 'fw_cadis'} + The weight window generation methodology applied during an update. max_realizations : int The upper limit for number of tally realizations when generating weight windows. @@ -767,7 +765,7 @@ def method(self) -> str: @method.setter def method(self, m: str): cv.check_type('generation method', m, str) - cv.check_value('generation method', m, {'magic'}) + cv.check_value('generation method', m, ('magic', 'fw_cadis')) self._method = m if self._update_parameters is not None: try: @@ -800,7 +798,7 @@ def update_parameters(self) -> dict: return self._update_parameters def _check_update_parameters(self, params: dict): - if self.method == 'magic': + if self.method == 'magic' or self.method == 'fw_cadis': check_params = self._MAGIC_PARAMS for key, val in params.items(): @@ -843,7 +841,7 @@ def _sanitize_update_parameters(cls, method: str, update_parameters: dict): update_parameters : dict The update parameters as-read from the XML node (keys: str, values: str) """ - if method == 'magic': + if method == 'magic' or method == 'fw_cadis': check_params = cls._MAGIC_PARAMS for param, param_type in check_params.items(): diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 3cc0532d2b1..e967afd7995 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -28,12 +28,12 @@ void collision_mg(Particle& p) // Add to the collision counter for the particle p.n_collision()++; - if (settings::weight_window_checkpoint_collision) - apply_weight_windows(p); - // Sample the reaction type sample_reaction(p); + if (settings::weight_window_checkpoint_collision) + apply_weight_windows(p); + // Display information about collision if ((settings::verbosity >= 10) || p.trace()) { write_message(fmt::format(" Energy Group = {}", p.g()), 1); diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index a9180c68e7b..a8e6d1b8229 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -173,6 +173,7 @@ void validate_random_ray_inputs() case FilterType::MATERIAL: case FilterType::MESH: case FilterType::UNIVERSE: + case FilterType::PARTICLE: break; default: fatal_error("Invalid filter specified. Only cell, cell_instance, " diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index e798a2f7abd..419b7620372 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -22,6 +22,7 @@ #include "openmc/particle.h" #include "openmc/particle_data.h" #include "openmc/physics_common.h" +#include "openmc/random_ray/flat_source_domain.h" #include "openmc/search.h" #include "openmc/settings.h" #include "openmc/tallies/filter_energy.h" @@ -482,8 +483,8 @@ void WeightWindows::set_bounds( upper_ww_ *= ratio; } -void WeightWindows::update_magic( - const Tally* tally, const std::string& value, double threshold, double ratio) +void WeightWindows::update_weights(const Tally* tally, const std::string& value, + double threshold, double ratio, WeightWindowUpdateMethod method) { /////////////////////////// // Setup and checks @@ -624,20 +625,44 @@ void WeightWindows::update_magic( auto mesh_vols = this->mesh()->volumes(); int e_bins = new_bounds.shape()[0]; - for (int e = 0; e < e_bins; e++) { - // select all - auto group_view = xt::view(new_bounds, e); - // divide by volume of mesh elements - for (int i = 0; i < group_view.size(); i++) { - group_view[i] /= mesh_vols[i]; + if (method == WeightWindowUpdateMethod::MAGIC) { + // If we are computing weight windows with forward fluxes derived from a + // Monte Carlo or forward random ray solve, we use the MAGIC algorithm. + for (int e = 0; e < e_bins; e++) { + // select all + auto group_view = xt::view(new_bounds, e); + + // divide by volume of mesh elements + for (int i = 0; i < group_view.size(); i++) { + group_view[i] /= mesh_vols[i]; + } + + double group_max = + *std::max_element(group_view.begin(), group_view.end()); + // normalize values in this energy group by the maximum value for this + // group + if (group_max > 0.0) + group_view /= 2.0 * group_max; } + } else { + // If we are computing weight windows with adjoint fluxes derived from an + // adjoint random ray solve, we use the FW-CADIS algorithm. + for (int e = 0; e < e_bins; e++) { + // select all + auto group_view = xt::view(new_bounds, e); + + // divide by volume of mesh elements + for (int i = 0; i < group_view.size(); i++) { + group_view[i] /= mesh_vols[i]; + } + } + + xt::noalias(new_bounds) = 1.0 / new_bounds; - double group_max = *std::max_element(group_view.begin(), group_view.end()); - // normalize values in this energy group by the maximum value for this - // group - if (group_max > 0.0) - group_view /= 2.0 * group_max; + auto max_val = xt::amax(new_bounds)(); + + xt::noalias(new_bounds) = new_bounds / (2.0 * max_val); } // make sure that values where the mean is zero are set s.t. the weight window @@ -760,37 +785,51 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node) e_bounds.push_back(data::energy_max[p_type]); } - // set method and parameters for updates - method_ = get_node_value(node, "method"); - if (method_ == "magic") { - // parse non-default update parameters if specified - if (check_for_node(node, "update_parameters")) { - pugi::xml_node params_node = node.child("update_parameters"); - if (check_for_node(params_node, "value")) - tally_value_ = get_node_value(params_node, "value"); - if (check_for_node(params_node, "threshold")) - threshold_ = std::stod(get_node_value(params_node, "threshold")); - if (check_for_node(params_node, "ratio")) { - ratio_ = std::stod(get_node_value(params_node, "ratio")); - } + // set method + std::string method_string = get_node_value(node, "method"); + if (method_string == "magic") { + method_ = WeightWindowUpdateMethod::MAGIC; + if (settings::solver_type == SolverType::RANDOM_RAY && + FlatSourceDomain::adjoint_) { + fatal_error("Random ray weight window generation with MAGIC cannot be " + "done in adjoint mode."); } - // check update parameter values - if (tally_value_ != "mean" && tally_value_ != "rel_err") { - fatal_error(fmt::format("Unsupported tally value '{}' specified for " - "weight window generation.", - tally_value_)); + } else if (method_string == "fw_cadis") { + method_ = WeightWindowUpdateMethod::FW_CADIS; + if (settings::solver_type != SolverType::RANDOM_RAY) { + fatal_error("FW-CADIS can only be run in random ray solver mode."); } - if (threshold_ <= 0.0) - fatal_error(fmt::format("Invalid relative error threshold '{}' (<= 0.0) " - "specified for weight window generation", - ratio_)); - if (ratio_ <= 1.0) - fatal_error(fmt::format("Invalid weight window ratio '{}' (<= 1.0) " - "specified for weight window generation")); + FlatSourceDomain::adjoint_ = true; } else { fatal_error(fmt::format( - "Unknown weight window update method '{}' specified", method_)); + "Unknown weight window update method '{}' specified", method_string)); + } + + // parse non-default update parameters if specified + if (check_for_node(node, "update_parameters")) { + pugi::xml_node params_node = node.child("update_parameters"); + if (check_for_node(params_node, "value")) + tally_value_ = get_node_value(params_node, "value"); + if (check_for_node(params_node, "threshold")) + threshold_ = std::stod(get_node_value(params_node, "threshold")); + if (check_for_node(params_node, "ratio")) { + ratio_ = std::stod(get_node_value(params_node, "ratio")); + } + } + + // check update parameter values + if (tally_value_ != "mean" && tally_value_ != "rel_err") { + fatal_error(fmt::format("Unsupported tally value '{}' specified for " + "weight window generation.", + tally_value_)); } + if (threshold_ <= 0.0) + fatal_error(fmt::format("Invalid relative error threshold '{}' (<= 0.0) " + "specified for weight window generation", + ratio_)); + if (ratio_ <= 1.0) + fatal_error(fmt::format("Invalid weight window ratio '{}' (<= 1.0) " + "specified for weight window generation")); // create a matching weight windows object auto wws = WeightWindows::create(); @@ -860,7 +899,7 @@ void WeightWindowsGenerator::update() const tally->n_realizations_ % update_interval_ != 0) return; - wws->update_magic(tally, tally_value_, threshold_, ratio_); + wws->update_weights(tally, tally_value_, threshold_, ratio_, method_); // if we're not doing on the fly generation, reset the tally results once // we're done with the update @@ -944,7 +983,7 @@ extern "C" int openmc_weight_windows_update_magic(int32_t ww_idx, // get the WeightWindows object const auto& wws = variance_reduction::weight_windows.at(ww_idx); - wws->update_magic(tally, value, threshold, ratio); + wws->update_weights(tally, value, threshold, ratio); return 0; } diff --git a/tests/regression_tests/weightwindows_fw_cadis/__init__.py b/tests/regression_tests/weightwindows_fw_cadis/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/weightwindows_fw_cadis/inputs_true.dat b/tests/regression_tests/weightwindows_fw_cadis/inputs_true.dat new file mode 100644 index 00000000000..586bbb12026 --- /dev/null +++ b/tests/regression_tests/weightwindows_fw_cadis/inputs_true.dat @@ -0,0 +1,260 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 12 12 12 + 0.0 0.0 0.0 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 + + + + + + + + + + fixed source + 90 + 10 + 5 + + + 100.0 1.0 + + + universe + 1 + + + multi-group + + + 1 + neutron + 10 + 1 + true + fw_cadis + + + + 6 6 6 + 0.0 0.0 0.0 + 30.0 30.0 30.0 + + + 500.0 + 100.0 + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + True + naive + + + + + 1 + + + 2 + + + 3 + + + 3 + flux + tracklength + + + 2 + flux + tracklength + + + 1 + flux + tracklength + + + diff --git a/tests/regression_tests/weightwindows_fw_cadis/results_true.dat b/tests/regression_tests/weightwindows_fw_cadis/results_true.dat new file mode 100644 index 00000000000..e39b2e23900 --- /dev/null +++ b/tests/regression_tests/weightwindows_fw_cadis/results_true.dat @@ -0,0 +1,442 @@ +RegularMesh + ID = 1 + Name = + Dimensions = 3 + Voxels = [6 6 6] + Lower left = [0. 0. 0.] + Upper Right = [np.float64(30.0), np.float64(30.0), np.float64(30.0)] + Width = [5. 5. 5.] +Lower Bounds +1.50e-01 +1.54e-01 +1.43e-01 +1.30e-01 +1.30e-01 +1.97e-01 +1.46e-01 +1.25e-01 +1.34e-01 +1.32e-01 +1.27e-01 +3.40e-01 +1.36e-01 +1.13e-01 +1.09e-01 +1.19e-01 +1.25e-01 +3.08e-01 +1.38e-01 +1.18e-01 +9.49e-02 +1.19e-01 +1.35e-01 +1.96e-01 +1.12e-01 +9.60e-02 +9.33e-02 +9.74e-02 +1.37e-01 +2.76e-01 +4.53e-02 +1.31e-01 +2.47e-01 +7.36e-02 +1.72e-01 +1.38e-01 +1.55e-01 +1.54e-01 +1.48e-01 +1.36e-01 +1.29e-01 +1.88e-01 +1.40e-01 +1.51e-01 +1.28e-01 +1.18e-01 +1.29e-01 +2.86e-01 +1.42e-01 +1.26e-01 +1.28e-01 +1.16e-01 +1.22e-01 +2.81e-01 +1.38e-01 +1.22e-01 +1.04e-01 +9.81e-02 +1.08e-01 +1.97e-01 +1.34e-01 +1.32e-01 +1.10e-01 +8.69e-02 +1.03e-01 +1.33e-01 +2.51e-01 +3.39e-01 +1.13e-01 +6.92e-02 +1.54e-01 +5.71e-02 +1.60e-01 +1.52e-01 +1.34e-01 +1.32e-01 +1.28e-01 +2.53e-01 +1.33e-01 +1.35e-01 +1.15e-01 +1.32e-01 +1.37e-01 +3.13e-01 +1.27e-01 +1.36e-01 +1.35e-01 +1.10e-01 +1.08e-01 +8.81e-02 +1.51e-01 +1.32e-01 +1.13e-01 +9.89e-02 +9.25e-02 +5.04e-02 +1.40e-01 +1.56e-01 +1.40e-01 +9.19e-02 +6.54e-02 +2.05e-02 +3.35e-01 +1.88e-01 +3.16e-01 +1.80e-01 +5.78e-02 +5.36e-03 +1.43e-01 +1.35e-01 +1.47e-01 +1.33e-01 +1.47e-01 +4.95e-01 +1.42e-01 +1.45e-01 +1.27e-01 +1.31e-01 +1.26e-01 +1.43e-01 +1.49e-01 +1.23e-01 +1.24e-01 +1.13e-01 +9.69e-02 +1.07e-01 +1.65e-01 +1.50e-01 +1.24e-01 +1.07e-01 +8.90e-02 +8.81e-02 +1.48e-01 +1.58e-01 +1.20e-01 +7.93e-02 +7.19e-02 +3.94e-02 +2.29e-01 +1.95e-01 +2.66e-01 +1.51e-01 +3.52e-02 +3.68e-03 +1.48e-01 +1.18e-01 +1.47e-01 +1.41e-01 +1.37e-01 +2.71e-01 +1.70e-01 +1.70e-01 +1.29e-01 +1.33e-01 +1.28e-01 +2.79e-01 +1.84e-01 +1.42e-01 +1.18e-01 +1.25e-01 +1.26e-01 +1.64e-01 +1.73e-01 +1.37e-01 +1.12e-01 +9.01e-02 +1.13e-01 +1.68e-01 +1.60e-01 +1.52e-01 +1.30e-01 +8.71e-02 +7.72e-02 +3.76e-02 +1.57e-01 +9.31e-02 +1.63e-01 +6.80e-02 +6.22e-02 +3.37e-03 +2.47e-01 +2.17e-01 +1.48e-01 +2.01e-01 +1.54e-01 +5.69e-02 +2.89e-01 +3.59e-01 +1.17e-01 +8.64e-02 +3.34e-01 +1.50e-01 +5.00e-01 +2.16e-01 +1.22e-01 +1.66e-01 +1.44e-01 +4.98e-02 +4.28e-01 +1.17e-01 +8.73e-02 +1.11e-01 +1.09e-01 +9.69e-02 +3.26e-01 +2.28e-01 +2.49e-01 +3.55e-02 +1.69e-02 +3.95e-03 +1.03e-01 +4.89e-02 +1.90e-01 +4.61e-02 +9.75e-03 +3.00e-04 +Upper Bounds +7.49e-01 +7.71e-01 +7.13e-01 +6.51e-01 +6.48e-01 +9.85e-01 +7.31e-01 +6.26e-01 +6.69e-01 +6.59e-01 +6.33e-01 +1.70e+00 +6.82e-01 +5.67e-01 +5.45e-01 +5.94e-01 +6.24e-01 +1.54e+00 +6.92e-01 +5.91e-01 +4.75e-01 +5.96e-01 +6.77e-01 +9.82e-01 +5.59e-01 +4.80e-01 +4.67e-01 +4.87e-01 +6.84e-01 +1.38e+00 +2.26e-01 +6.55e-01 +1.24e+00 +3.68e-01 +8.59e-01 +6.88e-01 +7.73e-01 +7.70e-01 +7.42e-01 +6.78e-01 +6.46e-01 +9.40e-01 +7.00e-01 +7.57e-01 +6.39e-01 +5.92e-01 +6.43e-01 +1.43e+00 +7.12e-01 +6.29e-01 +6.40e-01 +5.78e-01 +6.09e-01 +1.40e+00 +6.90e-01 +6.10e-01 +5.21e-01 +4.91e-01 +5.42e-01 +9.87e-01 +6.71e-01 +6.60e-01 +5.52e-01 +4.35e-01 +5.14e-01 +6.65e-01 +1.26e+00 +1.70e+00 +5.63e-01 +3.46e-01 +7.72e-01 +2.86e-01 +8.00e-01 +7.58e-01 +6.70e-01 +6.59e-01 +6.41e-01 +1.26e+00 +6.65e-01 +6.76e-01 +5.76e-01 +6.60e-01 +6.86e-01 +1.57e+00 +6.35e-01 +6.81e-01 +6.75e-01 +5.51e-01 +5.39e-01 +4.40e-01 +7.56e-01 +6.62e-01 +5.67e-01 +4.95e-01 +4.63e-01 +2.52e-01 +6.98e-01 +7.78e-01 +6.99e-01 +4.59e-01 +3.27e-01 +1.02e-01 +1.67e+00 +9.38e-01 +1.58e+00 +9.00e-01 +2.89e-01 +2.68e-02 +7.16e-01 +6.75e-01 +7.35e-01 +6.63e-01 +7.34e-01 +2.47e+00 +7.11e-01 +7.23e-01 +6.34e-01 +6.54e-01 +6.32e-01 +7.14e-01 +7.45e-01 +6.15e-01 +6.19e-01 +5.64e-01 +4.85e-01 +5.37e-01 +8.26e-01 +7.51e-01 +6.21e-01 +5.36e-01 +4.45e-01 +4.40e-01 +7.40e-01 +7.92e-01 +5.99e-01 +3.96e-01 +3.60e-01 +1.97e-01 +1.14e+00 +9.76e-01 +1.33e+00 +7.55e-01 +1.76e-01 +1.84e-02 +7.42e-01 +5.90e-01 +7.37e-01 +7.07e-01 +6.84e-01 +1.36e+00 +8.52e-01 +8.49e-01 +6.46e-01 +6.63e-01 +6.42e-01 +1.40e+00 +9.19e-01 +7.08e-01 +5.89e-01 +6.24e-01 +6.28e-01 +8.21e-01 +8.64e-01 +6.87e-01 +5.58e-01 +4.50e-01 +5.64e-01 +8.40e-01 +7.99e-01 +7.62e-01 +6.51e-01 +4.36e-01 +3.86e-01 +1.88e-01 +7.85e-01 +4.66e-01 +8.13e-01 +3.40e-01 +3.11e-01 +1.68e-02 +1.24e+00 +1.09e+00 +7.40e-01 +1.00e+00 +7.72e-01 +2.85e-01 +1.44e+00 +1.79e+00 +5.85e-01 +4.32e-01 +1.67e+00 +7.50e-01 +2.50e+00 +1.08e+00 +6.10e-01 +8.32e-01 +7.19e-01 +2.49e-01 +2.14e+00 +5.87e-01 +4.36e-01 +5.57e-01 +5.47e-01 +4.84e-01 +1.63e+00 +1.14e+00 +1.25e+00 +1.78e-01 +8.43e-02 +1.97e-02 +5.17e-01 +2.44e-01 +9.49e-01 +2.31e-01 +4.87e-02 +1.50e-03 \ No newline at end of file diff --git a/tests/regression_tests/weightwindows_fw_cadis/test.py b/tests/regression_tests/weightwindows_fw_cadis/test.py new file mode 100644 index 00000000000..7a4718ed5b2 --- /dev/null +++ b/tests/regression_tests/weightwindows_fw_cadis/test.py @@ -0,0 +1,33 @@ +import os + +import openmc +from openmc.examples import random_ray_three_region_cube + +from tests.testing_harness import WeightWindowPyAPITestHarness + + +class MGXSTestHarness(WeightWindowPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + + +def test_random_ray_adjoint_fixed_source(): + model = random_ray_three_region_cube() + + ww_mesh = openmc.RegularMesh() + n = 6 + width = 30.0 + ww_mesh.dimension = (n, n, n) + ww_mesh.lower_left = (0.0, 0.0, 0.0) + ww_mesh.upper_right = (width, width, width) + + wwg = openmc.WeightWindowGenerator( + method="fw_cadis", mesh=ww_mesh, max_realizations=model.settings.batches) + model.settings.weight_window_generators = wwg + model.settings.random_ray['volume_estimator'] = 'naive' + + harness = MGXSTestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/testing_harness.py b/tests/testing_harness.py index 81527452b84..ddae89e02e9 100644 --- a/tests/testing_harness.py +++ b/tests/testing_harness.py @@ -441,6 +441,40 @@ def _compare_results(self): assert compare, 'Results do not agree' +class WeightWindowPyAPITestHarness(PyAPITestHarness): + def _get_results(self): + """Digest info in the weight window file and return as a string.""" + ww = openmc.hdf5_to_wws()[0] + + # Access the weight window bounds + lower_bound = ww.lower_ww_bounds + upper_bound = ww.upper_ww_bounds + + # Flatten both arrays + flattened_lower_bound = lower_bound.flatten() + flattened_upper_bound = upper_bound.flatten() + + # Convert each element to a string in scientific notation with 2 decimal places + formatted_lower_bound = [f'{x:.2e}' for x in flattened_lower_bound] + formatted_upper_bound = [f'{x:.2e}' for x in flattened_upper_bound] + + # Concatenate the formatted arrays + concatenated_strings = ["Lower Bounds"] + formatted_lower_bound + \ + ["Upper Bounds"] + formatted_upper_bound + + # Join the concatenated strings into a single string with newline characters + final_string = '\n'.join(concatenated_strings) + + # Prepend the mesh text description and return final string + return str(ww.mesh) + final_string + + def _cleanup(self): + super()._cleanup() + f = 'weight_windows.h5' + if os.path.exists(f): + os.remove(f) + + class PlotTestHarness(TestHarness): """Specialized TestHarness for running OpenMC plotting tests.""" def __init__(self, plot_names, voxel_convert_checks=[]):