diff --git a/docs/source/capi/index.rst b/docs/source/capi/index.rst index d9ac0d1e019..1cdce08b003 100644 --- a/docs/source/capi/index.rst +++ b/docs/source/capi/index.rst @@ -84,6 +84,17 @@ Functions :return: Return status (negative if an error occurred) :rtype: int +.. c:function:: int openmc_cell_get_density(int32_t index, const int32_t* instance, double* rho) + + Get the density of a cell + + :param int32_t index: Index in the cells array + :param int32_t* instance: Which instance of the cell. If a null pointer is passed, the density + multiplier of the first instance is returned. + :param double* rho: density of the cell in g/cc + :return: Return status (negative if an error occurred) + :rtype: int + .. c:function:: int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices) Set the fill for a cell @@ -119,6 +130,20 @@ Functions :return: Return status (negative if an error occurred) :rtype: int +.. c:function:: int openmc_cell_set_density(index index, double rho, const int32_t* instance, bool set_contained) + + Set the density of a cell. + + :param int32_t index: Index in the cells array + :param double rho: Density of the cell in g/cc + :param instance: Which instance of the cell. To set the density multiplier for all + instances, pass a null pointer. + :param set_contained: If the cell is not filled by a material, whether to set the density multiplier + of all filled cells + :type instance: const int32_t* + :return: Return status (negative if an error occurred) + :rtype: int + .. c:function:: int openmc_energy_filter_get_bins(int32_t index, double** energies, int32_t* n) Return the bounding energies for an energy filter diff --git a/docs/source/io_formats/properties.rst b/docs/source/io_formats/properties.rst index 5030e78f35d..e01fe14f3e4 100644 --- a/docs/source/io_formats/properties.rst +++ b/docs/source/io_formats/properties.rst @@ -25,6 +25,7 @@ The current version of the properties file format is 1.0. **/geometry/cells/cell /** :Datasets: - **temperature** (*double[]*) -- Temperature of the cell in [K]. + - **density** (*double[]*) -- Density of the cell in [g/cm3]. **/materials/** diff --git a/docs/source/io_formats/summary.rst b/docs/source/io_formats/summary.rst index 7d3ab94d9f4..78a8876f58d 100644 --- a/docs/source/io_formats/summary.rst +++ b/docs/source/io_formats/summary.rst @@ -38,6 +38,7 @@ The current version of the summary file format is 6.0. is an array if the cell uses distributed materials, otherwise it is a scalar. - **temperature** (*double[]*) -- Temperature of the cell in Kelvin. + - **density** (*double[]*) -- Density of the cell in [g/cm3]. - **translation** (*double[3]*) -- Translation applied to the fill universe. This dataset is present only if fill_type is set to 'universe'. diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 54257d09385..d8041ef4149 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -17,6 +17,8 @@ int openmc_cell_get_fill( int openmc_cell_get_id(int32_t index, int32_t* id); int openmc_cell_get_temperature( int32_t index, const int32_t* instance, double* T); +int openmc_cell_get_density( + int32_t index, const int32_t* instance, double* rho); int openmc_cell_get_translation(int32_t index, double xyz[]); int openmc_cell_get_rotation(int32_t index, double rot[], size_t* n); int openmc_cell_get_name(int32_t index, const char** name); @@ -27,6 +29,8 @@ int openmc_cell_set_fill( int openmc_cell_set_id(int32_t index, int32_t id); int openmc_cell_set_temperature( int32_t index, double T, const int32_t* instance, bool set_contained = false); +int openmc_cell_set_density(int32_t index, double rho, const int32_t* instance, + bool set_contained = false); int openmc_cell_set_translation(int32_t index, const double xyz[]); int openmc_cell_set_rotation(int32_t index, const double rot[], size_t rot_len); int openmc_dagmc_universe_get_cell_ids( diff --git a/include/openmc/cell.h b/include/openmc/cell.h index e9fcfe3911b..cea88b20fe6 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -216,6 +216,18 @@ class Cell { //! \return Temperature in [K] double temperature(int32_t instance = -1) const; + //! Get the density multiplier of a cell instance + //! \param[in] instance Instance index. If -1 is given, the density multiplier + //! for the first instance is returned. + //! \return Density multiplier [-] + double density_mult(int32_t instance = -1) const; + + //! Get the density of a cell instance in g/cm3 + //! \param[in] instance Instance index. If -1 is given, the density + //! for the first instance is returned. + //! \return Density multiplier g/cm3 + double density(int32_t instance = -1) const; + //! Set the temperature of a cell instance //! \param[in] T Temperature in [K] //! \param[in] instance Instance index. If -1 is given, the temperature for @@ -226,6 +238,16 @@ class Cell { void set_temperature( double T, int32_t instance = -1, bool set_contained = false); + //! Set the density of a cell instance + //! \param[in] rho Density [g/cm3] + //! \param[in] instance Instance index. If -1 is given, the density + //! for all instances is set. + //! \param[in] set_contained If this cell is not filled with a material, + //! collect all contained cells with material fills and set their + //! densities. + void set_density( + double rho, int32_t instance = -1, bool set_contained = false); + int32_t n_instances() const; //! Set the rotation matrix of a cell instance @@ -334,6 +356,9 @@ class Cell { //! T. The units are sqrt(eV). vector sqrtkT_; + //! \brief Unitless density multiplier(s) within this cell. + vector rho_mult_ = {1.0}; + //! \brief Neighboring cells in the same universe. NeighborList neighbors_; @@ -351,6 +376,10 @@ class Cell { // Right now, either CSG or DAGMC cells are used. virtual GeometryType geom_type() const = 0; + + //! \brief A flag to indicate if this cell has it's density set in the + //! xml file. + bool xml_set_density_ = false; }; struct CellInstanceItem { diff --git a/include/openmc/constants.h b/include/openmc/constants.h index ae705607958..df13da3707d 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -28,11 +28,11 @@ constexpr int HDF5_VERSION[] {3, 0}; constexpr array VERSION_STATEPOINT {18, 1}; constexpr array VERSION_PARTICLE_RESTART {2, 0}; constexpr array VERSION_TRACK {3, 0}; -constexpr array VERSION_SUMMARY {6, 0}; +constexpr array VERSION_SUMMARY {6, 1}; constexpr array VERSION_VOLUME {1, 0}; constexpr array VERSION_VOXEL {2, 0}; constexpr array VERSION_MGXS_LIBRARY {1, 0}; -constexpr array VERSION_PROPERTIES {1, 0}; +constexpr array VERSION_PROPERTIES {1, 1}; constexpr array VERSION_WEIGHT_WINDOWS {1, 0}; // ============================================================================ diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index f60f2d649ee..4dafdea5c2c 100644 --- a/include/openmc/geometry_aux.h +++ b/include/openmc/geometry_aux.h @@ -37,6 +37,12 @@ void adjust_indices(); void assign_temperatures(); +//============================================================================== +//! Finalize densities (compute density multipliers). +//============================================================================== + +void finalize_cell_densities(); + //============================================================================== //! \brief Obtain a list of temperatures that each nuclide/thermal scattering //! table appears at in the model. Later, this list is used to determine the diff --git a/include/openmc/material.h b/include/openmc/material.h index fe587a86f9e..e36946c71b6 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -99,6 +99,13 @@ class Material { //---------------------------------------------------------------------------- // Accessors + //! Get the atom density in [atom/b-cm] + //! \return Density in [atom/b-cm] + double atom_density(int32_t i, double rho_multiplier = 1.0) const + { + return atom_density_(i) * rho_multiplier; + } + //! Get density in [atom/b-cm] //! \return Density in [atom/b-cm] double density() const { return density_; } diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index ec393845f81..9f50c53dbe4 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -389,6 +389,11 @@ class GeometryState { const double& sqrtkT() const { return sqrtkT_; } double& sqrtkT_last() { return sqrtkT_last_; } + // density multiplier of the current and last cell + double& rho_mult() { return rho_mult_; } + const double& rho_mult() const { return rho_mult_; } + double& rho_mult_last() { return rho_mult_last_; } + private: int64_t id_ {-1}; //!< Unique ID @@ -417,6 +422,9 @@ class GeometryState { double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV double sqrtkT_last_ {0.0}; //!< last temperature + double rho_mult_ {1.0}; //!< density multiplier + double rho_mult_last_ {1.0}; //!< last density multiplier + #ifdef OPENMC_DAGMC_ENABLED moab::DagMC::RayHistory history_; Direction last_dir_; diff --git a/openmc/cell.py b/openmc/cell.py index 88aaebc8f27..1e77fdce629 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -72,6 +72,10 @@ class Cell(IDManagerMixin): temperature : float or iterable of float Temperature of the cell in Kelvin. Multiple temperatures can be given to give each distributed cell instance a unique temperature. + density : float or iterable of float + Density of the cell in g/cm3. Multiple densities can be given to give + each distributed cell instance a unique density. Densities set here will + override the density set on materials used to fill the cell. translation : Iterable of float If the cell is filled with a universe, this array specifies a vector that is used to translate (shift) the universe. @@ -109,6 +113,7 @@ def __init__(self, cell_id=None, name='', fill=None, region=None): self._rotation = None self._rotation_matrix = None self._temperature = None + self._density = None self._translation = None self._paths = None self._num_instances = None @@ -141,6 +146,8 @@ def __repr__(self): if self.fill_type == 'material': string += '\t{0: <15}=\t{1}\n'.format('Temperature', self.temperature) + string += '\t{0: <15}=\t{1}\n'.format('Density', + self.density) string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation) string += '{: <16}=\t{}\n'.format('\tVolume', self.volume) @@ -257,6 +264,30 @@ def temperature(self, temperature): else: self._temperature = temperature + @property + def density(self): + return self._density + + @density.setter + def density(self, density): + # Make sure densities are greater than zero + cv.check_type('cell density', density, (Iterable, Real), none_ok=True) + if isinstance(density, Iterable): + cv.check_type('cell density', density, Iterable, Real) + for rho in density: + cv.check_greater_than('cell density', rho, 0.0, True) + elif isinstance(density, Real): + cv.check_greater_than('cell density', density, 0.0, True) + + # If this cell is filled with a universe or lattice, propagate + # densities to all cells contained. Otherwise, simply assign it. + if self.fill_type in ('universe', 'lattice'): + for c in self.get_all_cells().values(): + if c.fill_type == 'material': + c._density = density + else: + self._density = density + @property def translation(self): return self._translation @@ -525,6 +556,8 @@ def clone(self, clone_materials=True, clone_regions=True, memo=None): clone.volume = self.volume if self.temperature is not None: clone.temperature = self.temperature + if self.density is not None: + clone.density = self.density if self.translation is not None: clone.translation = self.translation if self.rotation is not None: @@ -650,6 +683,13 @@ def create_surface_elements(node, element, memo=None): else: element.set("temperature", str(self.temperature)) + if self.density is not None: + if isinstance(self.density, Iterable): + element.set("density", ' '.join( + str(t) for t in self.density)) + else: + element.set("density", str(self.density)) + if self.translation is not None: element.set("translation", ' '.join(map(str, self.translation))) @@ -711,10 +751,16 @@ def from_xml_element(cls, elem, surfaces, materials, get_universe): c.temperature = temperature else: c.temperature = temperature[0] + density = get_elem_list(elem, 'density', float) + if density is not None: + if len(density) > 1: + c.density = density + else: + c.density = density[0] v = get_text(elem, 'volume') if v is not None: c.volume = float(v) - for key in ('temperature', 'rotation', 'translation'): + for key in ('temperature', 'density', 'rotation', 'translation'): values = get_elem_list(elem, key, float) if values is not None: if key == 'rotation' and len(values) == 9: diff --git a/openmc/lib/cell.py b/openmc/lib/cell.py index 971a24cba91..2d672dbf239 100644 --- a/openmc/lib/cell.py +++ b/openmc/lib/cell.py @@ -34,6 +34,10 @@ c_int32, POINTER(c_int32), POINTER(c_double)] _dll.openmc_cell_get_temperature.restype = c_int _dll.openmc_cell_get_temperature.errcheck = _error_handler +_dll.openmc_cell_get_density.argtypes = [ + c_int32, POINTER(c_int32), POINTER(c_double)] +_dll.openmc_cell_get_density.restype = c_int +_dll.openmc_cell_get_density.errcheck = _error_handler _dll.openmc_cell_get_name.argtypes = [c_int32, POINTER(c_char_p)] _dll.openmc_cell_get_name.restype = c_int _dll.openmc_cell_get_name.errcheck = _error_handler @@ -58,6 +62,10 @@ c_int32, c_double, POINTER(c_int32), c_bool] _dll.openmc_cell_set_temperature.restype = c_int _dll.openmc_cell_set_temperature.errcheck = _error_handler +_dll.openmc_cell_set_density.argtypes = [ + c_int32, c_double, POINTER(c_int32), c_bool] +_dll.openmc_cell_set_density.restype = c_int +_dll.openmc_cell_set_density.errcheck = _error_handler _dll.openmc_cell_set_translation.argtypes = [c_int32, POINTER(c_double)] _dll.openmc_cell_set_translation.restype = c_int _dll.openmc_cell_set_translation.errcheck = _error_handler @@ -236,6 +244,43 @@ def set_temperature(self, T, instance=None, set_contained=False): _dll.openmc_cell_set_temperature(self._index, T, instance, set_contained) + def get_density(self, instance=None): + """Get the density of a cell (g/cc) + + Parameters + ---------- + instance: int or None + Which instance of the cell + + """ + + if instance is not None: + instance = c_int32(instance) + + rho = c_double() + _dll.openmc_cell_get_density(self._index, instance, rho) + return rho.value + + def set_density(self, rho, instance=None, set_contained=False): + """Set the density of a cell + + Parameters + ---------- + rho : float + Density of the cell (g/cc) + instance : int or None + Which instance of the cell + set_contained: bool + If cell is not filled by a material, whether to set the density + of all filled cells + + """ + + if instance is not None: + instance = c_int32(instance) + + _dll.openmc_cell_set_density(self._index, rho, instance, set_contained) + @property def translation(self): translation = np.zeros(3) diff --git a/openmc/model/model.py b/openmc/model/model.py index 03fdcda1fc8..59e6fa511e5 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -633,7 +633,7 @@ def import_properties(self, filename: PathLike): raise ValueError("Number of cells in properties file doesn't " "match current model.") - # Update temperatures for cells filled with materials + # Update temperatures and densities for cells filled with materials for name, group in cells_group.items(): cell_id = int(name.split()[1]) cell = cells[cell_id] @@ -648,6 +648,20 @@ def import_properties(self, filename: PathLike): else: lib_cell.set_temperature(temperature[0]) + if group['density']: + density = group['density'][()] + if density.size > 1: + cell.density = [rho for rho in density] + else: + cell.density = density + if self.is_initialized: + lib_cell = openmc.lib.cells[cell_id] + if density.size > 1: + for i, rho in enumerate(density): + lib_cell.set_density(rho, i) + else: + lib_cell.set_density(density[0]) + # Make sure number of materials matches mats_group = fh['materials'] n_cells = mats_group.attrs['n_materials'] diff --git a/src/cell.cpp b/src/cell.cpp index d838bfbb4d1..f482ba507b3 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -101,6 +101,27 @@ double Cell::temperature(int32_t instance) const } } +double Cell::density_mult(int32_t instance) const +{ + if (instance >= 0) { + double rho = + rho_mult_.size() == 1 ? rho_mult_.at(0) : rho_mult_.at(instance); + return rho; + } else { + return rho_mult_[0]; + } +} + +double Cell::density(int32_t instance) const +{ + const int32_t mat_index = material(instance); + if (mat_index == MATERIAL_VOID) + return 0.0; + + double rho_mult = density_mult(instance); + return rho_mult * model::materials[mat_index]->density_gpcc(); +} + void Cell::set_temperature(double T, int32_t instance, bool set_contained) { if (settings::temperature_method == TemperatureMethod::INTERPOLATION) { @@ -151,6 +172,47 @@ void Cell::set_temperature(double T, int32_t instance, bool set_contained) } } +void Cell::set_density(double rho, int32_t instance, bool set_contained) +{ + if (type_ != Fill::MATERIAL && !set_contained) { + fatal_error( + fmt::format("Attempted to set the density multiplier of cell {} " + "which is not filled by a material.", + id_)); + } + + if (type_ == Fill::MATERIAL) { + const int32_t mat_index = material(instance); + if (mat_index == MATERIAL_VOID) + return; + + if (instance >= 0) { + // If density multiplier vector is not big enough, resize it first + if (rho_mult_.size() != n_instances()) + rho_mult_.resize(n_instances(), rho_mult_[0]); + + // Set density multiplier for the corresponding instance + rho_mult_.at(instance) = + rho / model::materials[mat_index]->density_gpcc(); + } else { + // Set density multiplier for all instances + for (auto& Rho_ : rho_mult_) { + Rho_ = rho / model::materials[mat_index]->density_gpcc(); + } + } + } else { + auto contained_cells = this->get_contained_cells(instance); + for (const auto& entry : contained_cells) { + auto& cell = model::cells[entry.first]; + assert(cell->type_ == Fill::MATERIAL); + auto& instances = entry.second; + for (auto instance : instances) { + cell->set_density(rho, instance); + } + } + } +} + void Cell::export_properties_hdf5(hid_t group) const { // Create a group for this cell. @@ -162,6 +224,15 @@ void Cell::export_properties_hdf5(hid_t group) const temps.push_back(sqrtkT_val * sqrtkT_val / K_BOLTZMANN); write_dataset(cell_group, "temperature", temps); + // Write density for one or more cell instances + if (type_ == Fill::MATERIAL && material_.size() > 0) { + vector rho; + for (int32_t i = 0; i < rho_mult_.size(); ++i) + rho.push_back(density(i)); + + write_dataset(cell_group, "density", rho); + } + close_group(cell_group); } @@ -176,7 +247,7 @@ void Cell::import_properties_hdf5(hid_t group) // Ensure number of temperatures makes sense auto n_temps = temps.size(); if (n_temps > 1 && n_temps != n_instances()) { - throw std::runtime_error(fmt::format( + fatal_error(fmt::format( "Number of temperatures for cell {} doesn't match number of instances", id_)); } @@ -188,6 +259,24 @@ void Cell::import_properties_hdf5(hid_t group) this->set_temperature(temps[i], i); } + // Read densities + if (object_exists(cell_group, "density")) { + vector rho; + read_dataset(cell_group, "density", rho); + + // Ensure number of densities makes sense + auto n_rho = rho.size(); + if (n_rho > 1 && n_rho != n_instances()) { + fatal_error(fmt::format("Number of densities for cell {} " + "doesn't match number of instances", + id_)); + } + + // Set densities. + for (int32_t i = 0; i < n_rho; ++i) + set_density(rho[i], i); + } + close_group(cell_group); } @@ -227,6 +316,8 @@ void Cell::to_hdf5(hid_t cell_group) const temps.push_back(sqrtkT_val * sqrtkT_val / K_BOLTZMANN); write_dataset(group, "temperature", temps); + write_dataset(group, "density_mult", rho_mult_); + } else if (type_ == Fill::UNIVERSE) { write_dataset(group, "fill_type", "universe"); write_dataset(group, "fill", model::universes[fill_]->id_); @@ -344,6 +435,45 @@ CSGCell::CSGCell(pugi::xml_node cell_node) } } + // Read the density element which can be distributed similar to temperature. + // These get assigned to the density multiplier, requiring a division by + // the material density. + // Note: calculating the actual density multiplier is deferred until materials + // are finalized. rho_mult_ contains the true density in the meantime. + if (check_for_node(cell_node, "density")) { + rho_mult_ = get_node_array(cell_node, "density"); + rho_mult_.shrink_to_fit(); + xml_set_density_ = true; + + // Make sure this is a material-filled cell. + if (material_.size() == 0) { + fatal_error(fmt::format( + "Cell {} was specified with a density but no material. Density" + "specification is only valid for cells filled with a material.", + id_)); + } + + // Make sure this is a non-void material. + for (auto mat_id : material_) { + if (mat_id == MATERIAL_VOID) { + fatal_error(fmt::format( + "Cell {} was specified with a density, but contains a void " + "material. Density specification is only valid for cells " + "filled with a non-void material.", + id_)); + } + } + + // Make sure all densities are non-negative and greater than zero. + for (auto rho : rho_mult_) { + if (rho <= 0) { + fatal_error(fmt::format( + "Cell {} was specified with a density less than or equal to zero", + id_)); + } + } + } + // Read the region specification. std::string region_spec; if (check_for_node(cell_node, "region")) { @@ -1129,6 +1259,24 @@ extern "C" int openmc_cell_set_temperature( return 0; } +extern "C" int openmc_cell_set_density( + int32_t index, double rho, const int32_t* instance, bool set_contained) +{ + if (index < 0 || index >= model::cells.size()) { + strcpy(openmc_err_msg, "Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + int32_t instance_index = instance ? *instance : -1; + try { + model::cells[index]->set_density(rho, instance_index, set_contained); + } catch (const std::exception& e) { + set_errmsg(e.what()); + return OPENMC_E_UNASSIGNED; + } + return 0; +} + extern "C" int openmc_cell_get_temperature( int32_t index, const int32_t* instance, double* T) { @@ -1147,6 +1295,36 @@ extern "C" int openmc_cell_get_temperature( return 0; } +extern "C" int openmc_cell_get_density( + int32_t index, const int32_t* instance, double* rho) +{ + if (index < 0 || index >= model::cells.size()) { + strcpy(openmc_err_msg, "Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + int32_t instance_index = instance ? *instance : -1; + try { + if (model::cells[index]->type_ != Fill::MATERIAL) { + fatal_error( + fmt::format("Cell {}, instance {} is not filled with a material.", + model::cells[index]->id_, instance_index)); + } + + int32_t mat_index = model::cells[index]->material(instance_index); + if (mat_index == MATERIAL_VOID) { + *rho = 0.0; + } else { + *rho = model::cells[index]->density_mult(instance_index) * + model::materials[mat_index]->density_gpcc(); + } + } catch (const std::exception& e) { + set_errmsg(e.what()); + return OPENMC_E_UNASSIGNED; + } + return 0; +} + //! Get the bounding box of a cell extern "C" int openmc_cell_bounding_box( const int32_t index, double* llc, double* urc) diff --git a/src/geometry.cpp b/src/geometry.cpp index df485089163..f4a5b5cecaf 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -172,11 +172,13 @@ bool find_cell_inner( p.cell_instance() = cell_instance_at_level(p, p.n_coord() - 1); } - // Set the material and temperature. + // Set the material, temperature and density multiplier. p.material_last() = p.material(); p.material() = c.material(p.cell_instance()); p.sqrtkT_last() = p.sqrtkT(); p.sqrtkT() = c.sqrtkT(p.cell_instance()); + p.rho_mult_last() = p.rho_mult(); + p.rho_mult() = c.density_mult(p.cell_instance()); return true; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 51732cf8b90..5320146486d 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -195,6 +195,21 @@ void assign_temperatures() //============================================================================== +void finalize_cell_densities() +{ + for (auto& c : model::cells) { + // Convert to density multipliers. + if (c->xml_set_density_) { + for (int32_t instance = 0; instance < c->rho_mult_.size(); ++instance) { + c->rho_mult_[instance] /= + model::materials[c->material(instance)]->density_gpcc(); + } + } + } +} + +//============================================================================== + void get_temperatures( vector>& nuc_temps, vector>& thermal_temps) { @@ -362,6 +377,17 @@ void prepare_distribcell(const std::vector* user_distribcells) c.id_, c.sqrtkT_.size(), c.n_instances())); } } + + if (c.rho_mult_.size() > 1) { + if (c.rho_mult_.size() != c.n_instances()) { + fatal_error(fmt::format("Cell {} was specified with {} density " + "multipliers but has {} distributed " + "instances. The number of density multipliers " + "must equal one or the number " + "of instances.", + c.id_, c.rho_mult_.size(), c.n_instances())); + } + } } // Search through universes for material cells and assign each one a diff --git a/src/initialize.cpp b/src/initialize.cpp index 2da78137d10..c77fcd02ddb 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -401,6 +401,10 @@ bool read_model_xml() // Finalize cross sections having assigned temperatures finalize_cross_sections(); + // Compute cell density multipliers now that material densities + // have been finalized (from geometry_aux.h) + finalize_cell_densities(); + if (check_for_node(root, "tallies")) read_tallies_xml(root.child("tallies")); diff --git a/src/material.cpp b/src/material.cpp index 32384ddf12b..a6a4cb27ec6 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -890,7 +890,7 @@ void Material::calculate_neutron_xs(Particle& p) const // ADD TO MACROSCOPIC CROSS SECTION // Copy atom density of nuclide in material - double atom_density = atom_density_(i); + double atom_density = this->atom_density(i, p.rho_mult()); // Add contributions to cross sections p.macro_xs().total += atom_density * micro.total; @@ -925,7 +925,7 @@ void Material::calculate_photon_xs(Particle& p) const // ADD TO MACROSCOPIC CROSS SECTION // Copy atom density of nuclide in material - double atom_density = atom_density_(i); + double atom_density = this->atom_density(i, p.rho_mult()); // Add contributions to material macroscopic cross sections p.macro_xs().total += atom_density * micro.total; diff --git a/src/mgxs.cpp b/src/mgxs.cpp index a88d2c196b6..71386a3fd31 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -617,10 +617,11 @@ void Mgxs::calculate_xs(Particle& p) } int temperature = p.mg_xs_cache().t; int angle = p.mg_xs_cache().a; - p.macro_xs().total = xs[temperature].total(angle, p.g()); - p.macro_xs().absorption = xs[temperature].absorption(angle, p.g()); + p.macro_xs().total = xs[temperature].total(angle, p.g()) * p.rho_mult(); + p.macro_xs().absorption = + xs[temperature].absorption(angle, p.g()) * p.rho_mult(); p.macro_xs().nu_fission = - fissionable ? xs[temperature].nu_fission(angle, p.g()) : 0.; + fissionable ? xs[temperature].nu_fission(angle, p.g()) * p.rho_mult() : 0.; } //============================================================================== diff --git a/src/particle.cpp b/src/particle.cpp index fb41d82e950..2a22b0ae89f 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -200,7 +200,8 @@ void Particle::event_calculate_xs() // Calculate microscopic and macroscopic cross sections if (material() != MATERIAL_VOID) { if (settings::run_CE) { - if (material() != material_last() || sqrtkT() != sqrtkT_last()) { + if (material() != material_last() || sqrtkT() != sqrtkT_last() || + rho_mult() != rho_mult_last()) { // If the material is the same as the last material and the // temperature hasn't changed, we don't need to lookup cross // sections again. @@ -558,9 +559,10 @@ void Particle::cross_surface(const Surface& surf) int32_t i_cell = next_cell(surface_index(), cell_last(n_coord() - 1), lowest_coord().universe()) - 1; - // save material and temp + // save material, temp and density multiplier material_last() = material(); sqrtkT_last() = sqrtkT(); + rho_mult_last() = rho_mult(); // set new cell value lowest_coord().cell() = i_cell; auto& cell = model::cells[i_cell]; @@ -571,6 +573,7 @@ void Particle::cross_surface(const Surface& surf) material() = cell->material(cell_instance()); sqrtkT() = cell->sqrtkT(cell_instance()); + rho_mult() = cell->density_mult(cell_instance()); return; } #endif diff --git a/src/physics.cpp b/src/physics.cpp index 3c06e543de1..7a4949e2a26 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -496,7 +496,7 @@ int sample_nuclide(Particle& p) for (int i = 0; i < n; ++i) { // Get atom density int i_nuclide = mat->nuclide_[i]; - double atom_density = mat->atom_density_[i]; + double atom_density = mat->atom_density(i, p.rho_mult()); // Increment probability to compare to cutoff prob += atom_density * p.neutron_xs(i_nuclide).total; @@ -521,7 +521,7 @@ int sample_element(Particle& p) for (int i = 0; i < mat->element_.size(); ++i) { // Find atom density int i_element = mat->element_[i]; - double atom_density = mat->atom_density_[i]; + double atom_density = mat->atom_density(i, p.rho_mult()); // Determine microscopic cross section double sigma = atom_density * p.photon_xs(i_element).total; diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index e73fb90f311..f57efd66405 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -233,7 +233,7 @@ double score_fission_q(const Particle& p, int score_bin, const Tally& tally, double score {0.0}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); + auto atom_density = material.atom_density(i, p.rho_mult()); const Nuclide& nuc {*data::nuclides[j_nuclide]}; score += get_nuc_fission_q(nuc, p, score_bin) * atom_density * p.neutron_xs(j_nuclide).fission; @@ -696,7 +696,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); + auto atom_density = material.atom_density(i, p.rho_mult()); score += p.neutron_xs(j_nuclide).fission * data::nuclides[j_nuclide]->nu( E, ReactionProduct::EmissionMode::prompt) * @@ -743,7 +743,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); + auto atom_density = material.atom_density(i, p.rho_mult()); // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin]; @@ -763,7 +763,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); + auto atom_density = material.atom_density(i, p.rho_mult()); score += p.neutron_xs(j_nuclide).fission * data::nuclides[j_nuclide]->nu( E, ReactionProduct::EmissionMode::delayed) * @@ -824,7 +824,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); + auto atom_density = material.atom_density(i, p.rho_mult()); const auto& nuc {*data::nuclides[j_nuclide]}; if (nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; @@ -849,7 +849,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); + auto atom_density = material.atom_density(i, p.rho_mult()); const auto& nuc {*data::nuclides[j_nuclide]}; if (nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; @@ -893,7 +893,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); + auto atom_density = material.atom_density(i, p.rho_mult()); const auto& nuc {*data::nuclides[j_nuclide]}; if (nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; @@ -924,7 +924,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); + auto atom_density = material.atom_density(i, p.rho_mult()); if (p.neutron_xs(j_nuclide).elastic == CACHE_INVALID) data::nuclides[j_nuclide]->calculate_elastic_xs(p); score += p.neutron_xs(j_nuclide).elastic * atom_density * flux; @@ -1025,7 +1025,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); + auto atom_density = material.atom_density(i, p.rho_mult()); score += p.neutron_xs(j_nuclide).reaction[m] * atom_density * flux; } } @@ -1079,7 +1079,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); + auto atom_density = material.atom_density(i, p.rho_mult()); score += get_nuclide_xs(p, j_nuclide, score_bin) * atom_density * flux; } @@ -2383,7 +2383,8 @@ void score_analog_tally_mg(Particle& p) model::materials[p.material()]->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; - atom_density = model::materials[p.material()]->atom_density_(j); + atom_density = + model::materials[p.material()]->atom_density(j, p.rho_mult()); } score_general_mg(p, i_tally, i * tally.scores_.size(), filter_index, @@ -2451,8 +2452,9 @@ void score_tracklength_tally(Particle& p, double distance) atom_density = 1.0; } } else { - atom_density = - tally.multiply_density() ? mat->atom_density_(j) : 1.0; + atom_density = tally.multiply_density() + ? mat->atom_density(j, p.rho_mult()) + : 1.0; } } } @@ -2530,8 +2532,9 @@ void score_collision_tally(Particle& p) atom_density = 1.0; } } else { - atom_density = - tally.multiply_density() ? mat->atom_density_(j) : 1.0; + atom_density = tally.multiply_density() + ? mat->atom_density(j, p.rho_mult()) + : 1.0; } } diff --git a/tests/regression_tests/cpp_driver/driver.cpp b/tests/regression_tests/cpp_driver/driver.cpp index 48ed7f3171b..a99c97b64e4 100644 --- a/tests/regression_tests/cpp_driver/driver.cpp +++ b/tests/regression_tests/cpp_driver/driver.cpp @@ -15,14 +15,16 @@ using namespace openmc; -int main(int argc, char** argv) { +int main(int argc, char** argv) +{ #ifdef OPENMC_MPI MPI_Comm world {MPI_COMM_WORLD}; int err = openmc_init(argc, argv, &world); #else int err = openmc_init(argc, argv, nullptr); #endif - if (err) fatal_error(openmc_err_msg); + if (err) + fatal_error(openmc_err_msg); // create a new cell filter auto cell_filter = Filter::create(); @@ -30,7 +32,7 @@ int main(int argc, char** argv) { // add all cells to the cell filter std::vector cell_indices; for (auto& entry : openmc::model::cell_map) { - cell_indices.push_back(entry.second); + cell_indices.push_back(entry.second); } // enable distribcells offsets for all cells prepare_distribcell(&cell_indices); @@ -39,7 +41,6 @@ int main(int argc, char** argv) { std::sort(cell_indices.begin(), cell_indices.end()); cell_filter->set_cells(cell_indices); - // create a new tally auto tally = Tally::create(); std::vector filters = {cell_filter}; @@ -60,14 +61,19 @@ int main(int argc, char** argv) { } } - // set a higher temperature for only one of the lattice cells (ID is 4 in the model) + // set a higher temperature for only one of the lattice cells (ID is 4 in the + // model) model::cells[model::cell_map[4]]->set_temperature(400.0, 3, true); + // set the density of another lattice cell to 2 + model::cells[model::cell_map[4]]->set_density(2.0, 2, true); + // the summary file will be used to check that // temperatures were set correctly so clear // error output can be provided #ifdef OPENMC_MPI - if (openmc::mpi::master) openmc::write_summary(); + if (openmc::mpi::master) + openmc::write_summary(); #else openmc::write_summary(); #endif diff --git a/tests/regression_tests/cpp_driver/results_true.dat b/tests/regression_tests/cpp_driver/results_true.dat index c2b0b8d6f36..26f07842cd4 100644 --- a/tests/regression_tests/cpp_driver/results_true.dat +++ b/tests/regression_tests/cpp_driver/results_true.dat @@ -1,13 +1,13 @@ k-combined: -1.933305E+00 1.300360E-02 +1.953962E+00 1.828426E-02 tally 1: -9.552846E+01 -1.019358E+03 -2.887973E+01 -9.308509E+01 -9.732441E+01 -1.059022E+03 -2.217326E+02 -5.486892E+03 -2.217326E+02 -5.486892E+03 +9.607953E+01 +1.031898E+03 +2.853683E+01 +9.085469E+01 +9.745011E+01 +1.058928E+03 +2.220665E+02 +5.496813E+03 +2.220665E+02 +5.496813E+03 diff --git a/tests/regression_tests/lattice_distribrho/__init__.py b/tests/regression_tests/lattice_distribrho/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/lattice_distribrho/inputs_true.dat b/tests/regression_tests/lattice_distribrho/inputs_true.dat new file mode 100644 index 00000000000..5031bea6e28 --- /dev/null +++ b/tests/regression_tests/lattice_distribrho/inputs_true.dat @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + 1.0 1.0 + 2 2 + -1.0 -1.0 + +1 1 +1 1 + + + + + + + + + eigenvalue + 1000 + 10 + 5 + + diff --git a/tests/regression_tests/lattice_distribrho/results_true.dat b/tests/regression_tests/lattice_distribrho/results_true.dat new file mode 100644 index 00000000000..d7094f4e373 --- /dev/null +++ b/tests/regression_tests/lattice_distribrho/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.904471E+00 5.255549E-03 diff --git a/tests/regression_tests/lattice_distribrho/test.py b/tests/regression_tests/lattice_distribrho/test.py new file mode 100644 index 00000000000..c43ee213ff7 --- /dev/null +++ b/tests/regression_tests/lattice_distribrho/test.py @@ -0,0 +1,45 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + +@pytest.fixture +def model(): + model = openmc.model.Model() + + uo2 = openmc.Material(name='UO2') + uo2.set_density('g/cm3', 10.0) + uo2.add_nuclide('U235', 1.0) + uo2.add_nuclide('O16', 2.0) + water = openmc.Material(name='light water') + water.add_nuclide('H1', 2.0) + water.add_nuclide('O16', 1.0) + water.set_density('g/cm3', 1.0) + water.add_s_alpha_beta('c_H_in_H2O') + model.materials.extend([uo2, water]) + + cyl = openmc.ZCylinder(r=0.4) + pin = openmc.model.pin([cyl], [uo2, water]) + d = 1.0 + + lattice = openmc.RectLattice() + lattice.lower_left = (-d, -d) + lattice.pitch = (d, d) + lattice.universes = [[pin, pin], + [pin, pin]] + box = openmc.model.RectangularPrism(2.0 * d, 2.0 * d, origin=(0.0,0.0), boundary_type='reflective') + + pin.cells[1].density = [10.0, 20.0, 10.0, 20.0] + + model.geometry = openmc.Geometry([openmc.Cell(fill=lattice, region = -box)]) + model.geometry.merge_surfaces = True + + model.settings.batches = 10 + model.settings.inactive = 5 + model.settings.particles = 1000 + + return model + +def test_lattice_checkerboard(model): + harness = PyAPITestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/regression_tests/multipole/results_true.dat b/tests/regression_tests/multipole/results_true.dat index 8e82655ce50..8eb75b894bc 100644 --- a/tests/regression_tests/multipole/results_true.dat +++ b/tests/regression_tests/multipole/results_true.dat @@ -38,5 +38,6 @@ Cell Region = -1 Rotation = None Temperature = [500. 700. 0. 800.] + Density = None Translation = None Volume = None diff --git a/tests/unit_tests/test_cell.py b/tests/unit_tests/test_cell.py index 95c8249bb33..60b20581868 100644 --- a/tests/unit_tests/test_cell.py +++ b/tests/unit_tests/test_cell.py @@ -126,6 +126,29 @@ def test_temperature(cell_with_lattice): c.temperature = (300., 600., 900.) +def test_densities(cell_with_lattice): + # Make sure density propagates through universes + m = openmc.Material() + s = openmc.XPlane() + c1 = openmc.Cell(fill=m, region=+s) + c2 = openmc.Cell(fill=m, region=-s) + u1 = openmc.Universe(cells=[c1, c2]) + c = openmc.Cell(fill=u1) + + c.density = 1. + assert c1.density == 1. + assert c2.density == 1. + with pytest.raises(ValueError): + c.density = -1. + c.density = None + assert c1.density == None + assert c2.density == None + + # distributed density + cells, _, _, _ = cell_with_lattice + c = cells[0] + c.density = (1., 2., 3.) + def test_rotation(): u = openmc.Universe() c = openmc.Cell(fill=u) diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index 8ab35335fc0..8080ce89483 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -159,6 +159,29 @@ def test_properties_temperature(lib_init): assert cell.get_temperature() == pytest.approx(200.0) +def test_cell_density(lib_init): + cell = openmc.lib.cells[1] + cell.set_density(1.5, 0) + assert cell.get_density(0) == pytest.approx(1.5) + cell.set_density(2.0) + assert cell.get_density() == pytest.approx(2.0) + + +def test_properties_cell_density(lib_init): + # Cell density should be 2.0 from above test + cell = openmc.lib.cells[1] + assert cell.get_density() == pytest.approx(2.0) + + # Export properties and change density + openmc.lib.export_properties('properties.h5') + cell.set_density(3.0) + assert cell.get_density() == pytest.approx(3.0) + + # Import properties and check that density is restored + openmc.lib.import_properties('properties.h5') + assert cell.get_density() == pytest.approx(2.0) + + def test_new_cell(lib_init): with pytest.raises(exc.AllocationError): openmc.lib.Cell(1) diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py index 6e4dec00fff..12f0df6c9e5 100644 --- a/tests/unit_tests/test_model.py +++ b/tests/unit_tests/test_model.py @@ -251,10 +251,11 @@ def test_import_properties(run_in_tmpdir, mpi_intracomm): model = openmc.examples.pwr_pin_cell() model.init_lib(output=False, intracomm=mpi_intracomm) - # Change fuel temperature and density and export properties + # Change cell fuel temperature, density, material density and export properties cell = openmc.lib.cells[1] cell.set_temperature(600.0) cell.fill.set_density(5.0, 'g/cm3') + cell.set_density(10.0) openmc.lib.export_properties(output=False) # Import properties to existing model @@ -264,9 +265,11 @@ def test_import_properties(run_in_tmpdir, mpi_intracomm): # First python cell = model.geometry.get_all_cells()[1] assert cell.temperature == [600.0] + assert cell.density == [pytest.approx(10.0, 1e-5)] assert cell.fill.get_mass_density() == pytest.approx(5.0) # Now C assert openmc.lib.cells[1].get_temperature() == 600. + assert openmc.lib.cells[1].get_density() == pytest.approx(10.0, 1e-5) assert openmc.lib.materials[1].get_density('g/cm3') == pytest.approx(5.0) # Clear the C API @@ -283,6 +286,7 @@ def test_import_properties(run_in_tmpdir, mpi_intracomm): ) cell = model_with_properties.geometry.get_all_cells()[1] assert cell.temperature == [600.0] + assert cell.density == [pytest.approx(10.0, 1e-5)] assert cell.fill.get_mass_density() == pytest.approx(5.0)