Fixed a bug in distribcell offsets logic (#3424)

This commit is contained in:
GuySten 2025-07-30 20:10:34 +03:00 committed by GitHub
parent 6f5e1347d4
commit a64cc96ed5
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
16 changed files with 325 additions and 87 deletions

View file

@ -226,6 +226,8 @@ public:
void set_temperature(
double T, int32_t instance = -1, bool set_contained = false);
int32_t n_instances() const;
//! Set the rotation matrix of a cell instance
//! \param[in] rot The rotation matrix of length 3 or 9
void set_rotation(const vector<double>& rot);
@ -312,12 +314,11 @@ public:
//----------------------------------------------------------------------------
// Data members
int32_t id_; //!< Unique ID
std::string name_; //!< User-defined name
Fill type_; //!< Material, universe, or lattice
int32_t universe_; //!< Universe # this cell is in
int32_t fill_; //!< Universe # filling this cell
int32_t n_instances_ {0}; //!< Number of instances of this cell
int32_t id_; //!< Unique ID
std::string name_; //!< User-defined name
Fill type_; //!< Material, universe, or lattice
int32_t universe_; //!< Universe # this cell is in
int32_t fill_; //!< Universe # filling this cell
//! \brief Index corresponding to this cell in distribcell arrays
int distribcell_index_ {C_NONE};

View file

@ -15,8 +15,6 @@
namespace openmc {
namespace model {
extern std::unordered_map<int32_t, std::unordered_map<int32_t, int32_t>>
universe_cell_counts;
extern std::unordered_map<int32_t, int32_t> universe_level_counts;
} // namespace model
@ -80,15 +78,13 @@ void prepare_distribcell(
const std::vector<int32_t>* user_distribcells = nullptr);
//==============================================================================
//! Recursively search through the geometry and count cell instances.
//! Recursively search through the geometry and count universe instances.
//!
//! This function will update the Cell::n_instances value for each cell in the
//! geometry.
//! \param univ_indx The index of the universe to begin searching from (probably
//! the root universe).
//! This function will update Universe.n_instances_ for each
//! universe in the geometry.
//==============================================================================
void count_cell_instances(int32_t univ_indx);
void count_universe_instances();
//==============================================================================
//! Recursively search through universes and count universe instances.

View file

@ -76,7 +76,7 @@ public:
}
//! Populate the distribcell offset tables.
int32_t fill_offset_table(int32_t offset, int32_t target_univ_id, int map,
int32_t fill_offset_table(int32_t target_univ_id, int map,
std::unordered_map<int32_t, int32_t>& univ_count_memo);
//! \brief Check lattice indices.

View file

@ -29,6 +29,7 @@ class Universe {
public:
int32_t id_; //!< Unique ID
vector<int32_t> cells_; //!< Cells within this universe
int32_t n_instances_; //!< Number of instances of this universe
//! \brief Write universe information to an HDF5 group.
//! \param group_id An HDF5 group id.

View file

@ -40,6 +40,11 @@ vector<unique_ptr<Cell>> cells;
// Cell implementation
//==============================================================================
int32_t Cell::n_instances() const
{
return model::universes[universe_]->n_instances_;
}
void Cell::set_rotation(const vector<double>& rot)
{
if (fill_ == C_NONE) {
@ -115,8 +120,8 @@ void Cell::set_temperature(double T, int32_t instance, bool set_contained)
if (type_ == Fill::MATERIAL) {
if (instance >= 0) {
// If temperature vector is not big enough, resize it first
if (sqrtkT_.size() != n_instances_)
sqrtkT_.resize(n_instances_, sqrtkT_[0]);
if (sqrtkT_.size() != n_instances())
sqrtkT_.resize(n_instances(), sqrtkT_[0]);
// Set temperature for the corresponding instance
sqrtkT_.at(instance) = std::sqrt(K_BOLTZMANN * T);
@ -170,7 +175,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_) {
if (n_temps > 1 && n_temps != n_instances()) {
throw std::runtime_error(fmt::format(
"Number of temperatures for cell {} doesn't match number of instances",
id_));
@ -1618,7 +1623,7 @@ extern "C" int openmc_cell_get_num_instances(
set_errmsg("Index in cells array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
*num_instances = model::cells[index]->n_instances_;
*num_instances = model::cells[index]->n_instances();
return 0;
}

View file

@ -185,7 +185,6 @@ int openmc_finalize()
int openmc_reset()
{
model::universe_cell_counts.clear();
model::universe_level_counts.clear();
for (auto& t : model::tallies) {

View file

@ -25,21 +25,9 @@
namespace openmc {
namespace model {
std::unordered_map<int32_t, std::unordered_map<int32_t, int32_t>>
universe_cell_counts;
std::unordered_map<int32_t, int32_t> universe_level_counts;
} // namespace model
// adds the cell counts of universe b to universe a
void update_universe_cell_count(int32_t a, int32_t b)
{
auto& universe_a_counts = model::universe_cell_counts[a];
const auto& universe_b_counts = model::universe_cell_counts[b];
for (const auto& it : universe_b_counts) {
universe_a_counts[it.first] += it.second;
}
}
void read_geometry_xml()
{
// Display output message
@ -263,7 +251,7 @@ void finalize_geometry()
{
// Perform some final operations to set up the geometry
adjust_indices();
count_cell_instances(model::root_universe);
count_universe_instances();
partition_universes();
// Assign temperatures to cells that don't have temperatures already assigned
@ -356,35 +344,38 @@ void prepare_distribcell(const std::vector<int32_t>* user_distribcells)
Cell& c {*model::cells[i]};
if (c.material_.size() > 1) {
if (c.material_.size() != c.n_instances_) {
if (c.material_.size() != c.n_instances()) {
fatal_error(fmt::format(
"Cell {} was specified with {} materials but has {} distributed "
"instances. The number of materials must equal one or the number "
"of instances.",
c.id_, c.material_.size(), c.n_instances_));
c.id_, c.material_.size(), c.n_instances()));
}
}
if (c.sqrtkT_.size() > 1) {
if (c.sqrtkT_.size() != c.n_instances_) {
if (c.sqrtkT_.size() != c.n_instances()) {
fatal_error(fmt::format(
"Cell {} was specified with {} temperatures but has {} distributed "
"instances. The number of temperatures must equal one or the number "
"of instances.",
c.id_, c.sqrtkT_.size(), c.n_instances_));
c.id_, c.sqrtkT_.size(), c.n_instances()));
}
}
}
// Search through universes for material cells and assign each one a
// unique distribcell array index.
int distribcell_index = 0;
// distribcell array index according to the containing universe.
vector<int32_t> target_univ_ids;
for (const auto& u : model::universes) {
for (auto idx : u->cells_) {
if (distribcells.find(idx) != distribcells.end()) {
model::cells[idx]->distribcell_index_ = distribcell_index++;
target_univ_ids.push_back(u->id_);
if (!contains(target_univ_ids, u->id_)) {
target_univ_ids.push_back(u->id_);
}
model::cells[idx]->distribcell_index_ =
std::find(target_univ_ids.begin(), target_univ_ids.end(), u->id_) -
target_univ_ids.begin();
}
}
}
@ -419,8 +410,7 @@ void prepare_distribcell(const std::vector<int32_t>* user_distribcells)
} else if (c.type_ == Fill::LATTICE) {
c.offset_[map] = offset;
Lattice& lat = *model::lattices[c.fill_];
offset +=
lat.fill_offset_table(offset, target_univ_id, map, univ_count_memo);
offset += lat.fill_offset_table(target_univ_id, map, univ_count_memo);
}
}
}
@ -429,32 +419,12 @@ void prepare_distribcell(const std::vector<int32_t>* user_distribcells)
//==============================================================================
void count_cell_instances(int32_t univ_indx)
void count_universe_instances()
{
const auto univ_counts = model::universe_cell_counts.find(univ_indx);
if (univ_counts != model::universe_cell_counts.end()) {
for (const auto& it : univ_counts->second) {
model::cells[it.first]->n_instances_ += it.second;
}
} else {
for (int32_t cell_indx : model::universes[univ_indx]->cells_) {
Cell& c = *model::cells[cell_indx];
++c.n_instances_;
model::universe_cell_counts[univ_indx][cell_indx] += 1;
if (c.type_ == Fill::UNIVERSE) {
// This cell contains another universe. Recurse into that universe.
count_cell_instances(c.fill_);
update_universe_cell_count(univ_indx, c.fill_);
} else if (c.type_ == Fill::LATTICE) {
// This cell contains a lattice. Recurse into the lattice universes.
Lattice& lat = *model::lattices[c.fill_];
for (auto it = lat.begin(); it != lat.end(); ++it) {
count_cell_instances(*it);
update_universe_cell_count(univ_indx, *it);
}
}
}
for (auto& univ : model::universes) {
std::unordered_map<int32_t, int32_t> univ_count_memo;
univ->n_instances_ = count_universe_instances(
model::root_universe, univ->id_, univ_count_memo);
}
}
@ -528,19 +498,16 @@ std::string distribcell_path_inner(int32_t target_cell, int32_t map,
// Material cells don't contain other cells so ignore them.
if (c.type_ != Fill::MATERIAL) {
int32_t temp_offset;
if (c.type_ == Fill::UNIVERSE) {
temp_offset =
offset + c.offset_[map]; // TODO: should also apply to lattice fills?
} else {
int32_t temp_offset = offset + c.offset_[map];
if (c.type_ == Fill::LATTICE) {
Lattice& lat = *model::lattices[c.fill_];
int32_t indx = lat.universes_.size() * map + lat.begin().indx_;
temp_offset = offset + lat.offsets_[indx];
temp_offset += lat.offsets_[indx];
}
// The desired cell is the first cell that gives an offset smaller or
// equal to the target offset.
if (temp_offset <= target_offset - c.offset_[map])
if (temp_offset <= target_offset)
break;
}
}
@ -570,12 +537,12 @@ std::string distribcell_path_inner(int32_t target_cell, int32_t map,
path << "l" << lat.id_;
for (ReverseLatticeIter it = lat.rbegin(); it != lat.rend(); ++it) {
int32_t indx = lat.universes_.size() * map + it.indx_;
int32_t temp_offset = offset + lat.offsets_[indx];
if (temp_offset <= target_offset - c.offset_[map]) {
int32_t temp_offset = offset + lat.offsets_[indx] + c.offset_[map];
if (temp_offset <= target_offset) {
offset = temp_offset;
path << "(" << lat.index_to_string(it.indx_) << ")->";
path << distribcell_path_inner(target_cell, map, target_offset,
*model::universes[*it], offset + c.offset_[map]);
path << distribcell_path_inner(
target_cell, map, target_offset, *model::universes[*it], offset);
return path.str();
}
}

View file

@ -103,8 +103,8 @@ void Lattice::adjust_indices()
//==============================================================================
int32_t Lattice::fill_offset_table(int32_t offset, int32_t target_univ_id,
int map, std::unordered_map<int32_t, int32_t>& univ_count_memo)
int32_t Lattice::fill_offset_table(int32_t target_univ_id, int map,
std::unordered_map<int32_t, int32_t>& univ_count_memo)
{
// If the offsets have already been determined for this "map", don't bother
// recalculating all of them and just return the total offset. Note that the
@ -119,6 +119,7 @@ int32_t Lattice::fill_offset_table(int32_t offset, int32_t target_univ_id,
count_universe_instances(last_univ, target_univ_id, univ_count_memo);
}
int32_t offset = 0;
for (LatticeIter it = begin(); it != end(); ++it) {
offsets_[map * universes_.size() + it.indx_] = offset;
offset += count_universe_instances(*it, target_univ_id, univ_count_memo);

View file

@ -46,7 +46,7 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
source_region_offsets_.push_back(-1);
} else {
source_region_offsets_.push_back(base_source_regions);
base_source_regions += c->n_instances_;
base_source_regions += c->n_instances();
}
}
@ -61,7 +61,7 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
for (int i = 0; i < model::cells.size(); i++) {
Cell& cell = *model::cells[i];
if (cell.type_ == Fill::MATERIAL) {
for (int j = 0; j < cell.n_instances_; j++) {
for (int j = 0; j < cell.n_instances(); j++) {
source_regions_.material(source_region_id++) = cell.material(j);
}
}
@ -1014,7 +1014,7 @@ void FlatSourceDomain::apply_external_source_to_cell_and_children(
Cell& cell = *model::cells[i_cell];
if (cell.type_ == Fill::MATERIAL) {
vector<int> instances(cell.n_instances_);
vector<int> instances(cell.n_instances());
std::iota(instances.begin(), instances.end(), 0);
apply_external_source_to_cell_instances(
i_cell, discrete, strength_factor, target_material_id, instances);
@ -1343,12 +1343,12 @@ void FlatSourceDomain::apply_mesh_to_cell_and_children(int32_t i_cell,
Cell& cell = *model::cells[i_cell];
if (cell.type_ == Fill::MATERIAL) {
vector<int> instances(cell.n_instances_);
vector<int> instances(cell.n_instances());
std::iota(instances.begin(), instances.end(), 0);
apply_mesh_to_cell_instances(
i_cell, mesh_idx, target_material_id, instances, is_target_void);
} else if (target_material_id == C_NONE && !is_target_void) {
for (int j = 0; j < cell.n_instances_; j++) {
for (int j = 0; j < cell.n_instances(); j++) {
std::unordered_map<int32_t, vector<int32_t>> cell_instance_list =
cell.get_contained_cells(j, nullptr);
for (const auto& pair : cell_instance_list) {

View file

@ -34,7 +34,7 @@ void DistribcellFilter::set_cell(int32_t cell)
assert(cell >= 0);
assert(cell < model::cells.size());
cell_ = cell;
n_bins_ = model::cells[cell]->n_instances_;
n_bins_ = model::cells[cell]->n_instances();
}
void DistribcellFilter::get_all_bins(

View file

@ -0,0 +1,74 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material depletable="true" id="1" name="UO2">
<density units="g/cm3" value="10.0"/>
<nuclide ao="1.0" name="U235"/>
<nuclide ao="2.0" name="O16"/>
</material>
<material id="2" name="light water">
<density units="g/cm3" value="1.0"/>
<nuclide ao="2.0" name="H1"/>
<nuclide ao="1.0" name="O16"/>
<sab name="c_H_in_H2O"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-1" universe="1"/>
<cell id="2" material="2" region="1" universe="1"/>
<cell fill="2" id="3" region="2 -3 4 -5" universe="6"/>
<cell fill="3" id="4" region="2 -3 5 -9" universe="6"/>
<cell fill="4" id="5" region="3 -11 4 -5" universe="6"/>
<cell fill="5" id="6" region="3 -11 5 -9" universe="6"/>
<cell fill="6" id="7" region="18 -19 20 -21" universe="7"/>
<lattice id="2">
<pitch>1.0 1.0</pitch>
<outer>1</outer>
<dimension>1 1</dimension>
<lower_left>-1.0 -1.0</lower_left>
<universes>
1 </universes>
</lattice>
<lattice id="3">
<pitch>1.0 1.0</pitch>
<outer>1</outer>
<dimension>1 1</dimension>
<lower_left>-1.0 0</lower_left>
<universes>
1 </universes>
</lattice>
<lattice id="4">
<pitch>1.0 1.0</pitch>
<outer>1</outer>
<dimension>1 1</dimension>
<lower_left>0 -1.0</lower_left>
<universes>
1 </universes>
</lattice>
<lattice id="5">
<pitch>1.0 1.0</pitch>
<outer>1</outer>
<dimension>1 1</dimension>
<lower_left>0 0</lower_left>
<universes>
1 </universes>
</lattice>
<surface coeffs="0.0 0.0 0.4" id="1" type="z-cylinder"/>
<surface coeffs="-1.0" id="2" name="minimum x" type="x-plane"/>
<surface coeffs="0.0" id="3" name="maximum x" type="x-plane"/>
<surface coeffs="-1.0" id="4" name="minimum y" type="y-plane"/>
<surface coeffs="0.0" id="5" name="maximum y" type="y-plane"/>
<surface coeffs="1.0" id="9" name="maximum y" type="y-plane"/>
<surface coeffs="1.0" id="11" name="maximum x" type="x-plane"/>
<surface boundary="reflective" coeffs="-1.0" id="18" name="minimum x" type="x-plane"/>
<surface boundary="reflective" coeffs="1.0" id="19" name="maximum x" type="x-plane"/>
<surface boundary="reflective" coeffs="-1.0" id="20" name="minimum y" type="y-plane"/>
<surface boundary="reflective" coeffs="1.0" id="21" name="maximum y" type="y-plane"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>10</batches>
<inactive>5</inactive>
</settings>
</model>

View file

@ -0,0 +1,2 @@
k-combined:
1.863277E+00 1.289821E-02

View file

@ -0,0 +1,107 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material depletable="true" id="1" name="UO2">
<density units="g/cm3" value="10.0"/>
<nuclide ao="1.0" name="U235"/>
<nuclide ao="2.0" name="O16"/>
</material>
<material depletable="true" id="2" name="UO2">
<density units="g/cm3" value="10.0"/>
<nuclide ao="1.0" name="U235"/>
<nuclide ao="2.0" name="O16"/>
</material>
<material depletable="true" id="3" name="UO2">
<density units="g/cm3" value="10.0"/>
<nuclide ao="1.0" name="U235"/>
<nuclide ao="2.0" name="O16"/>
</material>
<material depletable="true" id="4" name="UO2">
<density units="g/cm3" value="10.0"/>
<nuclide ao="1.0" name="U235"/>
<nuclide ao="2.0" name="O16"/>
</material>
<material id="5" name="light water">
<density units="g/cm3" value="1.0"/>
<nuclide ao="2.0" name="H1"/>
<nuclide ao="1.0" name="O16"/>
<sab name="c_H_in_H2O"/>
</material>
<material id="6" name="light water">
<density units="g/cm3" value="1.0"/>
<nuclide ao="2.0" name="H1"/>
<nuclide ao="1.0" name="O16"/>
<sab name="c_H_in_H2O"/>
</material>
<material id="7" name="light water">
<density units="g/cm3" value="1.0"/>
<nuclide ao="2.0" name="H1"/>
<nuclide ao="1.0" name="O16"/>
<sab name="c_H_in_H2O"/>
</material>
<material id="8" name="light water">
<density units="g/cm3" value="1.0"/>
<nuclide ao="2.0" name="H1"/>
<nuclide ao="1.0" name="O16"/>
<sab name="c_H_in_H2O"/>
</material>
</materials>
<geometry>
<cell id="8" material="1 2 3 4" region="-6" universe="8"/>
<cell id="9" material="5 6 7 8" region="6" universe="8"/>
<cell fill="9" id="10" region="7 -8 10 -12" universe="13"/>
<cell fill="10" id="11" region="7 -8 12 -16" universe="13"/>
<cell fill="11" id="12" region="8 -22 10 -12" universe="13"/>
<cell fill="12" id="13" region="8 -22 12 -16" universe="13"/>
<cell fill="13" id="14" region="29 -30 31 -32" universe="14"/>
<lattice id="9">
<pitch>1.0 1.0</pitch>
<outer>8</outer>
<dimension>1 1</dimension>
<lower_left>-1.0 -1.0</lower_left>
<universes>
8 </universes>
</lattice>
<lattice id="10">
<pitch>1.0 1.0</pitch>
<outer>8</outer>
<dimension>1 1</dimension>
<lower_left>-1.0 0</lower_left>
<universes>
8 </universes>
</lattice>
<lattice id="11">
<pitch>1.0 1.0</pitch>
<outer>8</outer>
<dimension>1 1</dimension>
<lower_left>0 -1.0</lower_left>
<universes>
8 </universes>
</lattice>
<lattice id="12">
<pitch>1.0 1.0</pitch>
<outer>8</outer>
<dimension>1 1</dimension>
<lower_left>0 0</lower_left>
<universes>
8 </universes>
</lattice>
<surface coeffs="0.0 0.0 0.4" id="6" type="z-cylinder"/>
<surface coeffs="-1.0" id="7" name="minimum x" type="x-plane"/>
<surface coeffs="0.0" id="8" name="maximum x" type="x-plane"/>
<surface coeffs="-1.0" id="10" name="minimum y" type="y-plane"/>
<surface coeffs="0.0" id="12" name="maximum y" type="y-plane"/>
<surface coeffs="1.0" id="16" name="maximum y" type="y-plane"/>
<surface coeffs="1.0" id="22" name="maximum x" type="x-plane"/>
<surface boundary="reflective" coeffs="-1.0" id="29" name="minimum x" type="x-plane"/>
<surface boundary="reflective" coeffs="1.0" id="30" name="maximum x" type="x-plane"/>
<surface boundary="reflective" coeffs="-1.0" id="31" name="minimum y" type="y-plane"/>
<surface boundary="reflective" coeffs="1.0" id="32" name="maximum y" type="y-plane"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>10</batches>
<inactive>5</inactive>
</settings>
</model>

View file

@ -0,0 +1,2 @@
k-combined:
1.863277E+00 1.289821E-02

View file

@ -0,0 +1,83 @@
import numpy as np
import openmc
from openmc.utility_funcs import change_directory
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
lattice00 = openmc.RectLattice()
lattice00.lower_left = (-d, -d)
lattice00.pitch = (d, d)
lattice00.outer = pin
lattice00.universes = [[pin]]
box00 = openmc.model.RectangularPrism(d, d, origin=(-d/2,-d/2))
lattice01 = openmc.RectLattice()
lattice01.lower_left = (-d, 0)
lattice01.pitch = (d, d)
lattice01.outer = pin
lattice01.universes = [[pin]]
box01 = openmc.model.RectangularPrism(d, d, origin=(-d/2,d/2))
lattice10 = openmc.RectLattice()
lattice10.lower_left = (0, -d)
lattice10.pitch = (d, d)
lattice10.outer = pin
lattice10.universes = [[pin]]
box10 = openmc.model.RectangularPrism(d, d, origin=(d/2,-d/2))
lattice11 = openmc.RectLattice()
lattice11.lower_left = (0, 0)
lattice11.pitch = (d, d)
lattice11.outer = pin
lattice11.universes = [[pin]]
box11 = openmc.model.RectangularPrism(d, d, origin=(d/2,d/2))
cell00 = openmc.Cell(fill=lattice00, region = -box00)
cell01 = openmc.Cell(fill=lattice01, region = -box01)
cell10 = openmc.Cell(fill=lattice10, region = -box10)
cell11 = openmc.Cell(fill=lattice11, region = -box11)
univ = openmc.Universe(cells=[cell00, cell01, cell10, cell11])
box = openmc.model.RectangularPrism(2*d, 2*d, boundary_type='reflective')
main_cell = openmc.Cell(fill=univ, region=-box)
model.geometry = openmc.Geometry([main_cell])
model.geometry.merge_surfaces = True
model.settings.batches = 10
model.settings.inactive = 5
model.settings.particles = 1000
return model
@pytest.mark.parametrize("distribmat", [False, True])
def test_lattice(model, distribmat):
with change_directory(str(distribmat)):
openmc.reset_auto_ids()
if distribmat:
model.differentiate_mats(depletable_only=False)
harness = PyAPITestHarness('statepoint.10.h5', model)
harness.main()