diff --git a/.gitmodules b/.gitmodules index d492a3188..f776a8fa6 100644 --- a/.gitmodules +++ b/.gitmodules @@ -40,3 +40,6 @@ [submodule "contrib/embree"] path = contrib/embree url = https://github.com/RenderKit/embree.git +[submodule "contrib/nuclear_data"] + path = contrib/nuclear_data + url = https://github.com/OTU-Centre-for-SMRs/nuclear_data.git diff --git a/Makefile b/Makefile index e799a132d..5db165036 100644 --- a/Makefile +++ b/Makefile @@ -186,6 +186,8 @@ OPENMC_INCLUDES := -I$(OPENMC_INSTALL_DIR)/include OPENMC_LIBDIR := $(OPENMC_INSTALL_DIR)/lib OPENMC_LIB := $(OPENMC_LIBDIR)/libopenmc.so +ADDITIONAL_INCLUDES += -I$(CONTRIB_DIR)/nuclear_data + # This is used in $(FRAMEWORK_DIR)/build.mk HDF5_INCLUDES := -I$(HDF5_INCLUDE_DIR) -I$(HDF5_ROOT)/include diff --git a/config/check_deps.mk b/config/check_deps.mk index 72ff1e495..af36ca502 100644 --- a/config/check_deps.mk +++ b/config/check_deps.mk @@ -6,6 +6,7 @@ endef # Set default values for all third party dependencies NEKRS_DIR ?= $(CONTRIB_DIR)/nekRS OPENMC_DIR ?= $(CONTRIB_DIR)/openmc +NUCLEARDATA_DIR ?= $(CONTRIB_DIR)/nuclear_data DAGMC_DIR ?= $(CONTRIB_DIR)/DAGMC DOUBLEDOWN_DIR ?= $(CONTRIB_DIR)/double-down EMBREE_DIR ?= $(CONTRIB_DIR)/embree @@ -23,6 +24,7 @@ IAPWS95_DIR ?= $(CONTRIB_DIR)/iapws95 MOOSE_CONTENT := $(shell ls $(MOOSE_DIR) 2> /dev/null) NEKRS_CONTENT := $(shell ls $(NEKRS_DIR) 2> /dev/null) OPENMC_CONTENT := $(shell ls $(OPENMC_DIR) 2> /dev/null) +NUCLEARDATA_CONTENT:= $(shell ls $(NUCLEARDATA_DIR) 2> /dev/null) DAGMC_CONTENT := $(shell ls $(DAGMC_DIR) 2> /dev/null) DOUBLEDOWN_CONTENT := $(shell ls $(DOUBLEDOWN_DIR) 2> /dev/null) EMBREE_CONTENT := $(shell ls $(EMBREE_DIR) 2> /dev/null) @@ -90,6 +92,12 @@ ifeq ($(ENABLE_OPENMC), yes) $(info Cardinal is using OpenMC from $(OPENMC_DIR)) endif + ifeq ($(NUCLEARDATA_CONTENT),) + $(error $n"nuclear_data does not seem to be available, but ENABLE_OPENMC is enabled. Make sure that the submodule is checked out.$n$nTo fetch the nuclear_data submodule, use ./scripts/get-dependencies.sh") + else + # we dont print out anything about where nuclear_data is coming from because it's a minor dependency and we don't expect anyone to be using different versions of it (and we don't want to confuse them on the cross section library data) + endif + openmc_status := $(shell git -C $(CONTRIB_DIR) submodule status 2>/dev/null | grep openmc | cut -c1) ifneq (,$(findstring +,$(openmc_status))) $(warning $n"***WARNING***: Your OpenMC submodule is not pointing to the commit tied to Cardinal.$n To fetch the paired commit, use ./scripts/get-dependencies.sh"$n) diff --git a/contrib/nuclear_data b/contrib/nuclear_data new file mode 160000 index 000000000..a2ec1bcaa --- /dev/null +++ b/contrib/nuclear_data @@ -0,0 +1 @@ +Subproject commit a2ec1bcaa21ee888e18ef0a99267595ed4438c68 diff --git a/doc/content/publications.md b/doc/content/publications.md index 4df5083d2..11e573937 100644 --- a/doc/content/publications.md +++ b/doc/content/publications.md @@ -145,6 +145,10 @@ H. Park, Y. Yu, E. Shemon, and A. Novak, # Sodium Fast Reactors +A.J. Novak, C. Bourdot Dutra, D. Shaver, and E. Merzari, +["CFD Simulations of Interassembly Bypass Flow in Sodium Fast Reactors"](https://www.sciencedirect.com/science/article/pii/S0029549325002213) +*Nuclear Engineering and Design* (2025) + A.J. Novak, C. Bourdot Dutra, D. Shaver, and E. Merzari, ["NekRS CFD Simulations of Interassembly Flow in Sodium Fast Reactors"](https://www.researchgate.net/publication/373331042_NekRS_CFD_Simulations_of_Interassembly_Flow_in_Sodium_Fast_Reactors) *Proceedings of Nureth* (2023) diff --git a/doc/content/source/actions/AddCriticalitySearchAction.md b/doc/content/source/actions/AddCriticalitySearchAction.md new file mode 100644 index 000000000..2f14b96d8 --- /dev/null +++ b/doc/content/source/actions/AddCriticalitySearchAction.md @@ -0,0 +1,56 @@ +# AddCriticalitySearchAction + +## Overview + +The `AddCriticalitySearchAction` is responsible for performing a complete criticality +search during the OpenMC solve. Examples of parameters which can be used for the +criticality search include: + +- [OpenMCMaterialDensity](OpenMCMaterialDensity.md), to change a material's total density +- [BoratedWater](BoratedWater.md), to change the boron weight ppm in water + +The converged value of the criticality search will automatically be populated into +a postprocessor named `critical_value`. + +!alert warning +If the convergence `tolerance` is too small (on the order of or smaller than the statistical standard deviation), +convergence may not be possible. If the search fails to converge, use a looser `tolerance` or +increase the number of particles. + +## Example Input File Syntax + +As an example, a [OpenMCMaterialDensity](OpenMCMaterialDensity.md) object is used to perform +a criticality search based on the total density of material 1. + +!listing /tests/criticality/material_density/openmc.i + block=Problem + +When running this case, Cardinal will print out a table showing a summary of the results of +the criticality search. + +``` +--------------------------------------------------------------------------- +| Iteration | material 1 density [kg/m3] | k (mean) | k (std dev) | +--------------------------------------------------------------------------- +| 0 | 1.000000e+03 | 6.620000e-02 | 1.098139e-03 | +| 1 | 3.000000e+04 | 1.621666e+00 | 2.077345e-02 | +| 2 | 1.840971e+04 | 1.102440e+00 | 1.108761e-02 | +| 3 | 1.634907e+04 | 9.942194e-01 | 1.174444e-02 | +| 4 | 1.645914e+04 | 9.961657e-01 | 1.040603e-02 | +| 5 | 1.666816e+04 | 1.008298e+00 | 1.162806e-02 | +| 6 | 1.652520e+04 | 1.013949e+00 | 1.304124e-02 | +| 7 | 1.647338e+04 | 9.847013e-01 | 1.192703e-02 | +| 8 | 1.650048e+04 | 9.890565e-01 | 1.018419e-02 | +| 9 | 1.651135e+04 | 1.010180e+00 | 1.071604e-02 | +| 10 | 1.650611e+04 | 1.001364e+00 | 1.163008e-02 | +| 11 | 1.650539e+04 | 9.971522e-01 | 1.117813e-02 | +| 12 | 1.650588e+04 | 1.020902e+00 | 1.023352e-02 | +| 13 | 1.650545e+04 | 9.994847e-01 | 1.192011e-02 | +| 14 | 1.650546e+04 | 1.002346e+00 | 1.093321e-02 | +| 15 | 1.650546e+04 | 9.943650e-01 | 7.631514e-03 | +--------------------------------------------------------------------------- +``` + +!syntax list /Problem/CriticalitySearch actions=false subsystems=false heading=Available CriticalitySearch Objects + +!syntax parameters /Problem/CriticalitySearch/AddCriticalitySearchAction diff --git a/doc/content/source/criticalitysearch/BoratedWater.md b/doc/content/source/criticalitysearch/BoratedWater.md new file mode 100644 index 000000000..d97afdc9c --- /dev/null +++ b/doc/content/source/criticalitysearch/BoratedWater.md @@ -0,0 +1,33 @@ +# BoratedWater + +## Description + +`BoratedWater` is a [CriticalitySearch](AddCriticalitySearchAction.md) that +performs a criticality search based on the natural boron weight ppm in water. +This object will clear the nuclides in the provided material, and overwrite +with natural water and natural boron, to the specific boron weight ppm. +The density of the material remains unchanged (but it may still be being set by a +coupled thermal-fluid application). + +!alert note +This object does not search in your cross section library to see what nuclides +are present (some libraries do not include O18, for instance). OpenMC will +throw an error if your library is missing a nuclide that `BoratedWater` is +trying to add. To circumvent this, please explicitly note which of the natural +nuclides from hydrogen (H1, H2), boron (B10, B11), and oxygen (O16, O17, O18) +are missing from your material with `absent_nuclides`. Currently, only O18 is supported. + +!alert warning +This object does not currently add S(a,b) tables. To correctly account for thermal scattering in water, the starting material definition in OpenMC must include the correct thermal scattering library. + +## Example Input File Syntax + +The following input will perform a criticality search for each OpenMC solve +by searching via the boron ppm of material 2. + +!listing test/tests/criticality/borated_water/openmc.i + block=Problem + +!syntax parameters /Problem/CriticalitySearch/BoratedWater + +!syntax inputs /Problem/CriticalitySearch/BoratedWater diff --git a/doc/content/source/criticalitysearch/OpenMCMaterialDensity.md b/doc/content/source/criticalitysearch/OpenMCMaterialDensity.md new file mode 100644 index 000000000..cc6d7323f --- /dev/null +++ b/doc/content/source/criticalitysearch/OpenMCMaterialDensity.md @@ -0,0 +1,18 @@ +# OpenMCMaterialDensity + +## Description + +`OpenMCMaterialDensity` is a [CriticalitySearch](AddCriticalitySearchAction.md) that +performs a criticality search based on the total density of a provided material. + +## Example Input File Syntax + +The following input will perform a criticality search for each OpenMC solve +by searching via the density of material ID 1. + +!listing test/tests/criticality/material_density/openmc.i + block=Problem + +!syntax parameters /Problem/CriticalitySearch/OpenMCMaterialDensity + +!syntax inputs /Problem/CriticalitySearch/OpenMCMaterialDensity diff --git a/doc/content/source/postprocessors/NekPressureSurfaceForce.md b/doc/content/source/postprocessors/NekPressureSurfaceForce.md index e36bdde47..ad8830543 100644 --- a/doc/content/source/postprocessors/NekPressureSurfaceForce.md +++ b/doc/content/source/postprocessors/NekPressureSurfaceForce.md @@ -9,8 +9,7 @@ stress tensor over a boundary, multiplied by negative 1 in order to compute the force that the fluid exerts ON the boundary, \begin{equation} -r_i=&\ -\int_{\Gamma}-Pn_i\ d\Gamma\\ -=&\ \int_{\Gamma}Pn_i\ d\Gamma\\ +r_i=-\int_{\Gamma}-Pn_i\ d\Gamma\rightarrow\int_{\Gamma}Pn_i\ d\Gamma \end{equation} where $\Gamma$ is the boundary of the NekRS mesh and diff --git a/doc/content/source/postprocessors/NekViscousSurfaceForce.md b/doc/content/source/postprocessors/NekViscousSurfaceForce.md index 7286819ac..30bceeb1e 100644 --- a/doc/content/source/postprocessors/NekViscousSurfaceForce.md +++ b/doc/content/source/postprocessors/NekViscousSurfaceForce.md @@ -9,7 +9,7 @@ stress tensor over a boundary, multiplied by negative 1 in order to compute the force that the fluid exerts ON the boundary, \begin{equation} -r_i=&\ -\int_{\Gamma}2\mu e_{ij}n_j\ d\Gamma\\ +r_i=-\int_{\Gamma}2\mu e_{ij}n_j\ d\Gamma \end{equation} where $\Gamma$ is the boundary of the NekRS mesh and diff --git a/doc/content/syntax/Problem/CriticalitySearch/index.md b/doc/content/syntax/Problem/CriticalitySearch/index.md new file mode 100644 index 000000000..e69de29bb diff --git a/include/actions/AddCriticalitySearchAction.h b/include/actions/AddCriticalitySearchAction.h new file mode 100644 index 000000000..16aa3159d --- /dev/null +++ b/include/actions/AddCriticalitySearchAction.h @@ -0,0 +1,31 @@ +/********************************************************************/ +/* SOFTWARE COPYRIGHT NOTIFICATION */ +/* Cardinal */ +/* */ +/* (c) 2021 UChicago Argonne, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by UChicago Argonne, LLC */ +/* Under Contract No. DE-AC02-06CH11357 */ +/* With the U. S. Department of Energy */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See LICENSE for full restrictions */ +/********************************************************************/ + +#pragma once + +#include "MooseObjectAction.h" + +class AddCriticalitySearchAction : public MooseObjectAction +{ +public: + static InputParameters validParams(); + + AddCriticalitySearchAction(const InputParameters & parameters); + + virtual void act() override; +}; diff --git a/include/base/OpenMCBase.h b/include/base/OpenMCBase.h index 064294091..903131875 100644 --- a/include/base/OpenMCBase.h +++ b/include/base/OpenMCBase.h @@ -44,6 +44,20 @@ class OpenMCBase */ Real stdev(const double & mean, const double & sum_sq, unsigned int realizations) const; + /** + * A function which computes the mean value of \f$k_{eff}\f$. + * @param[in] estimator type of estimator + * @return the mean value of the k-eigenvalue + */ + Real kMean(const eigenvalue::EigenvalueEnum estimator) const; + + /** + * A function which computes the standard deviation of \f$k_{eff}\f$. + * @param[in] estimator type of estimator + * @return the standard deviation of the k-eigenvalue + */ + Real kStandardDeviation(const eigenvalue::EigenvalueEnum estimator) const; + /// The OpenMCCellAverageProblem required by all objects which inherit from OpenMCBase. OpenMCCellAverageProblem * _openmc_problem; }; diff --git a/include/base/OpenMCProblemBase.h b/include/base/OpenMCProblemBase.h index 27b2f5bed..dd36f89fa 100644 --- a/include/base/OpenMCProblemBase.h +++ b/include/base/OpenMCProblemBase.h @@ -45,6 +45,7 @@ class OpenMCNuclideDensities; class OpenMCDomainFilterEditor; class OpenMCTallyEditor; +class CriticalitySearchBase; /** * Base class for all MOOSE wrappings of OpenMC @@ -65,13 +66,6 @@ class OpenMCProblemBase : public CardinalProblem, public PostprocessorInterface */ std::string subdomainName(const SubdomainID & id) const; - /** - * Print a full error message when catching errors from OpenMC - * @param[in] err OpenMC error code - * @param[in] descriptor descriptive message for error - */ - void catchOpenMCError(const int & err, const std::string descriptor) const; - /** * Whether the score is a reaction rate score * @return whether the tally from OpenMC has units of 1/src @@ -527,4 +521,7 @@ class OpenMCProblemBase : public CardinalProblem, public PostprocessorInterface /// ID used by OpenMC to indicate that a material fill is VOID static constexpr int MATERIAL_VOID{-1}; + + /// Object to use for a criticality search + CriticalitySearchBase * _criticality_search = nullptr; }; diff --git a/include/base/UserErrorChecking.h b/include/base/UserErrorChecking.h index 66cbbb7c8..cf9129377 100644 --- a/include/base/UserErrorChecking.h +++ b/include/base/UserErrorChecking.h @@ -75,3 +75,10 @@ void checkRequiredParam(const InputParameters & p, void checkJointParams(const InputParameters & p, const std::vector & name, const std::string & explanation); + +/** + * Print a full error message when catching errors from OpenMC + * @param[in] err OpenMC error code + * @param[in] descriptor descriptive message for error + */ +void catchOpenMCError(const int & err, const std::string descriptor); diff --git a/include/criticality/BoratedWater.h b/include/criticality/BoratedWater.h new file mode 100644 index 000000000..40479584e --- /dev/null +++ b/include/criticality/BoratedWater.h @@ -0,0 +1,64 @@ +/********************************************************************/ +/* SOFTWARE COPYRIGHT NOTIFICATION */ +/* Cardinal */ +/* */ +/* (c) 2021 UChicago Argonne, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by UChicago Argonne, LLC */ +/* Under Contract No. DE-AC02-06CH11357 */ +/* With the U. S. Department of Energy */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See LICENSE for full restrictions */ +/********************************************************************/ + +#pragma once + +#include "OpenMCMaterialSearch.h" + +/** + * Perform a criticality search based on the boron ppm in water + */ +class BoratedWater : public OpenMCMaterialSearch +{ +public: + static InputParameters validParams(); + + BoratedWater(const InputParameters & parameters); + + virtual void updateOpenMCModel(const Real & input) override; + +protected: + virtual std::string quantity() const override + { + return "material " + std::to_string(_material_id) + " boron"; + } + + virtual std::string units() const override { return "[ppm]"; } + + /// Natural isotopes of hydrogen with their abundances + std::vector> _hydrogen_natural; + + /// Natural isotopes of boron with their abundances + std::vector> _boron_natural; + + /// Natural isotopes of oxygen with their abundances + std::vector> _oxygen_natural; + + /// Molar mass of water + Real _M_H2O; + + /// Molar mass of boron + Real _M_B; + +private: + /** Adjust the natural abundances used for the criticality search if nuclides + * are missing from the cross section library. + * @param[in] allowable possible natural isotopes that will exist in the borated water + */ + void applyAbsentNuclides(const std::vector & allowable); +}; diff --git a/include/criticality/CriticalitySearchBase.h b/include/criticality/CriticalitySearchBase.h new file mode 100644 index 000000000..8205cca11 --- /dev/null +++ b/include/criticality/CriticalitySearchBase.h @@ -0,0 +1,78 @@ +/********************************************************************/ +/* SOFTWARE COPYRIGHT NOTIFICATION */ +/* Cardinal */ +/* */ +/* (c) 2021 UChicago Argonne, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by UChicago Argonne, LLC */ +/* Under Contract No. DE-AC02-06CH11357 */ +/* With the U. S. Department of Energy */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See LICENSE for full restrictions */ +/********************************************************************/ + +#pragma once + +#include "MooseObject.h" +#include "OpenMCBase.h" +#include "CardinalEnums.h" + +/** + * Applies a criticality search for each OpenMC solve. The critical configuration + * is then used for the multiphysics coupling. + */ +class CriticalitySearchBase : public MooseObject, public OpenMCBase +{ +public: + static InputParameters validParams(); + + CriticalitySearchBase(const InputParameters & parameters); + + /** + * Modify the OpenMC model + * @param[in] input value to apply to the new model; interpretation depends on derived classes + */ + virtual void updateOpenMCModel(const Real & input) = 0; + + /// Use Brent's method to search for criticality + virtual void searchForCriticality(); + +protected: + /// The quantity being varied in the search for criticality, for console prints + virtual std::string quantity() const = 0; + + /// Assumed units in the input quantities + virtual std::string units() const = 0; + + /// Target k_eff for search + const Real & _keff_target; + + /// Maximum range of value to explore + const Real & _maximum; + + /// Minimum range of value to explore + const Real & _minimum; + + /// Absolute tolerance for finding a critical configuration + const Real & _tolerance; + + /// Estimator to use for k + const eigenvalue::EigenvalueEnum _estimator; + + /// Values used in search + std::vector _inputs; + + /// Values obtained in search + std::vector _k_values; + + /// Standard deviation values obtained in search + std::vector _k_std_dev_values; + + /// Postprocessor that holds the result of the criticality search + const std::string _pp_name = "critical_value"; +}; diff --git a/include/criticality/OpenMCMaterialDensity.h b/include/criticality/OpenMCMaterialDensity.h new file mode 100644 index 000000000..68b220fff --- /dev/null +++ b/include/criticality/OpenMCMaterialDensity.h @@ -0,0 +1,42 @@ +/********************************************************************/ +/* SOFTWARE COPYRIGHT NOTIFICATION */ +/* Cardinal */ +/* */ +/* (c) 2021 UChicago Argonne, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by UChicago Argonne, LLC */ +/* Under Contract No. DE-AC02-06CH11357 */ +/* With the U. S. Department of Energy */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See LICENSE for full restrictions */ +/********************************************************************/ + +#pragma once + +#include "OpenMCMaterialSearch.h" + +/** + * Perform a criticality search based on a material total density + */ +class OpenMCMaterialDensity : public OpenMCMaterialSearch +{ +public: + static InputParameters validParams(); + + OpenMCMaterialDensity(const InputParameters & parameters); + + virtual void updateOpenMCModel(const Real & input) override; + +protected: + virtual std::string quantity() const override + { + return "material " + std::to_string(_material_id) + " density"; + } + + virtual std::string units() const override { return "[kg/m3]"; } +}; diff --git a/include/criticality/OpenMCMaterialSearch.h b/include/criticality/OpenMCMaterialSearch.h new file mode 100644 index 000000000..70dc80274 --- /dev/null +++ b/include/criticality/OpenMCMaterialSearch.h @@ -0,0 +1,39 @@ +/********************************************************************/ +/* SOFTWARE COPYRIGHT NOTIFICATION */ +/* Cardinal */ +/* */ +/* (c) 2021 UChicago Argonne, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by UChicago Argonne, LLC */ +/* Under Contract No. DE-AC02-06CH11357 */ +/* With the U. S. Department of Energy */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See LICENSE for full restrictions */ +/********************************************************************/ + +#pragma once + +#include "CriticalitySearchBase.h" + +/** + * Perform a criticality search based on a material + */ +class OpenMCMaterialSearch : public CriticalitySearchBase +{ +public: + static InputParameters validParams(); + + OpenMCMaterialSearch(const InputParameters & parameters); + +protected: + /// Material to be modified + const int32_t & _material_id; + + /// Material index corresponding to the ID + int32_t _material_index; +}; diff --git a/include/postprocessors/KEigenvalue.h b/include/postprocessors/KEigenvalue.h index 254572ccb..9d0517907 100644 --- a/include/postprocessors/KEigenvalue.h +++ b/include/postprocessors/KEigenvalue.h @@ -39,18 +39,6 @@ class KEigenvalue : public GeneralPostprocessor, public OpenMCBase virtual Real getValue() const override; protected: - /** - * A function which computes the mean value of k_{eff}. - * @return the mean value of the k-eigenvalue - */ - Real kMean() const; - - /** - * A function which computes the standard deviation of k_{eff}. - * @return the standard deviation of the k-eigenvalue - */ - Real KStandardDeviation() const; - /** * A function which computes the relative error of k_{eff}. * @return the relative error of the k-eigenvalue diff --git a/scripts/get-dependencies.sh b/scripts/get-dependencies.sh index 3666f8915..4ba909845 100755 --- a/scripts/get-dependencies.sh +++ b/scripts/get-dependencies.sh @@ -6,6 +6,7 @@ set -ex git submodule update --init --progress contrib/moose git submodule update --init --progress contrib/nekRS git submodule update --init --recursive --progress contrib/openmc +git submodule update --init --progress contrib/nuclear_data git submodule update --init --progress contrib/DAGMC git submodule update --init --progress contrib/moab git submodule update --init --progress test/tests/nek_ci diff --git a/src/actions/AddCriticalitySearchAction.C b/src/actions/AddCriticalitySearchAction.C new file mode 100644 index 000000000..b0a287748 --- /dev/null +++ b/src/actions/AddCriticalitySearchAction.C @@ -0,0 +1,54 @@ +/********************************************************************/ +/* SOFTWARE COPYRIGHT NOTIFICATION */ +/* Cardinal */ +/* */ +/* (c) 2021 UChicago Argonne, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by UChicago Argonne, LLC */ +/* Under Contract No. DE-AC02-06CH11357 */ +/* With the U. S. Department of Energy */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See LICENSE for full restrictions */ +/********************************************************************/ + +#ifdef ENABLE_OPENMC_COUPLING + +#include "AddCriticalitySearchAction.h" +#include "OpenMCCellAverageProblem.h" +#include "CriticalitySearchBase.h" + +registerMooseAction("CardinalApp", AddCriticalitySearchAction, "add_criticality_search"); + +InputParameters +AddCriticalitySearchAction::validParams() +{ + auto params = MooseObjectAction::validParams(); + params.addClassDescription("Adds a criticality search for OpenMC"); + return params; +} + +AddCriticalitySearchAction::AddCriticalitySearchAction(const InputParameters & parameters) + : MooseObjectAction(parameters) +{ +} + +void +AddCriticalitySearchAction::act() +{ + if (_current_task == "add_criticality_search") + { + auto openmc_problem = dynamic_cast(_problem.get()); + + if (!openmc_problem) + mooseError("The [CriticalitySearch] block can only be used with wrapped OpenMC cases! " + "You need to change the [Problem] block to 'OpenMCCellAverageProblem'."); + + openmc_problem->addObject(_type, _name, _moose_object_pars, false)[0]; + } +} +#endif diff --git a/src/actions/AddFieldTransferAction.C b/src/actions/AddFieldTransferAction.C index 5d1035e7c..398665538 100644 --- a/src/actions/AddFieldTransferAction.C +++ b/src/actions/AddFieldTransferAction.C @@ -48,13 +48,8 @@ AddFieldTransferAction::act() mooseError("The [FieldTransfers] block can only be used with wrapped Nek cases! " "You need to change the [Problem] block to 'NekRSProblem'."); - if (_type == "NekFieldVariable" || _type == "NekVolumetricSource" || - _type == "NekBoundaryFlux" || _type == "NekMeshDeformation") - { - _moose_object_pars.set("_nek_problem") = nek_problem; - auto transfer = - nek_problem->addObject(_type, _name, _moose_object_pars, false)[0]; - } + _moose_object_pars.set("_nek_problem") = nek_problem; + nek_problem->addObject(_type, _name, _moose_object_pars, false)[0]; } } #endif diff --git a/src/actions/AddScalarTransferAction.C b/src/actions/AddScalarTransferAction.C index d8c57cbd5..fe7f9f2e3 100644 --- a/src/actions/AddScalarTransferAction.C +++ b/src/actions/AddScalarTransferAction.C @@ -49,12 +49,8 @@ AddScalarTransferAction::act() mooseError("The [ScalarTransfers] block can only be used with wrapped Nek cases! " "You need to change the [Problem] block to 'NekRSProblem'."); - if (_type == "NekScalarValue" || _type == "NekPostprocessorValue") - { - _moose_object_pars.set("_nek_problem") = nek_problem; - auto transfer = - nek_problem->addObject(_type, _name, _moose_object_pars, false)[0]; - } + _moose_object_pars.set("_nek_problem") = nek_problem; + nek_problem->addObject(_type, _name, _moose_object_pars, false)[0]; } } #endif diff --git a/src/base/CardinalApp.C b/src/base/CardinalApp.C index b8bf44ca8..dd33353cd 100644 --- a/src/base/CardinalApp.C +++ b/src/base/CardinalApp.C @@ -219,6 +219,12 @@ CardinalApp::associateSyntaxInner(Syntax & syntax, ActionFactory & /* action_fac // Can only add external auxvars after the tallies have been added. addTaskDependency("add_external_aux_variables", "add_tallies"); + // Add the [Problem/CriticalitySearch] block + registerSyntax("AddCriticalitySearchAction", "Problem/CriticalitySearch"); + registerMooseObjectTask("add_criticality_search", CriticalitySearch, false); + registerTask("add_criticality_search", false /* is required */); + addTaskDependency("add_criticality_search", "init_problem"); + // Register a modify outputs task to enable variable hiding in the MGXS action. registerTask("modify_outputs", true /* is required */); addTaskDependency("modify_outputs", "common_output"); diff --git a/src/base/OpenMCBase.C b/src/base/OpenMCBase.C index 7643af823..bb083be89 100644 --- a/src/base/OpenMCBase.C +++ b/src/base/OpenMCBase.C @@ -19,7 +19,8 @@ #ifdef ENABLE_OPENMC_COUPLING #include "OpenMCBase.h" - +#include "openmc/eigenvalue.h" +#include "openmc/math_functions.h" #include "libmesh/elem.h" InputParameters @@ -46,4 +47,82 @@ OpenMCBase::stdev(const double & mean, const double & sum_sq, unsigned int reali : 0.0; } +Real +OpenMCBase::kMean(const eigenvalue::EigenvalueEnum estimator) const +{ + int n = openmc::simulation::n_realizations; + const auto & gt = openmc::simulation::global_tallies; + + switch (estimator) + { + case eigenvalue::collision: + return gt(openmc::GlobalTally::K_COLLISION, openmc::TallyResult::SUM) / n; + case eigenvalue::absorption: + return gt(openmc::GlobalTally::K_ABSORPTION, openmc::TallyResult::SUM) / n; + case eigenvalue::tracklength: + return gt(openmc::GlobalTally::K_TRACKLENGTH, openmc::TallyResult::SUM) / n; + case eigenvalue::combined: + { + if (n <= 3) + mooseError("Cannot compute combined k-effective estimate with fewer than 4 realizations!\n" + "Please change the estimator type to either 'collision', 'tracklength', or " + "'absorption'."); + + double k_eff[2]; + openmc::openmc_get_keff(k_eff); + return k_eff[0]; + } + default: + mooseError("Internal error: Unhandled EigenvalueEnum!"); + } +} + +Real +OpenMCBase::kStandardDeviation(const eigenvalue::EigenvalueEnum estimator) const +{ + const auto & gt = openmc::simulation::global_tallies; + int n = openmc::simulation::n_realizations; + + double t_n1 = 1.0; + if (openmc::settings::confidence_intervals) + { + double alpha = 1.0 - openmc::CONFIDENCE_LEVEL; + t_n1 = openmc::t_percentile(1.0 - alpha / 2.0, n - 1); + } + + switch (estimator) + { + case eigenvalue::collision: + { + double mean = gt(openmc::GlobalTally::K_COLLISION, openmc::TallyResult::SUM) / n; + double sum_sq = gt(openmc::GlobalTally::K_COLLISION, openmc::TallyResult::SUM_SQ); + return t_n1 * stdev(mean, sum_sq, openmc::simulation::n_realizations); + } + case eigenvalue::absorption: + { + double mean = gt(openmc::GlobalTally::K_ABSORPTION, openmc::TallyResult::SUM) / n; + double sum_sq = gt(openmc::GlobalTally::K_ABSORPTION, openmc::TallyResult::SUM_SQ); + return t_n1 * stdev(mean, sum_sq, openmc::simulation::n_realizations); + } + case eigenvalue::tracklength: + { + double mean = gt(openmc::GlobalTally::K_TRACKLENGTH, openmc::TallyResult::SUM) / n; + double sum_sq = gt(openmc::GlobalTally::K_TRACKLENGTH, openmc::TallyResult::SUM_SQ); + return t_n1 * stdev(mean, sum_sq, openmc::simulation::n_realizations); + } + case eigenvalue::combined: + { + if (n <= 3) + mooseError("Cannot compute combined k-effective standard deviation with fewer than 4 " + "realizations!"); + + double k_eff[2]; + openmc::openmc_get_keff(k_eff); + return k_eff[1]; + } + default: + mooseError("Internal error: Unhandled StandardDeviationEnum!"); + } +} + #endif diff --git a/src/base/OpenMCCellAverageProblem.C b/src/base/OpenMCCellAverageProblem.C index 82755d9f4..20a44643c 100644 --- a/src/base/OpenMCCellAverageProblem.C +++ b/src/base/OpenMCCellAverageProblem.C @@ -1485,7 +1485,7 @@ OpenMCCellAverageProblem::getMaterialFills() if (!hasDensityFeedback(cell_info)) { // TODO: this check should be extended for non-fluid cells which may contain - // lattices or universes + // lattices or universes if (is_material_cell) other_materials.insert(material_index); continue; diff --git a/src/base/OpenMCProblemBase.C b/src/base/OpenMCProblemBase.C index e84d1579b..e78999528 100644 --- a/src/base/OpenMCProblemBase.C +++ b/src/base/OpenMCProblemBase.C @@ -27,6 +27,7 @@ #include "OpenMCNuclideDensities.h" #include "OpenMCDomainFilterEditor.h" #include "OpenMCTallyEditor.h" +#include "CriticalitySearchBase.h" #include "openmc/random_lcg.h" @@ -232,14 +233,6 @@ OpenMCProblemBase::OpenMCProblemBase(const InputParameters & params) OpenMCProblemBase::~OpenMCProblemBase() { openmc_finalize(); } -void -OpenMCProblemBase::catchOpenMCError(const int & err, const std::string descriptor) const -{ - if (err) - mooseError("In attempting to ", descriptor, ", OpenMC reported:\n\n", - std::string(openmc_err_msg)); -} - void OpenMCProblemBase::fillElementalAuxVariable(const unsigned int & var_num, const std::vector & elem_ids, @@ -346,7 +339,7 @@ OpenMCProblemBase::externalSolve() return; } - _console << " Running OpenMC with " << nParticles() << " particles per batch..." << std::endl; + _console << "Running OpenMC with " << nParticles() << " particles per batch..." << std::endl; // apply a new starting fission source if (_reuse_source && !firstSolve()) @@ -365,9 +358,15 @@ OpenMCProblemBase::externalSolve() openmc_set_seed(_initial_seed); } - int err = openmc_run(); - if (err) - mooseError(openmc_err_msg); + int err; + if (_criticality_search) + _criticality_search->searchForCriticality(); + else + { + err = openmc_run(); + if (err) + mooseError(openmc_err_msg); + } _total_n_particles += nParticles(); @@ -385,7 +384,7 @@ OpenMCProblemBase::externalSolve() void OpenMCProblemBase::initialSetup() { - ExternalProblem::initialSetup(); + CardinalProblem::initialSetup(); // Initialize the IFP parameters tally. if (_calc_kinetics_params) @@ -395,6 +394,17 @@ OpenMCProblemBase::initialSetup() _ifp_tally->set_scores({"ifp-time-numerator", "ifp-beta-numerator", "ifp-denominator"}); _ifp_tally->estimator_ = openmc::TallyEstimator::COLLISION; } + + // Find a criticality search object + TheWarehouse::Query query = theWarehouse().query().condition("CriticalitySearch"); + std::vector objs; + query.queryInto(objs); + + if (objs.size() > 1) + mooseError("Cannot have more than one CriticalitySearch object"); + + if (objs.size()) + _criticality_search = objs[0]; } void @@ -895,6 +905,10 @@ void OpenMCProblemBase::executeFilterEditors() { executeControls(EXEC_FILTER_EDITORS); + + if (!_filter_editor_uos.size()) + return; + _console << "Executing filter editors..." << std::endl; for (const auto & fe : _filter_editor_uos) fe->execute(); @@ -904,6 +918,10 @@ void OpenMCProblemBase::executeTallyEditors() { executeControls(EXEC_TALLY_EDITORS); + + if (!_tally_editor_uos.size()) + return; + _console << "Executing tally editors..." << std::endl; for (const auto & te : _tally_editor_uos) te->execute(); diff --git a/src/base/UserErrorChecking.C b/src/base/UserErrorChecking.C index 209fb4f0f..3f25e281a 100644 --- a/src/base/UserErrorChecking.C +++ b/src/base/UserErrorChecking.C @@ -89,3 +89,13 @@ checkJointParams(const InputParameters & p, "be specified or ALL omitted; you have only provided a subset of parameters!"); } } + +#ifdef ENABLE_OPENMC_COUPLING +void +catchOpenMCError(const int & err, const std::string descriptor) +{ + if (err) + mooseError( + "In attempting to ", descriptor, ", OpenMC reported:\n\n", std::string(openmc_err_msg)); +} +#endif diff --git a/src/criticality/BoratedWater.C b/src/criticality/BoratedWater.C new file mode 100644 index 000000000..ec024ab98 --- /dev/null +++ b/src/criticality/BoratedWater.C @@ -0,0 +1,219 @@ +/********************************************************************/ +/* SOFTWARE COPYRIGHT NOTIFICATION */ +/* Cardinal */ +/* */ +/* (c) 2021 UChicago Argonne, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by UChicago Argonne, LLC */ +/* Under Contract No. DE-AC02-06CH11357 */ +/* With the U. S. Department of Energy */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See LICENSE for full restrictions */ +/********************************************************************/ + +#ifdef ENABLE_OPENMC_COUPLING + +#define NUCLEAR_DATA_IMPLEMENTATION + +#include "BoratedWater.h" +#include "UserErrorChecking.h" +#include "NuclearData.h" +#include "openmc/capi.h" +#include "openmc/constants.h" +#include "openmc/nuclide.h" + +registerMooseObject("CardinalApp", BoratedWater); + +InputParameters +BoratedWater::validParams() +{ + auto params = OpenMCMaterialSearch::validParams(); + params.addParam>( + "absent_nuclides", + "Natural nuclides of hydrogen, oxygen, or boron which are missing from your cross section " + "library; some cross section libraries do not have entries for O17 and O18. If your library " + "does not have these nuclides you will get an error from this object trying to add them. For " + "these missing nuclides, specify them here and their abundance will be applied to the main " + "isotope of each element. Currently, only O18 is supported."); + params.addClassDescription( + "Searches for criticality using natural boron ppm in water in units of weight ppm"); + return params; +} + +BoratedWater::BoratedWater(const InputParameters & parameters) : OpenMCMaterialSearch(parameters) +{ + // apply additional checks on the minimum and maximum + if (_minimum < 0) + paramError("minimum", + "The 'minimum' boron ppm (" + std::to_string(_minimum) + ") must be positive!"); + if (_maximum > 50000) + paramError("maximum", + "The borated water composition is computed using a dilute species approximation. " + "Results will not be accurate if the boron species is no longer dilute, which we " + "take as 5\% weight concentration or higher. Please decrease 'maximum' (" + + std::to_string(_maximum) + ") to an upper limit which is in the dilute regime."); + + // We take the provided material, retaining its density, but then overwriting + // any nuclides there to be H2O + boron weight ppm and add the S(a,b) tables + // for hydrogen. Therefore, we perform a check here if any nuclides in the + // material are *not* those we are going to add back when we wipe the material, + // to notify the user they probably provided the wrong material ID or they need + // to use a more general material modification object that allows full control + // over additional nuclides (e.g., if their water includes corrosion products + // that they do want to be there). + // TODO: add the S(a,b) tables + + const int * nuclides; + const double * densities; + int n; + int err = openmc_material_get_densities(_material_index, &nuclides, &densities, &n); + catchOpenMCError(err, "get nuclide densities from material " + std::to_string(_material_id)); + + // get all the natural isotopes of H, O, B + _hydrogen_natural = NuclearData::Nuclide::getAbundances("H"); + _boron_natural = NuclearData::Nuclide::getAbundances("B"); + _oxygen_natural = NuclearData::Nuclide::getAbundances("O"); + + std::vector allowable; + for (const auto & i : _hydrogen_natural) + allowable.push_back(i.first); + for (const auto & i : _boron_natural) + allowable.push_back(i.first); + for (const auto & i : _oxygen_natural) + allowable.push_back(i.first); + + applyAbsentNuclides(allowable); + + // check if any nuclides already defined on the material do not intersect with + // natural isotopes of H, O, B + std::string full_names = ""; + for (int i = 0; i < n; ++i) + { + auto idx = nuclides[i]; + std::string name = openmc::data::nuclides[idx]->name_; + if (std::find(allowable.begin(), allowable.end(), name) == allowable.end()) + full_names += name + " "; + } + + if (!full_names.empty()) + { + std::ostringstream msg; + msg << "The criticality search will clear out all nuclides in material " + << std::to_string(_material_id) + << " and replace them with the naturally-abundant nuclides in hydrogen, oxygen, and boron. " + "Any other nuclides which existed in the material will be deleted." + << std::endl; + msg << "\nThe material you provided contains nuclides which are not the natural isotopes of H, " + "B, and O: " + << full_names.substr(0, full_names.length() - 1) << ". These will be deleted from material " + << std::to_string(_material_id) << " when the boron concentration is changed." << std::endl; + msg << "\nFor general criticality searches based on material composition, please contact the " + "Cardinal developer team."; + mooseWarning(msg.str()); + } + + // Compute the molar mass of pure water + _M_H2O = 0.0; + for (const auto & h : _hydrogen_natural) + _M_H2O += 2.0 * h.second * NuclearData::Nuclide::getAtomicMass(h.first); + for (const auto & o : _oxygen_natural) + _M_H2O += 1.0 * o.second * NuclearData::Nuclide::getAtomicMass(o.first); + + // Compute the molar mass of pure boron + _M_B = 0.0; + for (const auto & b : _boron_natural) + _M_B += 1.0 * b.second * NuclearData::Nuclide::getAtomicMass(b.first); +} + +void +BoratedWater::updateOpenMCModel(const Real & ppm) +{ + _console << "Searching for boron = " << ppm << " [ppm] ..." << std::endl; + + // A coupled thermal-fluid app may set the material density just before we enter + // this routine; we will preserve that density which may be set and we interpret + // it to be the SOLUTION density to be fully consistent with energy conservation. + + // Compute the number fractions of each element + Real frac_H2O = (1 - ppm * 1e-6) / _M_H2O; + Real frac_H = 2 * frac_H2O; + Real frac_O = frac_H2O; + Real frac_B = ppm * 1e-6 / _M_B; + + double rho; + int err = openmc_material_get_density(_material_index, &rho); + catchOpenMCError(err, "get density for material " + std::to_string(_material_id)); + + std::vector names; + std::vector densities; + for (const auto & h : _hydrogen_natural) + { + names.push_back(h.first); + densities.push_back(h.second * frac_H * rho * openmc::N_AVOGADRO); + } + for (const auto & o : _oxygen_natural) + { + names.push_back(o.first); + densities.push_back(o.second * frac_O * rho * openmc::N_AVOGADRO); + } + for (const auto & b : _boron_natural) + { + names.push_back(b.first); + densities.push_back(b.second * frac_B * rho * openmc::N_AVOGADRO); + } + + openmc::model::materials[_material_index]->set_densities(names, densities /* atom/b-cm */); +} + +void +BoratedWater::applyAbsentNuclides(const std::vector & allowable) +{ + if (isParamValid("absent_nuclides")) + { + const auto & absent = getParam>("absent_nuclides"); + + for (const auto & a : absent) + { + // only missing nuclides for H, O, B are meaningful + if (std::find(allowable.begin(), allowable.end(), a) == allowable.end()) + paramWarning("absent_nuclides", + "Only absent isotopes of hydrogen, oxygen, or boron will be used to adjust " + "natural abundances. The entry '" + + a + "' will be unused."); + + // adjust the natural abundances if nuclides are missing from the cross section library + // TODO: implement in a general fashion, this only works for O18 which is also the only + // absent nuclide we expect in practice + if (a == "O18") + { + // find which entry in _oxygen_natural has O16 and O18; this allows these isotopes + // to appear in any order + unsigned int idx16; + unsigned int idx18; + unsigned int idx = 0; + for (const auto & o : _oxygen_natural) + { + if (o.first == "O16") + idx16 = idx; + if (o.first == "O18") + idx18 = idx; + idx++; + } + + _oxygen_natural[idx16].second += _oxygen_natural[idx18].second; + _oxygen_natural.erase(_oxygen_natural.begin() + idx18); + } + else + paramError( + "absent_nuclides", + "Cardinal currently only assumes that O18 may be missing from your cross section " + "library; please contact the Cardinal developer team to generalize this capability"); + } + } +} +#endif diff --git a/src/criticality/CriticalitySearchBase.C b/src/criticality/CriticalitySearchBase.C new file mode 100644 index 000000000..94d9f0bd0 --- /dev/null +++ b/src/criticality/CriticalitySearchBase.C @@ -0,0 +1,132 @@ +/********************************************************************/ +/* SOFTWARE COPYRIGHT NOTIFICATION */ +/* Cardinal */ +/* */ +/* (c) 2021 UChicago Argonne, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by UChicago Argonne, LLC */ +/* Under Contract No. DE-AC02-06CH11357 */ +/* With the U. S. Department of Energy */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See LICENSE for full restrictions */ +/********************************************************************/ + +#ifdef ENABLE_OPENMC_COUPLING + +#include "CriticalitySearchBase.h" +#include "UserErrorChecking.h" +#include "VariadicTable.h" +#include "BrentsMethod.h" + +InputParameters +CriticalitySearchBase::validParams() +{ + auto params = MooseObject::validParams(); + params += OpenMCBase::validParams(); + params.addRangeCheckParam( + "keff_target", + 1.0, + "keff_target > 0", + "Target value of k_eff for the criticality search."); + params.addRequiredParam( + "minimum", + "Minimum for values to search over; the root must occur at a value greater than the minimum"); + params.addRequiredParam( + "maximum", + "Maximum for values to search over; the root must occur at a value smaller than the maximum"); + params.addRangeCheckedParam( + "tolerance", + 1e-3, + "tolerance > 0", + "Absolute tolerance to converge multiplication factor; be aware that if too few particles " + "are used, statistical noise may require many criticality calculations to converge."); + params.addParam( + "estimator", getEigenvalueEnum(), "Type of eigenvalue estimator to use"); + params.addClassDescription( + "Base class for defining parameters used in a criticality search in OpenMC."); + params.registerBase("CriticalitySearch"); + params.registerSystemAttributeName("CriticalitySearch"); + params.addPrivateParam("_openmc_problem"); + return params; +} + +CriticalitySearchBase::CriticalitySearchBase(const InputParameters & parameters) + : MooseObject(parameters), + OpenMCBase(this, parameters), + _keff_target(getParam("keff_target")), + _maximum(getParam("maximum")), + _minimum(getParam("minimum")), + _tolerance(getParam("tolerance")), + _estimator(getParam("estimator").getEnum()) +{ + if (_minimum >= _maximum) + paramError("minimum", + "The 'minimum' value (" + std::to_string(_minimum) + + ") must be less than the 'maximum' value (" + std::to_string(_maximum) + ")."); + + auto pp_params = _factory.getValidParams("Receiver"); + _openmc_problem->addPostprocessor("Receiver", _pp_name, pp_params); +} + +void +CriticalitySearchBase::searchForCriticality() +{ + _console << "Running criticality search in OpenMC for " << quantity() << " in range " + << std::to_string(_minimum) << " - " << std::to_string(_maximum) << " " << units() << " " + << std::endl; + + VariadicTable vt( + {"Iteration", quantity() + " " + units(), " k (mean) ", " k (std dev) "}); + vt.setColumnFormat({VariadicTableColumnFormat::AUTO, + VariadicTableColumnFormat::SCIENTIFIC, + VariadicTableColumnFormat::SCIENTIFIC, + VariadicTableColumnFormat::SCIENTIFIC}); + + std::function func; + func = [&](Real x) + { + // update the OpenMC model with a new parameter + updateOpenMCModel(x); + _inputs.push_back(x); + + // re-run the model + int err = openmc_run(); + if (err) + mooseError(openmc_err_msg); + + // fetch k and return the residual, k-1 + Real k = kMean(_estimator); + Real k_std_dev = kStandardDeviation(_estimator); + _k_values.push_back(k); + _k_std_dev_values.push_back(k_std_dev); + + vt.addRow(_k_values.size() - 1, x, k, k_std_dev); + vt.print(_console); + + if (_tolerance < 3 * k_std_dev) + mooseDoOnce(mooseWarning( + "The 'tolerance' for the criticality search (" + std::to_string(_tolerance) + + ") is smaller than 3-sigma standard deviation in k (" + std::to_string(3 * k_std_dev) + + "); you may have to run a lot of criticality search points to converge to this " + "tolerance. You may want to loosen 'tolerance' or increase the number of particles.")); + + return k - _keff_target; + }; + + BrentsMethod::root(func, _minimum, _maximum, _tolerance); + + // check if the method converged + if (abs(kMean(_estimator) - _keff_target) >= _tolerance) + mooseError("Failed to converge criticality search! This may happen if your tolerance is too " + "tight given the statistical error in the computation of k."); + + // fill the converged value into a postprocessor + _openmc_problem->setPostprocessorValueByName(_pp_name, _inputs.back()); +} + +#endif diff --git a/src/criticality/OpenMCMaterialDensity.C b/src/criticality/OpenMCMaterialDensity.C new file mode 100644 index 000000000..d74a7b04d --- /dev/null +++ b/src/criticality/OpenMCMaterialDensity.C @@ -0,0 +1,59 @@ +/********************************************************************/ +/* SOFTWARE COPYRIGHT NOTIFICATION */ +/* Cardinal */ +/* */ +/* (c) 2021 UChicago Argonne, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by UChicago Argonne, LLC */ +/* Under Contract No. DE-AC02-06CH11357 */ +/* With the U. S. Department of Energy */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See LICENSE for full restrictions */ +/********************************************************************/ + +#ifdef ENABLE_OPENMC_COUPLING + +#include "OpenMCMaterialDensity.h" +#include "openmc/capi.h" + +registerMooseObject("CardinalApp", OpenMCMaterialDensity); + +InputParameters +OpenMCMaterialDensity::validParams() +{ + auto params = OpenMCMaterialSearch::validParams(); + params.addClassDescription("Searches for criticality using material density in units of kg/m3"); + return params; +} + +OpenMCMaterialDensity::OpenMCMaterialDensity(const InputParameters & parameters) + : OpenMCMaterialSearch(parameters) +{ + // a material in OpenMC must always have some nuclides in it or macroscopic data, + // so we don't need to have any checks on whether the material is void + + // apply additional checks on the minimum and maximum values - both must be positive. + // we don't need to check for negative 'maximum' because we already require maximum > minimum + // and if we enforce non-negative minimum this will require maximum to also be non-negative. + if (_minimum < 0.0) + paramError("minimum", + "The 'minimum' density (" + std::to_string(_minimum) + ") must be positive!"); +} + +void +OpenMCMaterialDensity::updateOpenMCModel(const Real & density) +{ + _console << "Searching for density = " << density << " [kg/m3] ..." << std::endl; + + const char * units = "g/cc"; + int err = openmc_material_set_density( + _material_index, density * _openmc_problem->densityConversionFactor(), units); + catchOpenMCError(err, "set material density to " + std::to_string(density)); +} + +#endif diff --git a/src/criticality/OpenMCMaterialSearch.C b/src/criticality/OpenMCMaterialSearch.C new file mode 100644 index 000000000..13e0c4a1d --- /dev/null +++ b/src/criticality/OpenMCMaterialSearch.C @@ -0,0 +1,42 @@ +/********************************************************************/ +/* SOFTWARE COPYRIGHT NOTIFICATION */ +/* Cardinal */ +/* */ +/* (c) 2021 UChicago Argonne, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by UChicago Argonne, LLC */ +/* Under Contract No. DE-AC02-06CH11357 */ +/* With the U. S. Department of Energy */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See LICENSE for full restrictions */ +/********************************************************************/ + +#ifdef ENABLE_OPENMC_COUPLING + +#include "OpenMCMaterialSearch.h" +#include "UserErrorChecking.h" +#include "openmc/capi.h" + +InputParameters +OpenMCMaterialSearch::validParams() +{ + auto params = CriticalitySearchBase::validParams(); + params.addRequiredParam("material_id", "Material ID to modify"); + params.addClassDescription( + "Base class for criticality searches using the properties of a material"); + return params; +} + +OpenMCMaterialSearch::OpenMCMaterialSearch(const InputParameters & parameters) + : CriticalitySearchBase(parameters), _material_id(getParam("material_id")) +{ + int err = openmc_get_material_index(_material_id, &_material_index); + catchOpenMCError(err, "get index for material with ID " + std::to_string(_material_id)); +} + +#endif diff --git a/src/postprocessors/KEigenvalue.C b/src/postprocessors/KEigenvalue.C index 9cf27d28f..288c7bfeb 100644 --- a/src/postprocessors/KEigenvalue.C +++ b/src/postprocessors/KEigenvalue.C @@ -19,9 +19,6 @@ #ifdef ENABLE_OPENMC_COUPLING #include "KEigenvalue.h" -#include "openmc/eigenvalue.h" -#include "openmc/math_functions.h" -#include "openmc/constants.h" registerMooseObject("CardinalApp", KEigenvalue); @@ -58,10 +55,10 @@ KEigenvalue::getValue() const switch (_output) { case statistics::OutputEnum::Mean: - return kMean(); + return kMean(_type); case statistics::OutputEnum::StDev: - return KStandardDeviation(); + return kStandardDeviation(_type); case statistics::OutputEnum::RelError: return kRelativeError(); @@ -72,88 +69,11 @@ KEigenvalue::getValue() const } } -Real -KEigenvalue::kMean() const -{ - int n = openmc::simulation::n_realizations; - const auto & gt = openmc::simulation::global_tallies; - - switch (_type) - { - case eigenvalue::collision: - return gt(openmc::GlobalTally::K_COLLISION, openmc::TallyResult::SUM) / n; - case eigenvalue::absorption: - return gt(openmc::GlobalTally::K_ABSORPTION, openmc::TallyResult::SUM) / n; - case eigenvalue::tracklength: - return gt(openmc::GlobalTally::K_TRACKLENGTH, openmc::TallyResult::SUM) / n; - case eigenvalue::combined: - { - if (n <= 3) - mooseError("Cannot compute combined k-effective estimate with fewer than 4 realizations!\n" - "Please change the 'value_type' to either 'collision', 'tracklength', or 'absorption'."); - - double k_eff[2]; - openmc::openmc_get_keff(k_eff); - return k_eff[0]; - } - default: - mooseError("Internal error: Unhandled EigenvalueEnum in KEigenvalue!"); - } -} - -Real -KEigenvalue::KStandardDeviation() const -{ - const auto & gt = openmc::simulation::global_tallies; - int n = openmc::simulation::n_realizations; - - double t_n1 = 1.0; - if (openmc::settings::confidence_intervals) - { - double alpha = 1.0 - openmc::CONFIDENCE_LEVEL; - t_n1 = openmc::t_percentile(1.0 - alpha / 2.0, n - 1); - } - - switch (_type) - { - case eigenvalue::collision: - { - double mean = gt(openmc::GlobalTally::K_COLLISION, openmc::TallyResult::SUM) / n; - double sum_sq = gt(openmc::GlobalTally::K_COLLISION, openmc::TallyResult::SUM_SQ); - return t_n1 * stdev(mean, sum_sq, openmc::simulation::n_realizations); - } - case eigenvalue::absorption: - { - double mean = gt(openmc::GlobalTally::K_ABSORPTION, openmc::TallyResult::SUM) / n; - double sum_sq = gt(openmc::GlobalTally::K_ABSORPTION, openmc::TallyResult::SUM_SQ); - return t_n1 * stdev(mean, sum_sq, openmc::simulation::n_realizations); - } - case eigenvalue::tracklength: - { - double mean = gt(openmc::GlobalTally::K_TRACKLENGTH, openmc::TallyResult::SUM) / n; - double sum_sq = gt(openmc::GlobalTally::K_TRACKLENGTH, openmc::TallyResult::SUM_SQ); - return t_n1 * stdev(mean, sum_sq, openmc::simulation::n_realizations); - } - case eigenvalue::combined: - { - if (n <= 3) - mooseError("Cannot compute combined k-effective standard deviation with fewer than 4 " - "realizations!"); - - double k_eff[2]; - openmc::openmc_get_keff(k_eff); - return k_eff[1]; - } - default: - mooseError("Internal error: Unhandled StandardDeviationEnum in KEigenvalue!"); - } -} - Real KEigenvalue::kRelativeError() const { - const auto mean = kMean(); - return mean > 0.0 ? KStandardDeviation() / kMean() : 0.0; + const auto mean = kMean(_type); + return mean > 0.0 ? kStandardDeviation(_type) / mean : 0.0; } #endif diff --git a/src/postprocessors/LambdaEffective.C b/src/postprocessors/LambdaEffective.C index bbc17fcdb..8593d33db 100644 --- a/src/postprocessors/LambdaEffective.C +++ b/src/postprocessors/LambdaEffective.C @@ -64,7 +64,7 @@ LambdaEffective::getValue() const const auto den_ss = xt::view(ifp_tally.results_, xt::all(), 2, static_cast(openmc::TallyResult::SUM_SQ)); - const auto mean_k = kMean(); + const auto mean_k = kMean(_type); const auto k_rel = kRelativeError(); const Real lambda_eff = (num_sum[0] / n) / (den_sum[0] / n) / mean_k; diff --git a/src/userobjects/OpenMCNuclideDensities.C b/src/userobjects/OpenMCNuclideDensities.C index e781929bf..949cf55e4 100644 --- a/src/userobjects/OpenMCNuclideDensities.C +++ b/src/userobjects/OpenMCNuclideDensities.C @@ -19,6 +19,7 @@ #ifdef ENABLE_OPENMC_COUPLING #include "OpenMCNuclideDensities.h" +#include "UserErrorChecking.h" #include "openmc/material.h" registerMooseObject("CardinalApp", OpenMCNuclideDensities); @@ -45,9 +46,8 @@ OpenMCNuclideDensities::OpenMCNuclideDensities(const InputParameters & parameter _names(getParam>("names")), _densities(getParam>("densities")) { - _openmc_problem->catchOpenMCError(openmc_get_material_index(_material_id, &_material_index), - "get the material index for material with ID " + - std::to_string(_material_id)); + catchOpenMCError(openmc_get_material_index(_material_id, &_material_index), + "get the material index for material with ID " + std::to_string(_material_id)); } void diff --git a/test/tests/criticality/borated_water/geometry.xml b/test/tests/criticality/borated_water/geometry.xml new file mode 100644 index 000000000..0be48b700 --- /dev/null +++ b/test/tests/criticality/borated_water/geometry.xml @@ -0,0 +1,12 @@ + + + + + + + + + + + + diff --git a/test/tests/criticality/borated_water/gold/openmc_out.csv b/test/tests/criticality/borated_water/gold/openmc_out.csv new file mode 100644 index 000000000..df65100ff --- /dev/null +++ b/test/tests/criticality/borated_water/gold/openmc_out.csv @@ -0,0 +1,3 @@ +time,critical_value,k +0,0,0 +1,522.19791328636,1.0002525559801 diff --git a/test/tests/criticality/borated_water/make_openmc_model.py b/test/tests/criticality/borated_water/make_openmc_model.py new file mode 100644 index 000000000..326e19ebc --- /dev/null +++ b/test/tests/criticality/borated_water/make_openmc_model.py @@ -0,0 +1,68 @@ +#********************************************************************/ +#* SOFTWARE COPYRIGHT NOTIFICATION */ +#* Cardinal */ +#* */ +#* (c) 2021 UChicago Argonne, LLC */ +#* ALL RIGHTS RESERVED */ +#* */ +#* Prepared by UChicago Argonne, LLC */ +#* Under Contract No. DE-AC02-06CH11357 */ +#* With the U. S. Department of Energy */ +#* */ +#* Prepared by Battelle Energy Alliance, LLC */ +#* Under Contract No. DE-AC07-05ID14517 */ +#* With the U. S. Department of Energy */ +#* */ +#* See LICENSE for full restrictions */ +#********************************************************************/ + +import openmc + +def make_uo2_bwater_box(boron_ppm_in=3000): + model = openmc.Model() + + u = openmc.Material() + u.set_density('g/cc', 22.0) + u.add_nuclide('U234', 0.001) + u.add_nuclide('U235', 0.1) + u.add_nuclide('U238', 1.0) + u.add_nuclide('H3', 1e-4) + + bw = openmc.model.borated_water(boron_ppm=boron_ppm_in, density=0.8) + h1_percent = bw.nuclides[0].percent + for n in bw.nuclides: + print(n.name, n.percent / h1_percent) + + model.materials = openmc.Materials([u, bw]) + + # Create two cubes side by side + R = 10.0 + l = openmc.XPlane(x0=0, boundary_type='reflective') + r = openmc.XPlane(x0=R) + b = openmc.YPlane(y0=0, boundary_type='reflective') + t = openmc.YPlane(y0=R, boundary_type='reflective') + f = openmc.ZPlane(z0=0, boundary_type='reflective') + k = openmc.ZPlane(z0=R, boundary_type='reflective') + z = openmc.XPlane(x0=2*R, boundary_type='reflective') + + # create a box + box1 = openmc.Cell(region=+l & -r & +b & -t & +f & -k, fill=u) + box2 = openmc.Cell(region=+r & -z & +b & -t & +f & -k, fill=bw) + model.geometry = openmc.Geometry([box1, box2]) + + # Finally, define some run settings + model.settings.batches = 10 + model.settings.inactive = 3 + model.settings.particles = 1000 + + lower_left = (0, 0, 0) + upper_right = (R, R, R) + uniform_dist = openmc.stats.Box(lower_left, upper_right) + model.settings.source = openmc.source.IndependentSource(space=uniform_dist) + model.settings.temperature = {'default': 600.0, + 'method': 'nearest', + 'range': (294.0, 1600.0)} + +if __name__ == "__main__": + model = make_uo2_bwater_box() + model.export_to_xml() diff --git a/test/tests/criticality/borated_water/materials.xml b/test/tests/criticality/borated_water/materials.xml new file mode 100644 index 000000000..c9584899a --- /dev/null +++ b/test/tests/criticality/borated_water/materials.xml @@ -0,0 +1,20 @@ + + + + + + + + + + + + + + + + + + + + diff --git a/test/tests/criticality/borated_water/openmc.i b/test/tests/criticality/borated_water/openmc.i new file mode 100644 index 000000000..de84d1e6b --- /dev/null +++ b/test/tests/criticality/borated_water/openmc.i @@ -0,0 +1,53 @@ +[Mesh] + [g] + type = GeneratedMeshGenerator + dim = 3 + xmin = 10 + xmax = 20 + ymin = 0 + ymax = 10 + zmin = 0 + zmax = 10 + nx = 10 + ny = 1 + nz = 1 + [] +[] + +[ICs] + [density] + type = ConstantIC + variable = density + value = 1900.0 + [] +[] + +[Problem] + type = OpenMCCellAverageProblem + density_blocks = '0' + cell_level = 0 + verbose = true + + [CriticalitySearch] + type = BoratedWater + absent_nuclides = 'O18' + material_id = 2 + minimum = 0 + maximum = 1000 + tolerance = 5e-1 + [] +[] + +[Executioner] + type = Steady +[] + +[Postprocessors] + [k] + type = KEigenvalue + [] +[] + +[Outputs] + csv = true +[] diff --git a/test/tests/criticality/borated_water/search_keff.py b/test/tests/criticality/borated_water/search_keff.py new file mode 100644 index 000000000..dca90a0b6 --- /dev/null +++ b/test/tests/criticality/borated_water/search_keff.py @@ -0,0 +1,15 @@ +import openmc + +from make_openmc_model import make_uo2_bwater_box + +crit_boron, guesses, keffs = openmc.search_for_keff( + make_uo2_bwater_box, + bracketed_method="brentq", + bracket=[0, 1000], + tol=1e-3, + print_iterations=True) + +for guess, keff in zip(guesses, keffs): + print( guess, keff) + +print("final boron concentration ", crit_boron) \ No newline at end of file diff --git a/test/tests/criticality/borated_water/settings.xml b/test/tests/criticality/borated_water/settings.xml new file mode 100644 index 000000000..158347ae7 --- /dev/null +++ b/test/tests/criticality/borated_water/settings.xml @@ -0,0 +1,15 @@ + + + eigenvalue + 1000 + 10 + 3 + + + 0 0 0 10.0 10.0 10.0 + + + 600.0 + nearest + 294.0 1600.0 + diff --git a/test/tests/criticality/borated_water/tests b/test/tests/criticality/borated_water/tests new file mode 100644 index 000000000..c15c10c43 --- /dev/null +++ b/test/tests/criticality/borated_water/tests @@ -0,0 +1,45 @@ +[Tests] + design = 'BoratedWater.md' + issues = '#1158' + [negative_min] + type = RunException + input = openmc.i + cli_args = 'Problem/CriticalitySearch/minimum=-1' + expect_err = "The 'minimum' boron ppm \(-1.000000\) must be positive!" + requirement = 'The system shall error if an invalid minimum boron ppm is provided.' + capabilities = 'openmc' + [] + [non_dilute] + type = RunException + input = openmc.i + cli_args = 'Problem/CriticalitySearch/maximum=53000' + expect_err = "The borated water composition is computed using a dilute species approximation. Results will not be accurate if the boron species is no longer dilute, which we take as 5\% weight concentration or higher. Please decrease 'maximum' \(53000.000000\) to an upper limit which is in the dilute regime." + requirement = 'The system shall error if an invalid maximum boron ppm is provided.' + capabilities = 'openmc' + [] + [warn_nuclides] + type = RunException + input = openmc.i + cli_args = 'Problem/CriticalitySearch/material_id=1 --error' + expect_err = "The criticality search will clear out all nuclides in material 1 and replace them with the naturally-abundant nuclides in hydrogen, oxygen, and boron. Any other nuclides which existed in the material will be deleted.\n\nThe material you provided contains nuclides which are not the natural isotopes of H, B, and O: U234 U235 U238 H3. These will be deleted from material 1 when the boron concentration is changed." + requirement = "The system shall warn the user if nuclides will be erased by the borated water criticality search. This test covers both non-H,B,O elements as well as non-natural isotopes of these elements." + capabilities = 'openmc' + [] + [absent_nuclides] + type = RunException + input = openmc.i + cli_args = 'Problem/CriticalitySearch/absent_nuclides="W182" --error' + expect_err = "Only absent isotopes of hydrogen, oxygen, or boron will be used to adjust natural abundances. The entry 'W182' will be unused." + requirement = "The system shall warn the user if they are omitting a nuclide from their cross section library which is irrelevant for boric acid control." + capabilities = 'openmc' + [] + [borated_water] + type = CSVDiff + input = openmc.i + csvdiff = openmc_out.csv + requirement = "The system shall perform a criticality search based on boron weight ppm in water. This test changes the water density and then searches for a criticality point; the resulting ppm is compared to a standalone OpenMC model (the comparison is done when running with many more particles for both scenarios). The multiplication factors agree within statistics." + # the relative nuclide concentrations are also printed to the screen and shown to match against + # those added by OpenMC's borated_water function + capabilities = 'openmc' + [] +[] diff --git a/test/tests/criticality/material_density/geometry.xml b/test/tests/criticality/material_density/geometry.xml new file mode 100644 index 000000000..d72e56397 --- /dev/null +++ b/test/tests/criticality/material_density/geometry.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/test/tests/criticality/material_density/gold/openmc_keff0.95_out.csv b/test/tests/criticality/material_density/gold/openmc_keff0.95_out.csv new file mode 100644 index 000000000..b8e23a015 --- /dev/null +++ b/test/tests/criticality/material_density/gold/openmc_keff0.95_out.csv @@ -0,0 +1,3 @@ +time,critical_value,k +0,0,0 +1,16505.455090222,0.94836498988699 diff --git a/test/tests/criticality/material_density/gold/openmc_out.csv b/test/tests/criticality/material_density/gold/openmc_out.csv new file mode 100644 index 000000000..10e30b634 --- /dev/null +++ b/test/tests/criticality/material_density/gold/openmc_out.csv @@ -0,0 +1,3 @@ +time,critical_value,k +0,0,0 +1,16505.455090222,0.99436498988699 diff --git a/test/tests/criticality/material_density/make_openmc_model.py b/test/tests/criticality/material_density/make_openmc_model.py new file mode 100644 index 000000000..05cfa56a2 --- /dev/null +++ b/test/tests/criticality/material_density/make_openmc_model.py @@ -0,0 +1,60 @@ +#********************************************************************/ +#* SOFTWARE COPYRIGHT NOTIFICATION */ +#* Cardinal */ +#* */ +#* (c) 2021 UChicago Argonne, LLC */ +#* ALL RIGHTS RESERVED */ +#* */ +#* Prepared by UChicago Argonne, LLC */ +#* Under Contract No. DE-AC02-06CH11357 */ +#* With the U. S. Department of Energy */ +#* */ +#* Prepared by Battelle Energy Alliance, LLC */ +#* Under Contract No. DE-AC07-05ID14517 */ +#* With the U. S. Department of Energy */ +#* */ +#* See LICENSE for full restrictions */ +#********************************************************************/ + +import openmc + +def make_pu_cube(pu_density= 10.0): + model = openmc.Model() + + pu = openmc.Material() + pu.set_density('g/cc', pu_density) + pu.add_nuclide('Pu239', 1.0) + + model.materials = openmc.Materials([pu]) + + # Create cube + R = 10.0 + l = openmc.XPlane(x0=0, boundary_type='vacuum') + r = openmc.XPlane(x0=R, boundary_type='vacuum') + b = openmc.YPlane(y0=0, boundary_type='vacuum') + t = openmc.YPlane(y0=R, boundary_type='vacuum') + f = openmc.ZPlane(z0=0, boundary_type='vacuum') + k = openmc.ZPlane(z0=R, boundary_type='vacuum') + + # create a box + box = openmc.Cell(region=+l & -r & +b & -t & +f & -k, fill=pu) + model.geometry = openmc.Geometry([box]) + + # Finally, define some run settings + model.settings.batches = 50 + model.settings.inactive = 10 + model.settings.particles = 100 + + lower_left = (0, 0, 0) + upper_right = (R, R, R) + uniform_dist = openmc.stats.Box(lower_left, upper_right) + model.settings.source = openmc.source.IndependentSource(space=uniform_dist) + model.settings.temperature = {'default': 600.0, + 'method': 'nearest', + 'range': (294.0, 1600.0)} + + return model + +if __name__ == "__main__": + model = make_pu_cube() + model.export_to_xml() diff --git a/test/tests/criticality/material_density/materials.xml b/test/tests/criticality/material_density/materials.xml new file mode 100644 index 000000000..12d19c3e2 --- /dev/null +++ b/test/tests/criticality/material_density/materials.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/test/tests/criticality/material_density/openmc.i b/test/tests/criticality/material_density/openmc.i new file mode 100644 index 000000000..9d38f694d --- /dev/null +++ b/test/tests/criticality/material_density/openmc.i @@ -0,0 +1,38 @@ +[Mesh] + [g] + type = GeneratedMeshGenerator + dim = 3 + xmin = 0 + xmax = 10 + ymin = 0 + ymax = 10 + zmin = 0 + zmax = 10 + [] +[] + +[Problem] + type = OpenMCCellAverageProblem + + [CriticalitySearch] + type = OpenMCMaterialDensity + material_id = 1 + minimum = 1000 + maximum = 30000 + tolerance = 1e-2 + [] +[] + +[Executioner] + type = Steady +[] + +[Postprocessors] + [k] + type = KEigenvalue + [] +[] + +[Outputs] + csv = true +[] diff --git a/test/tests/criticality/material_density/search_keff.py b/test/tests/criticality/material_density/search_keff.py new file mode 100644 index 000000000..8223c4c48 --- /dev/null +++ b/test/tests/criticality/material_density/search_keff.py @@ -0,0 +1,22 @@ +import openmc +import argparse + +from make_openmc_model import make_pu_cube + +parser = argparse.ArgumentParser() +parser.add_argument("-t", "--target", type=float, default=1.0) +args = parser.parse_args() + +crit_density, guesses, keffs = openmc.search_for_keff( + make_pu_cube, + target=args.target, + bracketed_method="brentq", + bracket=[1, 30], + tol=1e-5, + print_iterations=True) + +for guess, keff in zip(guesses, keffs): + print( guess, keff) + +print("final density ", crit_density) + diff --git a/test/tests/criticality/material_density/settings.xml b/test/tests/criticality/material_density/settings.xml new file mode 100644 index 000000000..d01a1f523 --- /dev/null +++ b/test/tests/criticality/material_density/settings.xml @@ -0,0 +1,15 @@ + + + eigenvalue + 100 + 50 + 10 + + + 0 0 0 10.0 10.0 10.0 + + + 600.0 + nearest + 294.0 1600.0 + diff --git a/test/tests/criticality/material_density/tests b/test/tests/criticality/material_density/tests new file mode 100644 index 000000000..64e347be6 --- /dev/null +++ b/test/tests/criticality/material_density/tests @@ -0,0 +1,66 @@ +[Tests] + design = 'AddCriticalitySearchAction.md' + issues = '#1158' + [wrong_problem] + type = RunException + input = wrong_problem.i + expect_err = "The \[CriticalitySearch\] block can only be used with wrapped OpenMC cases! You need to change the \[Problem\] block to 'OpenMCCellAverageProblem'." + requirement = 'The system shall error if a criticality search is not paired with the correct problem class.' + capabilities = 'openmc' + [] + [negative_keff_target] + type = RunException + input = openmc.i + cli_args = 'Problem/CriticalitySearch/keff_target=-1.0' + expect_err = "The 'keff_target' value \(-1.000000\) must be greater than 0." + requirement = 'The system shall error if negative values are given for the target keff for the criticality search.' + capabilities = 'openmc' + [] + [min_max] + type = RunException + input = openmc.i + cli_args = 'Problem/CriticalitySearch/minimum=1.0 Problem/CriticalitySearch/maximum=-2.0' + expect_err = "The 'minimum' value \(1.000000\) must be less than the 'maximum' value \(-2.000000\)." + requirement = 'The system shall error if invalid values are provided for the range of values to consider for the criticality search.' + capabilities = 'openmc' + [] + [too_tight_tolerance] + type = RunException + input = openmc.i + cli_args = 'Problem/CriticalitySearch/tolerance=1e-6 --error' + expect_err = "The 'tolerance' for the criticality search \(0.000001\) is smaller than 3-sigma standard deviation in k \(0.003294\); you may have to run a lot of criticality search points to converge to this tolerance." + requirement = 'The system shall warn if the selected tolerance might cause failure to converge due to high statistical noise.' + capabilities = 'openmc' + [] + [fail_converge] + type = RunException + input = openmc.i + cli_args = 'Problem/CriticalitySearch/tolerance=1e-6' + expect_err = "Failed to converge criticality search!" + requirement = 'The system shall error if the criticality search does not converge.' + capabilities = 'openmc' + [] + [search] + type = CSVDiff + input = openmc.i + csvdiff = openmc_out.csv + requirement = 'The system shall conduct a criticality search based on material density.' + capabilities = 'openmc' + [] + [search_nonunity] + type = CSVDiff + input = openmc.i + cli_args = 'Problem/CriticalitySearch/keff_target=0.95' + csvdiff = openmc_keff0.95_out.csv + requirement = 'The system shall conduct an eigenvalue search for a k_eff = 0.95 based on material density.' + capabilities = 'openmc' + [] + [negative_min] + type = RunException + input = openmc.i + cli_args = 'Problem/CriticalitySearch/minimum=-1' + expect_err = "The 'minimum' density \(-1.000000\) must be positive!" + requirement = 'The system shall error if an invalid minimum density is provided.' + capabilities = 'openmc' + [] +[] diff --git a/test/tests/criticality/material_density/wrong_problem.i b/test/tests/criticality/material_density/wrong_problem.i new file mode 100644 index 000000000..0b91f739f --- /dev/null +++ b/test/tests/criticality/material_density/wrong_problem.i @@ -0,0 +1,30 @@ +[Mesh] + [g] + type = GeneratedMeshGenerator + dim = 3 + xmin = 0 + xmax = 10 + ymin = 0 + ymax = 10 + zmin = 0 + zmax = 10 + [] +[] + +[Problem] + type = FEProblem + solve = false + + [CriticalitySearch] + type = OpenMCMaterialDensity + material_id = 1 + minimum = 1000 + maximum = 30000 + tolerance = 1e-2 + [] +[] + +[Executioner] + type = Steady +[] + diff --git a/test/tests/neutronics/dagmc/with_csg/dag_cell_not_root/tests b/test/tests/neutronics/dagmc/with_csg/dag_cell_not_root/tests index bc1923bfa..2dfd27a73 100644 --- a/test/tests/neutronics/dagmc/with_csg/dag_cell_not_root/tests +++ b/test/tests/neutronics/dagmc/with_csg/dag_cell_not_root/tests @@ -1,6 +1,6 @@ [Tests] issues = '#1177' - design = 'MOABSkinner.md OpenMCCellAverageProblem.md' + design = 'MoabSkinner.md OpenMCCellAverageProblem.md' [dag_cell_not_root] type = RunException diff --git a/test/tests/neutronics/dagmc/with_csg/dag_lattice/tests b/test/tests/neutronics/dagmc/with_csg/dag_lattice/tests index 85340c3d0..c64ce6eed 100644 --- a/test/tests/neutronics/dagmc/with_csg/dag_lattice/tests +++ b/test/tests/neutronics/dagmc/with_csg/dag_lattice/tests @@ -1,6 +1,6 @@ [Tests] issues = '#1177' - design = 'MOABSkinner.md OpenMCCellAverageProblem.md' + design = 'MoabSkinner.md OpenMCCellAverageProblem.md' [dag_in_lattice] type = RunException diff --git a/test/tests/neutronics/dagmc/with_csg/dag_lattice_outer/tests b/test/tests/neutronics/dagmc/with_csg/dag_lattice_outer/tests index 9fce9c453..ec9a5198a 100644 --- a/test/tests/neutronics/dagmc/with_csg/dag_lattice_outer/tests +++ b/test/tests/neutronics/dagmc/with_csg/dag_lattice_outer/tests @@ -1,6 +1,6 @@ [Tests] issues = '#1177' - design = 'MOABSkinner.md OpenMCCellAverageProblem.md' + design = 'MoabSkinner.md OpenMCCellAverageProblem.md' [dag_lattice_outer] type = RunException diff --git a/test/tests/neutronics/dagmc/with_csg/mixed_feedback/tests b/test/tests/neutronics/dagmc/with_csg/mixed_feedback/tests index b030055af..a76e8b820 100644 --- a/test/tests/neutronics/dagmc/with_csg/mixed_feedback/tests +++ b/test/tests/neutronics/dagmc/with_csg/mixed_feedback/tests @@ -1,6 +1,6 @@ [Tests] issues = '#1177' - design = 'MOABSkinner.md OpenMCCellAverageProblem.md' + design = 'MoabSkinner.md OpenMCCellAverageProblem.md' [csg_with_dag_feedback] type = RunException diff --git a/test/tests/neutronics/dagmc/with_csg/multi_cell_dag/tests b/test/tests/neutronics/dagmc/with_csg/multi_cell_dag/tests index bb00a63b6..ecff8b546 100644 --- a/test/tests/neutronics/dagmc/with_csg/multi_cell_dag/tests +++ b/test/tests/neutronics/dagmc/with_csg/multi_cell_dag/tests @@ -1,6 +1,6 @@ [Tests] issues = '#1177' - design = 'MOABSkinner.md OpenMCCellAverageProblem.md' + design = 'MoabSkinner.md OpenMCCellAverageProblem.md' [multi_dag_uni_cells] type = RunException diff --git a/test/tests/neutronics/dagmc/with_csg/multiple_dag/tests b/test/tests/neutronics/dagmc/with_csg/multiple_dag/tests index 70f44680c..8444def0c 100644 --- a/test/tests/neutronics/dagmc/with_csg/multiple_dag/tests +++ b/test/tests/neutronics/dagmc/with_csg/multiple_dag/tests @@ -1,6 +1,6 @@ [Tests] issues = '#1177' - design = 'MOABSkinner.md OpenMCCellAverageProblem.md' + design = 'MoabSkinner.md OpenMCCellAverageProblem.md' [multi_dag_uni] type = RunException diff --git a/test/tests/neutronics/dagmc/with_csg/tests b/test/tests/neutronics/dagmc/with_csg/tests index 2c1341bdc..ae4d98a33 100644 --- a/test/tests/neutronics/dagmc/with_csg/tests +++ b/test/tests/neutronics/dagmc/with_csg/tests @@ -1,6 +1,6 @@ [Tests] issues = '#1176 #1177' - design = 'MOABSkinner.md OpenMCCellAverageProblem.md' + design = 'MoabSkinner.md OpenMCCellAverageProblem.md' [allows_csg_with_skinner] type = CSVDiff