From 282ba43eb5d41b6155a538c0155c9779cfceeb63 Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Mon, 15 Dec 2025 10:58:47 -0600 Subject: [PATCH 01/18] First prototype for temperature field --- CMakeLists.txt | 1 + include/openmc/field.h | 30 ++++++++++++ include/openmc/mesh.h | 16 +++++++ include/openmc/particle.h | 3 +- include/openmc/settings.h | 1 + include/openmc/simulation.h | 3 ++ openmc/__init__.py | 1 + openmc/field.py | 58 +++++++++++++++++++++++ openmc/settings.py | 49 +++++++++++++++++++ src/field.cpp | 50 ++++++++++++++++++++ src/finalize.cpp | 3 ++ src/geometry_aux.cpp | 13 ++++++ src/mesh.cpp | 52 +++++++++++++++++++++ src/particle.cpp | 93 ++++++++++++++++++++++++++++++++++--- src/settings.cpp | 39 ++++++++++++++++ src/simulation.cpp | 50 ++++++++++++++++++-- 16 files changed, 450 insertions(+), 12 deletions(-) create mode 100644 include/openmc/field.h create mode 100644 openmc/field.py create mode 100644 src/field.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 87b8789d101..7f53af98d0d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -350,6 +350,7 @@ list(APPEND libopenmc_SOURCES src/endf.cpp src/error.cpp src/event.cpp + src/field.cpp src/file_utils.cpp src/finalize.cpp src/geometry.cpp diff --git a/include/openmc/field.h b/include/openmc/field.h new file mode 100644 index 00000000000..ad64c419f08 --- /dev/null +++ b/include/openmc/field.h @@ -0,0 +1,30 @@ +#ifndef OPENMC_FIELD_H +#define OPENMC_FIELD_H + +#include "openmc/vector.h" +#include "openmc/mesh.h" + +namespace openmc { + +class TemperatureField{ + +public: + Mesh* mesh_ptr; + vector values; + + TemperatureField(){}; + + TemperatureField(Mesh* mesh_ptr, vector values); + + double distance_to_next_cell(Position r, Direction d); + + double get_temperature(Position r); + + double get_sqrtkT(Position r); + + // GeometryState, particle attribute: temperature mesh bin + +}; + +} // namespace openmc +#endif // OPENMC_FIELD_H diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 5c9272e93b9..73722b2d4af 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -173,6 +173,18 @@ class Mesh { virtual void surface_bins_crossed( Position r0, Position r1, const Direction& u, vector& bins) const = 0; + //! Distance to the next boundary + //! If the initial position is outside the mesh, the distance + //! will be from the initial position to the external boundary + //! of the mesh if hit. If the initial position is inside the + //! mesh, the distance will be from the initial position to + //! the boundary of the next cell. + // + //! \param[in] r Position of the particle + //! \param[in] u Particle direction + //! \return Distance to the next boundary + virtual double distance_to_next_boundary(Position r, Position u) const = 0; + //! Get bin at a given position in space // //! \param[in] r Position to get bin for @@ -302,6 +314,8 @@ class StructuredMesh : public Mesh { void surface_bins_crossed(Position r0, Position r1, const Direction& u, vector& bins) const override; + double distance_to_next_boundary(Position r, Position u) const override; + //! Determine which cell or surface bins were crossed by a particle // //! \param[in] r0 Previous position of the particle @@ -684,6 +698,8 @@ class UnstructuredMesh : public Mesh { UnstructuredMesh(pugi::xml_node node); UnstructuredMesh(hid_t group); + double distance_to_next_boundary(Position r, Position u) const override; + static const std::string mesh_type; virtual std::string get_mesh_type() const override; diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 0f37719b94f..2d97cf92c11 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -66,11 +66,12 @@ class Particle : public ParticleData { // Coarse-grained particle events void event_calculate_xs(); - void event_advance(); + int event_advance(); void event_cross_surface(); void event_collide(); void event_revive_from_secondary(); void event_death(); + void event_cross_temperature_mesh(); //! pulse-height recording void pht_collision_energy(); diff --git a/include/openmc/settings.h b/include/openmc/settings.h index b369c99fef8..6f9979d9e42 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -70,6 +70,7 @@ extern bool extern "C" bool entropy_on; //!< calculate Shannon entropy? extern "C" bool event_based; //!< use event-based mode (instead of history-based) +extern bool temperature_field_on ; //!< Is there a temperature field defined? extern bool ifp_on; //!< Use IFP for kinetics parameters? extern bool legendre_to_tabular; //!< convert Legendre distributions to tabular? extern bool material_cell_offsets; //!< create material cells offsets? diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 9a6cf1b2131..0bfe2e40308 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -4,6 +4,7 @@ #ifndef OPENMC_SIMULATION_H #define OPENMC_SIMULATION_H +#include "openmc/field.h" #include "openmc/mesh.h" #include "openmc/particle.h" #include "openmc/vector.h" @@ -46,6 +47,8 @@ extern int64_t work_per_rank; //!< number of particles per MPI rank extern const RegularMesh* entropy_mesh; extern const RegularMesh* ufs_mesh; +extern TemperatureField temperature_field; + extern vector k_generation; extern vector work_index; diff --git a/openmc/__init__.py b/openmc/__init__.py index bb972b4e6ad..85cfdc5fb3b 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -19,6 +19,7 @@ from openmc.source import * from openmc.settings import * from openmc.lattice import * +from openmc.field import * from openmc.filter import * from openmc.filter_expansion import * from openmc.trigger import * diff --git a/openmc/field.py b/openmc/field.py new file mode 100644 index 00000000000..eb5fc951082 --- /dev/null +++ b/openmc/field.py @@ -0,0 +1,58 @@ +from abc import ABC, abstractmethod + +import openmc + +from .mixin import IDManagerMixin + + +class ScalarField(IDManagerMixin, ABC): + """ + """ + def to_xml_element(self): + """Return XML representation of the field + + Returns + ------- + element : lxml.etree._Element + XML element containing mesh data + + """ + elem = ET.Element("field") + + elem.set("id", str(self._id)) + if self.name: + elem.set("name", self.name) + + return elem + + +class TemperatureField(ScalarField): + """ + """ + def __init__(self, mesh, values): + """ + """ + self.mesh = mesh + self.values = values + + @classmethod + def from_mesh_and_values(cls, mesh, values): + """ + """ + return cls(mesh, values) + + def to_xml_element(self): + """Return XML representation of the mesh + + Returns + ------- + element : lxml.etree._Element + XML element containing mesh data + + """ + element = super().to_xml_element() + + subelement = ET.SubElement(element, "mesh") + subelement.text = self.mesh.to_xml_element() + + return element diff --git a/openmc/settings.py b/openmc/settings.py index 43c1fe06982..ab121980f54 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -14,6 +14,7 @@ from ._xml import clean_indentation, get_elem_list, get_text from .mesh import _read_meshes, RegularMesh, MeshBase from .source import SourceBase, MeshSource, IndependentSource +from .field import TemperatureField from .utility_funcs import input_path from .volume import VolumeCalculation from .weight_windows import WeightWindows, WeightWindowGenerator, WeightWindowsList @@ -400,6 +401,9 @@ def __init__(self, **kwargs): # Shannon entropy mesh self._entropy_mesh = None + # Temperature field + self._temperature_field = None + # Trigger subelement self._trigger_active = None self._trigger_max_batches = None @@ -712,6 +716,15 @@ def entropy_mesh(self) -> RegularMesh: def entropy_mesh(self, entropy: RegularMesh): cv.check_type('entropy mesh', entropy, RegularMesh) self._entropy_mesh = entropy + + @property + def temperature_field(self) -> TemperatureField: + return self._temperature_field + + @temperature_field.setter + def temperature_field(self, temperature_field: TemperatureField): + cv.check_type('temperature field', temperature_field, TemperatureField) + self._temperature_field = temperature_field @property def trigger_active(self) -> bool: @@ -1650,6 +1663,31 @@ def _create_entropy_mesh_subelement(self, root, mesh_memo=None): if mesh_memo is not None: mesh_memo.add(self.entropy_mesh.id) + def _create_temperature_field_subelement(self, root, mesh_memo=None): + if self.temperature_field is None: + return + + # add mesh ID to this element + element = ET.SubElement(root, "temperature_field") + subelement = ET.SubElement(element, "mesh") + subelement.text = str(self.temperature_field.mesh.id) + + # If this mesh has already been written outside the + # settings element, skip writing it again + if mesh_memo and self.temperature_field.mesh.id in mesh_memo: + return + + # See if a element already exists -- if not, add it + path = f"./mesh[@id='{self.temperature_field.mesh.id}']" + if root.find(path) is None: + root.append(self.temperature_field.mesh.to_xml_element()) + if mesh_memo is not None: + mesh_memo.add(self.temperature_field.mesh.id) + + # Add temperature values + subelement = ET.SubElement(element, "values") + subelement.text = " ".join([str(i) for i in self.temperature_field.values]) + def _create_trigger_subelement(self, root): if self._trigger_active is not None: trigger_element = ET.SubElement(root, "trigger") @@ -2135,6 +2173,15 @@ def _entropy_mesh_from_xml_element(self, root, meshes): raise ValueError(f'Could not locate mesh with ID "{mesh_id}"') self.entropy_mesh = meshes[mesh_id] + def _temperature_field_from_xml_element(self, root, meshes): + text = get_text(root, 'temperature_field') + if text is None: + return + mesh_id = int(text) + if mesh_id not in meshes: + raise ValueError(f'Could not locate mesh with ID "{mesh_id}"') + self.temperature_field.mesh = meshes[mesh_id] + def _trigger_from_xml_element(self, root): elem = root.find('trigger') if elem is not None: @@ -2413,6 +2460,7 @@ def to_xml_element(self, mesh_memo=None): self._create_survival_biasing_subelement(element) self._create_cutoff_subelement(element) self._create_entropy_mesh_subelement(element, mesh_memo) + self._create_temperature_field_subelement(element, mesh_memo) self._create_trigger_subelement(element) self._create_no_reduce_subelement(element) self._create_verbosity_subelement(element) @@ -2527,6 +2575,7 @@ def from_xml_element(cls, elem, meshes=None): settings._survival_biasing_from_xml_element(elem) settings._cutoff_from_xml_element(elem) settings._entropy_mesh_from_xml_element(elem, meshes) + settings._temperature_field_from_xml_element(elem, meshes) settings._trigger_from_xml_element(elem) settings._no_reduce_from_xml_element(elem) settings._verbosity_from_xml_element(elem) diff --git a/src/field.cpp b/src/field.cpp new file mode 100644 index 00000000000..49080c137c7 --- /dev/null +++ b/src/field.cpp @@ -0,0 +1,50 @@ +#include "openmc/field.h" +#include "openmc/mesh.h" +#include "openmc/vector.h" + +namespace openmc { + +TemperatureField::TemperatureField(Mesh* mesh_ptr, vector values) +{ + this->mesh_ptr = mesh_ptr; + for (double v: values) { + this->values.push_back(v); + } +} + +double TemperatureField::distance_to_next_cell(Position r, Direction u) +{ + if (this->mesh_ptr != nullptr) { + return this->mesh_ptr->distance_to_next_boundary(r, u); + } + fatal_error("TODO"); +} + +double TemperatureField::get_temperature(Position r) +{ + // Get bin from position + int i = this->mesh_ptr->get_bin(r); + + // If we have a bin then use it to locate the value + if (i >= 0 && i < this->values.size()) { + // Return the temperature + return this->values[i]; + } + + // Return -1.0 if no values were found (probably outside of the mesh) + return -1.0; +} + +double TemperatureField::get_sqrtkT(Position r) +{ + double temperature = this->get_temperature(r); + if (temperature >= 0) { + return sqrt(temperature * K_BOLTZMANN); + } else { + //TODO + fatal_error(""); + return temperature; + } +} + +} // namespace openmc diff --git a/src/finalize.cpp b/src/finalize.cpp index 9ee94340997..b9c000528d7 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -131,6 +131,7 @@ int openmc_finalize() settings::ssw_max_particles = 0; settings::ssw_max_files = 1; settings::survival_biasing = false; + settings::temperature_field_on = false; settings::temperature_default = 293.6; settings::temperature_method = TemperatureMethod::NEAREST; settings::temperature_multipole = false; @@ -158,6 +159,8 @@ int openmc_finalize() simulation::entropy_mesh = nullptr; simulation::ufs_mesh = nullptr; + simulation::temperature_field; + data::energy_max = {INFTY, INFTY, INFTY, INFTY}; data::energy_min = {0.0, 0.0, 0.0, 0.0}; data::temperature_min = 0.0; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 8a145fb1f1e..1c07ff7c745 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -18,6 +18,7 @@ #include "openmc/material.h" #include "openmc/settings.h" #include "openmc/surface.h" +#include "openmc/simulation.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/filter_cell_instance.h" #include "openmc/tallies/filter_distribcell.h" @@ -261,6 +262,18 @@ void get_temperatures( } } } + + // Add temperature if temperature field: + if (settings::temperature_field_on) { + for (auto t: simulation::temperature_field.values) { + for (size_t i = 0; i < nuc_temps.size() ; i++){ + if (!contains(nuc_temps[i], t)) { + nuc_temps[i].push_back(t); + } + } + } + } + // TODO: thermal scattering data } //============================================================================== diff --git a/src/mesh.cpp b/src/mesh.cpp index 610d057cf08..b51a6e2484c 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -720,6 +720,12 @@ UnstructuredMesh::UnstructuredMesh(hid_t group) : Mesh(group) } } +double UnstructuredMesh::distance_to_next_boundary(Position r, Position u) const +{ + fatal_error("Not implemented"); + return -1.0; +} + void UnstructuredMesh::determine_bounds() { double xmin = INFTY; @@ -1175,6 +1181,52 @@ void StructuredMesh::surface_bins_crossed( raytrace_mesh(r0, r1, u, SurfaceAggregator(this, bins)); } +double StructuredMesh::distance_to_next_boundary(Position r, Position u) const +{ + double distance = -1.0; + + // Algo adapted from mesh.cpp + Position global_r = r; + Position local_r = local_coords(r); + + const int n = 3; + + bool in_mesh; + + StructuredMesh::MeshIndex ijk = get_indices(global_r + TINY_BIT * u, in_mesh); + + // Calculate initial distances to next surfaces in all three dimensions + std::array distances; + for (int k = 0; k < n; ++k) { + distances[k] = distance_to_grid_boundary(ijk, k, local_r, u, 0.0); + } + + if (in_mesh) { + + // find surface with minimal distance to current position + const auto k = std::min_element(distances.begin(), distances.end()) - + distances.begin(); + + distance = distances[k].distance; + + } else { // not inside mesh + + // For all directions outside the mesh, find the distance that we need + // to travel to reach the next surface. Use the largest distance, as + // only this will cross all outer surfaces. + int k_max {-1}; + for (int k = 0; k < n; ++k) { + if ((ijk[k] < 1 || ijk[k] > shape_[k]) && + (distances[k].distance > distance)) { + distance = distances[k].distance; + k_max = k; + } + } + } + + return distance; +} + //============================================================================== // RegularMesh implementation //============================================================================== diff --git a/src/particle.cpp b/src/particle.cpp index 2d70d715e22..811110b1e12 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -193,6 +193,18 @@ void Particle::event_calculate_xs() cell_last(j) = coord(j).cell(); } n_coord_last() = n_coord(); + + // Update temperature of the particle if temperature field + if (settings::temperature_field_on) { + double temperature_field_sqrtkT = simulation::temperature_field.get_sqrtkT(r()); + // Update temperature if we are inside the mesh + if (temperature_field_sqrtkT >= 0.) { + sqrtkT() = temperature_field_sqrtkT; + } else { + //TODO + fatal_error(""); + } + } } // Write particle track. @@ -229,9 +241,11 @@ void Particle::event_calculate_xs() } } -void Particle::event_advance() +int Particle::event_advance() { - // Find the distance to the nearest boundary + int stop_condition = 0; + + // Find the distance to the nearest geometry boundary boundary() = distance_to_boundary(*this); // Sample a distance to collision @@ -243,14 +257,32 @@ void Particle::event_advance() collision_distance() = -std::log(prn(current_seed())) / macro_xs().total; } + // Find the distance to the nearest temperature mesh cell surface + double distance_tmesh = INFTY; + if (settings::temperature_field_on) { + distance_tmesh = simulation::temperature_field.distance_to_next_cell(r(), u()); + } + + // Calculate the distance corresponding to the time cutoff double speed = this->speed(); double time_cutoff = settings::time_cutoff[static_cast(type())]; double distance_cutoff = (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY; - // Select smaller of the three distances + // Select smaller of the four distances double distance = - std::min({boundary().distance(), collision_distance(), distance_cutoff}); + std::min({boundary().distance(), collision_distance(), distance_cutoff, distance_tmesh}); + + // Prepare the stop condition + if (distance == distance_cutoff) { + stop_condition = 4; + } else if (distance == boundary().distance()) { + stop_condition = 1; + } else if (distance == distance_tmesh) { + stop_condition = 3; + } else if (distance == collision_distance()) { + stop_condition = 2; + } // Advance particle in space and time this->move_distance(distance); @@ -279,10 +311,44 @@ void Particle::event_advance() score_track_derivative(*this, distance); } + //if (p.collision_distance() > p.boundary().distance()) { + // p.event_cross_surface(); + + + + //const int CROSS_SURFACE = 1; + //const int COLLIDE = 2; + //const int CROSS_TEMPERATURE_MESH = 3; + //const int CUTOFF = 4; + // Set particle weight to zero if it hit the time boundary - if (distance == distance_cutoff) { - wgt() = 0.0; + //if (distance == distance_cutoff) { + // wgt() = 0.0; + //} + + return stop_condition; +} + +void Particle::event_cross_temperature_mesh() +{ + // Update temperature of the particle + sqrtkT_last() = sqrtkT(); + + // Find the temperture from the mesh + double temperature_field_sqrtkT = simulation::temperature_field.get_sqrtkT(r() + u() * TINY_BIT); + + // Update temperature + // If inside the mesh + if (temperature_field_sqrtkT >= 0.) { + sqrtkT() = temperature_field_sqrtkT; + // If outside the mesh + } else { + //TODO + fatal_error(""); } + + // TODO: what do we do if we are outside the mesh? + // Probably need to use the cell instance temperature } void Particle::event_cross_surface() @@ -322,6 +388,21 @@ void Particle::event_cross_surface() } event() = TallyEvent::SURFACE; } + + // Update temperature of the particle if temperature field + // This is important here just on case we both crossed a cell + // from the geometry and from the temperature mesh + if (settings::temperature_field_on) { + double temperature_field_sqrtkT = simulation::temperature_field.get_sqrtkT(r() + u() * TINY_BIT); + // Update temperature if we are inside the mesh + if (temperature_field_sqrtkT >= 0.) { + sqrtkT() = temperature_field_sqrtkT; + } else { + //TODO + fatal_error(""); + } + } + // Score cell to cell partial currents if (!model::active_surface_tallies.empty()) { score_surface_tally(*this, model::active_surface_tallies); diff --git a/src/settings.cpp b/src/settings.cpp index 9dcf7c8dbc8..82a408d942e 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -19,6 +19,7 @@ #include "openmc/distribution_spatial.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" +#include "openmc/field.h" #include "openmc/file_utils.h" #include "openmc/mcpl_interface.h" #include "openmc/mesh.h" @@ -75,6 +76,7 @@ bool surf_mcpl_write {false}; bool surf_source_read {false}; bool survival_biasing {false}; bool survival_normalization {false}; +bool temperature_field_on {false}; bool temperature_multipole {false}; bool trigger_on {false}; bool trigger_predict {false}; @@ -779,6 +781,43 @@ void read_settings_xml(pugi::xml_node root) "it by specifying its ID in an element."); } } + + // Temperature field + if (check_for_node(root, "temperature_field")) { + temperature_field_on = true; + + // Get pointer to temperature_field node + auto node_tf = root.child("temperature_field"); + + // Mesh parameter + Mesh* tf_mesh_ptr; + if (check_for_node(node_tf, "mesh")) { + int temp = std::stoi(get_node_value(node_tf, "mesh")); + if (model::mesh_map.find(temp) == model::mesh_map.end()) { + fatal_error(fmt::format( + "Mesh {} specified for the temperature field does not exist.", temp)); + } + tf_mesh_ptr = model::meshes[model::mesh_map.at(temp)].get(); + } else { + fatal_error("A mesh should be given for the temperature field."); + } + + // Values parameter + vector tf_values; + if (check_for_node(node_tf, "values")) { + auto temp = get_node_array(node_tf, "values"); + for (const auto& b : temp) { + tf_values.push_back(b); + } + } else { + fatal_error("Temperature values should be given for the temperature field."); + } + + // Instantiate the temperature field + //TemperatureField tf = TemperatureField(tf_mesh_ptr, tf_values); + simulation::temperature_field = TemperatureField(tf_mesh_ptr, tf_values); + } + // Uniform fission source weighting mesh if (check_for_node(root, "ufs_mesh")) { auto temp = std::stoi(get_node_value(root, "ufs_mesh")); diff --git a/src/simulation.cpp b/src/simulation.cpp index b536ae5881f..7675cf04657 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -7,6 +7,7 @@ #include "openmc/eigenvalue.h" #include "openmc/error.h" #include "openmc/event.h" +#include "openmc/field.h" #include "openmc/geometry_aux.h" #include "openmc/ifp.h" #include "openmc/material.h" @@ -320,6 +321,8 @@ int64_t work_per_rank; const RegularMesh* entropy_mesh {nullptr}; const RegularMesh* ufs_mesh {nullptr}; +TemperatureField temperature_field; + vector k_generation; vector work_index; @@ -801,18 +804,55 @@ void free_memory_simulation() simulation::entropy.clear(); } +//void transport_history_based_single_particle(Particle& p) +//{ +// while (p.alive()) { +// p.event_calculate_xs(); +// if (p.alive()) { +// p.event_advance(); +// } +// if (p.alive()) { +// if (p.collision_distance() > p.boundary().distance()) { +// p.event_cross_surface(); +// } else if (p.alive()) { +// p.event_collide(); +// } +// } +// p.event_revive_from_secondary(); +// } +// p.event_death(); +//} + +const int CROSS_SURFACE = 1; +const int COLLIDE = 2; +const int CROSS_TEMPERATURE_MESH = 3; +const int CUTOFF = 4; + void transport_history_based_single_particle(Particle& p) { + int advance_stop_condition; while (p.alive()) { p.event_calculate_xs(); if (p.alive()) { - p.event_advance(); + advance_stop_condition = p.event_advance(); // TODO maybe use a particle attribute instead (next_event) } if (p.alive()) { - if (p.collision_distance() > p.boundary().distance()) { - p.event_cross_surface(); - } else if (p.alive()) { - p.event_collide(); + switch (advance_stop_condition) { + case CROSS_SURFACE: + p.event_cross_surface(); + break; + case COLLIDE: + p.event_collide(); + break; + case CROSS_TEMPERATURE_MESH: + p.event_cross_temperature_mesh(); + break; + case CUTOFF: + p.wgt() = 0.0; + break; + default: + fatal_error("Unknown case"); + break; } } p.event_revive_from_secondary(); From cf1970937c0c821a49a9e7258c78df1231aeaf4f Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Mon, 15 Dec 2025 11:17:08 -0600 Subject: [PATCH 02/18] Clean event identifiers for history-based transport --- include/openmc/constants.h | 8 ++++++++ src/particle.cpp | 23 ++++------------------- src/simulation.cpp | 34 +++++----------------------------- 3 files changed, 17 insertions(+), 48 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index bba260c3a4b..15d6bf122a7 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -376,6 +376,14 @@ enum class GeometryType { CSG, DAG }; // representations. This value represents no surface. constexpr int32_t SURFACE_NONE {0}; +//============================================================================== +// EVENT IDENTIFIER IN HISTORY-BASED TRANSPORT + +const int EVENT_CROSS_SURFACE = 1; +const int EVENT_COLLIDE = 2; +const int EVENT_CROSS_TEMPERATURE_MESH = 3; +const int EVENT_TIME_CUTOFF = 4; + } // namespace openmc #endif // OPENMC_CONSTANTS_H diff --git a/src/particle.cpp b/src/particle.cpp index 811110b1e12..7d7fa3bbd38 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -275,13 +275,13 @@ int Particle::event_advance() // Prepare the stop condition if (distance == distance_cutoff) { - stop_condition = 4; + stop_condition = EVENT_TIME_CUTOFF; } else if (distance == boundary().distance()) { - stop_condition = 1; + stop_condition = EVENT_CROSS_SURFACE; } else if (distance == distance_tmesh) { - stop_condition = 3; + stop_condition = EVENT_CROSS_TEMPERATURE_MESH; } else if (distance == collision_distance()) { - stop_condition = 2; + stop_condition = EVENT_COLLIDE; } // Advance particle in space and time @@ -311,21 +311,6 @@ int Particle::event_advance() score_track_derivative(*this, distance); } - //if (p.collision_distance() > p.boundary().distance()) { - // p.event_cross_surface(); - - - - //const int CROSS_SURFACE = 1; - //const int COLLIDE = 2; - //const int CROSS_TEMPERATURE_MESH = 3; - //const int CUTOFF = 4; - - // Set particle weight to zero if it hit the time boundary - //if (distance == distance_cutoff) { - // wgt() = 0.0; - //} - return stop_condition; } diff --git a/src/simulation.cpp b/src/simulation.cpp index 7675cf04657..7201f900ec6 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -804,30 +804,6 @@ void free_memory_simulation() simulation::entropy.clear(); } -//void transport_history_based_single_particle(Particle& p) -//{ -// while (p.alive()) { -// p.event_calculate_xs(); -// if (p.alive()) { -// p.event_advance(); -// } -// if (p.alive()) { -// if (p.collision_distance() > p.boundary().distance()) { -// p.event_cross_surface(); -// } else if (p.alive()) { -// p.event_collide(); -// } -// } -// p.event_revive_from_secondary(); -// } -// p.event_death(); -//} - -const int CROSS_SURFACE = 1; -const int COLLIDE = 2; -const int CROSS_TEMPERATURE_MESH = 3; -const int CUTOFF = 4; - void transport_history_based_single_particle(Particle& p) { int advance_stop_condition; @@ -838,20 +814,20 @@ void transport_history_based_single_particle(Particle& p) } if (p.alive()) { switch (advance_stop_condition) { - case CROSS_SURFACE: + case EVENT_CROSS_SURFACE: p.event_cross_surface(); break; - case COLLIDE: + case EVENT_COLLIDE: p.event_collide(); break; - case CROSS_TEMPERATURE_MESH: + case EVENT_CROSS_TEMPERATURE_MESH: p.event_cross_temperature_mesh(); break; - case CUTOFF: + case EVENT_TIME_CUTOFF: p.wgt() = 0.0; break; default: - fatal_error("Unknown case"); + fatal_error("Unknown event in history-based transport!"); break; } } From 032275a90de12221e14cf5e4064a57e93436bedc Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Mon, 15 Dec 2025 11:43:39 -0600 Subject: [PATCH 03/18] New particle attribute for the next event --- include/openmc/particle.h | 4 ++-- include/openmc/particle_data.h | 7 +++++++ src/particle.cpp | 14 +++++--------- src/simulation.cpp | 10 +++++----- 4 files changed, 19 insertions(+), 16 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 2d97cf92c11..8daf8290b12 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -66,12 +66,12 @@ class Particle : public ParticleData { // Coarse-grained particle events void event_calculate_xs(); - int event_advance(); + void event_advance(); void event_cross_surface(); + void event_cross_temperature_mesh(); void event_collide(); void event_revive_from_secondary(); void event_death(); - void event_cross_temperature_mesh(); //! pulse-height recording void pht_collision_energy(); diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index fdacfa765b3..829d503d0fa 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -565,6 +565,8 @@ class ParticleData : public GeometryState { int64_t n_progeny_ {0}; + int next_event_ {0}; + public: //---------------------------------------------------------------------------- // Constructors @@ -768,6 +770,11 @@ class ParticleData : public GeometryState { d = 0; } } + + // Next event identifier for history-based transport + int next_event() const { return next_event_; } + int& next_event() { return next_event_; } + }; } // namespace openmc diff --git a/src/particle.cpp b/src/particle.cpp index 7d7fa3bbd38..c465cae3526 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -241,10 +241,8 @@ void Particle::event_calculate_xs() } } -int Particle::event_advance() +void Particle::event_advance() { - int stop_condition = 0; - // Find the distance to the nearest geometry boundary boundary() = distance_to_boundary(*this); @@ -275,13 +273,13 @@ int Particle::event_advance() // Prepare the stop condition if (distance == distance_cutoff) { - stop_condition = EVENT_TIME_CUTOFF; + next_event() = EVENT_TIME_CUTOFF; } else if (distance == boundary().distance()) { - stop_condition = EVENT_CROSS_SURFACE; + next_event() = EVENT_CROSS_SURFACE; } else if (distance == distance_tmesh) { - stop_condition = EVENT_CROSS_TEMPERATURE_MESH; + next_event() = EVENT_CROSS_TEMPERATURE_MESH; } else if (distance == collision_distance()) { - stop_condition = EVENT_COLLIDE; + next_event() = EVENT_COLLIDE; } // Advance particle in space and time @@ -310,8 +308,6 @@ int Particle::event_advance() if (!model::active_tallies.empty()) { score_track_derivative(*this, distance); } - - return stop_condition; } void Particle::event_cross_temperature_mesh() diff --git a/src/simulation.cpp b/src/simulation.cpp index 7201f900ec6..7a1af7c2afe 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -806,14 +806,13 @@ void free_memory_simulation() void transport_history_based_single_particle(Particle& p) { - int advance_stop_condition; while (p.alive()) { p.event_calculate_xs(); if (p.alive()) { - advance_stop_condition = p.event_advance(); // TODO maybe use a particle attribute instead (next_event) + p.event_advance(); } if (p.alive()) { - switch (advance_stop_condition) { + switch (p.next_event()) { case EVENT_CROSS_SURFACE: p.event_cross_surface(); break; @@ -827,9 +826,10 @@ void transport_history_based_single_particle(Particle& p) p.wgt() = 0.0; break; default: - fatal_error("Unknown event in history-based transport!"); + fatal_error(fmt::format( + "Unknown event '{}' in history-based transport!", p.next_event())); break; - } + } } p.event_revive_from_secondary(); } From 8a138dae5332642cd2bf5fa828e3e819967ce7ef Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Mon, 15 Dec 2025 14:30:37 -0600 Subject: [PATCH 04/18] Cleaning --- include/openmc/field.h | 62 ++++++++++++++++++++++++++++++------------ src/field.cpp | 46 ++++++++++++++++++------------- src/geometry_aux.cpp | 2 +- src/particle.cpp | 39 ++++++++++++-------------- 4 files changed, 91 insertions(+), 58 deletions(-) diff --git a/include/openmc/field.h b/include/openmc/field.h index ad64c419f08..8154ec4c6d2 100644 --- a/include/openmc/field.h +++ b/include/openmc/field.h @@ -1,29 +1,57 @@ #ifndef OPENMC_FIELD_H #define OPENMC_FIELD_H -#include "openmc/vector.h" #include "openmc/mesh.h" +#include "openmc/vector.h" namespace openmc { -class TemperatureField{ +class TemperatureField { public: - Mesh* mesh_ptr; - vector values; - - TemperatureField(){}; - - TemperatureField(Mesh* mesh_ptr, vector values); - - double distance_to_next_cell(Position r, Direction d); - - double get_temperature(Position r); - - double get_sqrtkT(Position r); - - // GeometryState, particle attribute: temperature mesh bin - + //---------------------------------------------------------------------------- + // Constructors + TemperatureField() {}; + TemperatureField(Mesh* mesh_ptr, vector values); + + //---------------------------------------------------------------------------- + // Methods + + //! Returns the distance to the next mesh boundary gien a particle position + //! and direction. If the particle is initially outside, the distance will + //! correspond to the nearest distance to the outer boundaries of the mesh. + // + //! \param[in] r Position of the particle + //! \param[in] u Direction of the particle + //! \return The distance in cm to the next mesh boundary + double distance_to_next_boundary(const Position& r, const Direction& d); + + //! Returns the temperature in Kelvin corresponding to the given position. + // + //! \param[in] r Position of the particle + //! \return Temperature in Kelvin + double get_temperature(const Position& r); + + //! Returns the square root of the temperature multiplied by the Boltzmann + //! constant in eV for the given position. + // + //! \param[in] r Position of the particle + //! \return Sqrt(k_Boltzmann * temperature) in eV + double get_sqrtkT(const Position& r); + + //---------------------------------------------------------------------------- + // Accessors + + // Temperature values + double& value(int i) { return values_[i]; } + const double& value(int i) const { return values_[i]; } + const vector& values() const { return values_; } + +private: + //---------------------------------------------------------------------------- + // Data members + Mesh* mesh_ptr_; //!< Pointer to the geometric mesh + vector values_; //!< Temperature values }; } // namespace openmc diff --git a/src/field.cpp b/src/field.cpp index 49080c137c7..ee32ebf477d 100644 --- a/src/field.cpp +++ b/src/field.cpp @@ -1,47 +1,55 @@ #include "openmc/field.h" +#include "openmc/constants.h" #include "openmc/mesh.h" #include "openmc/vector.h" namespace openmc { -TemperatureField::TemperatureField(Mesh* mesh_ptr, vector values) +TemperatureField::TemperatureField(Mesh* mesh_ptr, vector values) { - this->mesh_ptr = mesh_ptr; - for (double v: values) { - this->values.push_back(v); + this->mesh_ptr_ = mesh_ptr; + for (double v : values) { + this->values_.push_back(v); } } -double TemperatureField::distance_to_next_cell(Position r, Direction u) +double TemperatureField::distance_to_next_boundary( + const Position& r, const Direction& u) { - if (this->mesh_ptr != nullptr) { - return this->mesh_ptr->distance_to_next_boundary(r, u); + if (this->mesh_ptr_ != nullptr) { + return this->mesh_ptr_->distance_to_next_boundary(r, u); } fatal_error("TODO"); } -double TemperatureField::get_temperature(Position r) +double TemperatureField::get_temperature(const Position& r) { - // Get bin from position - int i = this->mesh_ptr->get_bin(r); + if (this->mesh_ptr_ != nullptr) { - // If we have a bin then use it to locate the value - if (i >= 0 && i < this->values.size()) { - // Return the temperature - return this->values[i]; + // Get bin from position + int i = this->mesh_ptr_->get_bin(r); + + // If we have a bin then use it to locate the value + if (i >= 0 && i < this->values_.size()) { + // Return the temperature + return this->values_[i]; + } + + // Return -1.0 if no values were found (probably outside of the mesh) + return -1.0; + + } else { + fatal_error("TODO"); } - - // Return -1.0 if no values were found (probably outside of the mesh) - return -1.0; } -double TemperatureField::get_sqrtkT(Position r) +double TemperatureField::get_sqrtkT(const Position& r) { double temperature = this->get_temperature(r); if (temperature >= 0) { return sqrt(temperature * K_BOLTZMANN); } else { - //TODO + // TODO fatal_error(""); return temperature; } diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 1c07ff7c745..b0e4e959624 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -265,7 +265,7 @@ void get_temperatures( // Add temperature if temperature field: if (settings::temperature_field_on) { - for (auto t: simulation::temperature_field.values) { + for (auto t: simulation::temperature_field.values()) { for (size_t i = 0; i < nuc_temps.size() ; i++){ if (!contains(nuc_temps[i], t)) { nuc_temps[i].push_back(t); diff --git a/src/particle.cpp b/src/particle.cpp index c465cae3526..5c6932e7ff1 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -196,10 +196,10 @@ void Particle::event_calculate_xs() // Update temperature of the particle if temperature field if (settings::temperature_field_on) { - double temperature_field_sqrtkT = simulation::temperature_field.get_sqrtkT(r()); + double field_sqrtkT = simulation::temperature_field.get_sqrtkT(r()); // Update temperature if we are inside the mesh - if (temperature_field_sqrtkT >= 0.) { - sqrtkT() = temperature_field_sqrtkT; + if (field_sqrtkT >= 0.) { + sqrtkT() = field_sqrtkT; } else { //TODO fatal_error(""); @@ -258,7 +258,7 @@ void Particle::event_advance() // Find the distance to the nearest temperature mesh cell surface double distance_tmesh = INFTY; if (settings::temperature_field_on) { - distance_tmesh = simulation::temperature_field.distance_to_next_cell(r(), u()); + distance_tmesh = simulation::temperature_field.distance_to_next_boundary(r(), u()); } // Calculate the distance corresponding to the time cutoff @@ -312,24 +312,21 @@ void Particle::event_advance() void Particle::event_cross_temperature_mesh() { - // Update temperature of the particle + // Save current temperature sqrtkT_last() = sqrtkT(); - // Find the temperture from the mesh - double temperature_field_sqrtkT = simulation::temperature_field.get_sqrtkT(r() + u() * TINY_BIT); - - // Update temperature - // If inside the mesh - if (temperature_field_sqrtkT >= 0.) { - sqrtkT() = temperature_field_sqrtkT; - // If outside the mesh + // Find the field temperature from the mesh + double field_sqrtkT = + simulation::temperature_field.get_sqrtkT(r() + u() * TINY_BIT); + + // Update the particle temperature + if (field_sqrtkT >= 0.) { + sqrtkT() = field_sqrtkT; } else { - //TODO - fatal_error(""); + // TODO: If we leave the temperature mesh, we need to go back to the cell + // temperature + fatal_error("Not implemented yet"); } - - // TODO: what do we do if we are outside the mesh? - // Probably need to use the cell instance temperature } void Particle::event_cross_surface() @@ -374,10 +371,10 @@ void Particle::event_cross_surface() // This is important here just on case we both crossed a cell // from the geometry and from the temperature mesh if (settings::temperature_field_on) { - double temperature_field_sqrtkT = simulation::temperature_field.get_sqrtkT(r() + u() * TINY_BIT); + double field_sqrtkT = simulation::temperature_field.get_sqrtkT(r() + u() * TINY_BIT); // Update temperature if we are inside the mesh - if (temperature_field_sqrtkT >= 0.) { - sqrtkT() = temperature_field_sqrtkT; + if (field_sqrtkT >= 0.) { + sqrtkT() = field_sqrtkT; } else { //TODO fatal_error(""); From d29e7007956db87f9858a553d7c36efb290ce8b0 Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Mon, 15 Dec 2025 15:42:44 -0600 Subject: [PATCH 05/18] Common update function for the particle temperature --- include/openmc/field.h | 8 ++++++++ src/field.cpp | 20 ++++++++++++++++++++ src/particle.cpp | 42 +++++++----------------------------------- 3 files changed, 35 insertions(+), 35 deletions(-) diff --git a/include/openmc/field.h b/include/openmc/field.h index 8154ec4c6d2..842886b3fd0 100644 --- a/include/openmc/field.h +++ b/include/openmc/field.h @@ -39,6 +39,14 @@ class TemperatureField { //! \return Sqrt(k_Boltzmann * temperature) in eV double get_sqrtkT(const Position& r); + //! Update the temperature of a particle based on its position and direction. + //! If the particle is inside the temperature field, its temperature is + //! updated. If outside, the particle takes the temperature value + //! associated with the current cell instance. + // + //! \param[inout] p Particle + void update_particle_temperature(Particle& p); + //---------------------------------------------------------------------------- // Accessors diff --git a/src/field.cpp b/src/field.cpp index ee32ebf477d..4a4b1c96698 100644 --- a/src/field.cpp +++ b/src/field.cpp @@ -1,4 +1,5 @@ #include "openmc/field.h" +#include "openmc/cell.h" #include "openmc/constants.h" #include "openmc/mesh.h" #include "openmc/vector.h" @@ -55,4 +56,23 @@ double TemperatureField::get_sqrtkT(const Position& r) } } +void TemperatureField::update_particle_temperature(Particle& p) { + + // Save current temperature + p.sqrtkT_last() = p.sqrtkT(); + + // Determine the temperature based on the temperature field + double field_sqrtkT = this->get_sqrtkT(p.r() + p.u() * TINY_BIT); + + // If particle inside the mesh, we use the temperature field + if (field_sqrtkT >= 0.) { + p.sqrtkT() = field_sqrtkT; + + // If particle outside the mesh, go back to the cell instance temperature + } else { + Cell& c {*model::cells[p.lowest_coord().cell()]}; + p.sqrtkT() = c.sqrtkT(p.cell_instance()); + } +} + } // namespace openmc diff --git a/src/particle.cpp b/src/particle.cpp index 5c6932e7ff1..0f68dba912e 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -196,14 +196,7 @@ void Particle::event_calculate_xs() // Update temperature of the particle if temperature field if (settings::temperature_field_on) { - double field_sqrtkT = simulation::temperature_field.get_sqrtkT(r()); - // Update temperature if we are inside the mesh - if (field_sqrtkT >= 0.) { - sqrtkT() = field_sqrtkT; - } else { - //TODO - fatal_error(""); - } + simulation::temperature_field.update_particle_temperature(*this); } } @@ -258,7 +251,8 @@ void Particle::event_advance() // Find the distance to the nearest temperature mesh cell surface double distance_tmesh = INFTY; if (settings::temperature_field_on) { - distance_tmesh = simulation::temperature_field.distance_to_next_boundary(r(), u()); + distance_tmesh = + simulation::temperature_field.distance_to_next_boundary(r(), u()); } // Calculate the distance corresponding to the time cutoff @@ -267,7 +261,7 @@ void Particle::event_advance() double distance_cutoff = (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY; - // Select smaller of the four distances + // Select smaller distance double distance = std::min({boundary().distance(), collision_distance(), distance_cutoff, distance_tmesh}); @@ -312,21 +306,8 @@ void Particle::event_advance() void Particle::event_cross_temperature_mesh() { - // Save current temperature - sqrtkT_last() = sqrtkT(); - - // Find the field temperature from the mesh - double field_sqrtkT = - simulation::temperature_field.get_sqrtkT(r() + u() * TINY_BIT); - - // Update the particle temperature - if (field_sqrtkT >= 0.) { - sqrtkT() = field_sqrtkT; - } else { - // TODO: If we leave the temperature mesh, we need to go back to the cell - // temperature - fatal_error("Not implemented yet"); - } + // Update temperature of the particle + simulation::temperature_field.update_particle_temperature(*this); } void Particle::event_cross_surface() @@ -368,17 +349,8 @@ void Particle::event_cross_surface() } // Update temperature of the particle if temperature field - // This is important here just on case we both crossed a cell - // from the geometry and from the temperature mesh if (settings::temperature_field_on) { - double field_sqrtkT = simulation::temperature_field.get_sqrtkT(r() + u() * TINY_BIT); - // Update temperature if we are inside the mesh - if (field_sqrtkT >= 0.) { - sqrtkT() = field_sqrtkT; - } else { - //TODO - fatal_error(""); - } + simulation::temperature_field.update_particle_temperature(*this); } // Score cell to cell partial currents From 2e5e0c907c5cd35e0a907d5b2dab063b9d2f2909 Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Mon, 15 Dec 2025 15:58:39 -0600 Subject: [PATCH 06/18] Cleaning --- src/field.cpp | 23 ++++++++++------------- src/settings.cpp | 4 ++-- src/simulation.cpp | 1 + 3 files changed, 13 insertions(+), 15 deletions(-) diff --git a/src/field.cpp b/src/field.cpp index 4a4b1c96698..dab214efcf9 100644 --- a/src/field.cpp +++ b/src/field.cpp @@ -19,8 +19,9 @@ double TemperatureField::distance_to_next_boundary( { if (this->mesh_ptr_ != nullptr) { return this->mesh_ptr_->distance_to_next_boundary(r, u); + } else { + fatal_error("No mesh found for the temperature field!"); } - fatal_error("TODO"); } double TemperatureField::get_temperature(const Position& r) @@ -30,17 +31,16 @@ double TemperatureField::get_temperature(const Position& r) // Get bin from position int i = this->mesh_ptr_->get_bin(r); - // If we have a bin then use it to locate the value + // If we have a bin, we use it to locate the value if (i >= 0 && i < this->values_.size()) { - // Return the temperature return this->values_[i]; } - // Return -1.0 if no values were found (probably outside of the mesh) + // No values were found (outside the mesh) return -1.0; } else { - fatal_error("TODO"); + fatal_error("No mesh found for the temperature field!"); } } @@ -49,26 +49,23 @@ double TemperatureField::get_sqrtkT(const Position& r) double temperature = this->get_temperature(r); if (temperature >= 0) { return sqrt(temperature * K_BOLTZMANN); - } else { - // TODO - fatal_error(""); - return temperature; } + return -1.0; } -void TemperatureField::update_particle_temperature(Particle& p) { - +void TemperatureField::update_particle_temperature(Particle& p) +{ // Save current temperature p.sqrtkT_last() = p.sqrtkT(); // Determine the temperature based on the temperature field double field_sqrtkT = this->get_sqrtkT(p.r() + p.u() * TINY_BIT); - + // If particle inside the mesh, we use the temperature field if (field_sqrtkT >= 0.) { p.sqrtkT() = field_sqrtkT; - // If particle outside the mesh, go back to the cell instance temperature + // If particle outside the mesh, go back to the cell instance temperature } else { Cell& c {*model::cells[p.lowest_coord().cell()]}; p.sqrtkT() = c.sqrtkT(p.cell_instance()); diff --git a/src/settings.cpp b/src/settings.cpp index 82a408d942e..1d4fb8d3041 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -810,11 +810,11 @@ void read_settings_xml(pugi::xml_node root) tf_values.push_back(b); } } else { - fatal_error("Temperature values should be given for the temperature field."); + fatal_error( + "Temperature values should be given for the temperature field."); } // Instantiate the temperature field - //TemperatureField tf = TemperatureField(tf_mesh_ptr, tf_values); simulation::temperature_field = TemperatureField(tf_mesh_ptr, tf_values); } diff --git a/src/simulation.cpp b/src/simulation.cpp index 7a1af7c2afe..2255522ba2d 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -3,6 +3,7 @@ #include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/collision_track.h" +#include "openmc/constants.h" #include "openmc/container_util.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" From ba02d466ce7595b60bdb9c66d30af46b3a25ae89 Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Mon, 15 Dec 2025 16:13:08 -0600 Subject: [PATCH 07/18] Cleaning --- include/openmc/mesh.h | 6 ++---- src/mesh.cpp | 17 +++++++---------- 2 files changed, 9 insertions(+), 14 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 73722b2d4af..43988d0049a 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -173,12 +173,10 @@ class Mesh { virtual void surface_bins_crossed( Position r0, Position r1, const Direction& u, vector& bins) const = 0; - //! Distance to the next boundary + //! Distance to the next boundary. //! If the initial position is outside the mesh, the distance //! will be from the initial position to the external boundary - //! of the mesh if hit. If the initial position is inside the - //! mesh, the distance will be from the initial position to - //! the boundary of the next cell. + //! of the mesh if hit. // //! \param[in] r Position of the particle //! \param[in] u Particle direction diff --git a/src/mesh.cpp b/src/mesh.cpp index b51a6e2484c..0c696eed220 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1183,14 +1183,11 @@ void StructuredMesh::surface_bins_crossed( double StructuredMesh::distance_to_next_boundary(Position r, Position u) const { - double distance = -1.0; - - // Algo adapted from mesh.cpp Position global_r = r; Position local_r = local_coords(r); - const int n = 3; - + double distance = INFTY; + const int n = n_dimension_; bool in_mesh; StructuredMesh::MeshIndex ijk = get_indices(global_r + TINY_BIT * u, in_mesh); @@ -1203,10 +1200,10 @@ double StructuredMesh::distance_to_next_boundary(Position r, Position u) const if (in_mesh) { - // find surface with minimal distance to current position - const auto k = std::min_element(distances.begin(), distances.end()) - - distances.begin(); - + // Find surface with minimal distance to current position + const auto k = + std::min_element(distances.begin(), distances.end()) - distances.begin(); + distance = distances[k].distance; } else { // not inside mesh @@ -1224,7 +1221,7 @@ double StructuredMesh::distance_to_next_boundary(Position r, Position u) const } } - return distance; + return distance; } //============================================================================== From 5c52c5fbb23d05e663842e76013d7118ab82532f Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Wed, 17 Dec 2025 15:56:00 -0600 Subject: [PATCH 08/18] Add temperature values during the initialization for thermal scattering data --- src/geometry_aux.cpp | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index b0e4e959624..bdb340b963f 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -263,17 +263,27 @@ void get_temperatures( } } - // Add temperature if temperature field: + // Add temperatures from a temperature field. + // We assume that we do not know how geometric cells are impacted by the + // temperature field in advance. If we had access to this information, we + // could limit the declarations of temperature from the temperature field to + // impacted nuclides only. if (settings::temperature_field_on) { for (auto t: simulation::temperature_field.values()) { + // Nuclide temperatures for (size_t i = 0; i < nuc_temps.size() ; i++){ if (!contains(nuc_temps[i], t)) { nuc_temps[i].push_back(t); } } + // Thermal scattering temperatures + for (size_t i = 0; i < thermal_temps.size() ; i++){ + if (!contains(thermal_temps[i], t)) { + thermal_temps[i].push_back(t); + } + } } } - // TODO: thermal scattering data } //============================================================================== From 06f27c7f81f5eed20e445903f03a38e67bd74094 Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Fri, 19 Dec 2025 15:20:54 -0600 Subject: [PATCH 09/18] Clean Python API --- openmc/field.py | 85 +++++++++++++++++++++++++--------------------- openmc/settings.py | 9 ++++- 2 files changed, 55 insertions(+), 39 deletions(-) diff --git a/openmc/field.py b/openmc/field.py index eb5fc951082..463d997e602 100644 --- a/openmc/field.py +++ b/openmc/field.py @@ -2,57 +2,66 @@ import openmc -from .mixin import IDManagerMixin +class ScalarField(ABC): + """Scalar field defined on a geometric mesh. -class ScalarField(IDManagerMixin, ABC): - """ - """ - def to_xml_element(self): - """Return XML representation of the field - - Returns - ------- - element : lxml.etree._Element - XML element containing mesh data - - """ - elem = ET.Element("field") - - elem.set("id", str(self._id)) - if self.name: - elem.set("name", self.name) + Attributes + ---------- + mesh : Mesh + Spatial mesh associated with the field + values : iterable of float + List of values associated with each mesh cell - return elem - - -class TemperatureField(ScalarField): - """ """ def __init__(self, mesh, values): - """ + """Initialization. + + Parameters + ---------- + mesh : Mesh + Spatial mesh associated with the field + values : iterable of float + List of values associated with each mesh cell + """ self.mesh = mesh self.values = values @classmethod - def from_mesh_and_values(cls, mesh, values): - """ + def from_exodus_file(cls, filepath): + """Construct a ScalarField from an Exodus mesh file + + Parameters + ---------- + filepath : path-like or str + Path to the Exodus mesh file + """ - return cls(mesh, values) + #TODO + raise NotImplementedError("Constructor not yet implemented.") - def to_xml_element(self): - """Return XML representation of the mesh - Returns - ------- - element : lxml.etree._Element - XML element containing mesh data +class TemperatureField(ScalarField): + """Temperature field. - """ - element = super().to_xml_element() + Attributes + ---------- + mesh : Mesh + Spatial mesh associated with the field + values : iterable of float + List of values associated with each mesh cell - subelement = ET.SubElement(element, "mesh") - subelement.text = self.mesh.to_xml_element() + """ + def __init__(self, mesh, values): + """Initialization. - return element + Parameters + ---------- + mesh : Mesh + Spatial mesh associated with the field + values : iterable of float + List of values associated with each mesh cell + + """ + super().__init__(mesh, values) diff --git a/openmc/settings.py b/openmc/settings.py index ab121980f54..97fad1c03bb 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -4,6 +4,7 @@ from math import ceil from numbers import Integral, Real from pathlib import Path +import textwrap import lxml.etree as ET import warnings @@ -309,6 +310,10 @@ class Settings: sections be loaded at all temperatures within the range. 'multipole' is a boolean indicating whether or not the windowed multipole method should be used to evaluate resolved resonance cross sections. + temperature_field : openmc.TemperatureField + Temperature field based on a geometric mesh used to specify temperatures + in a model. Temperatures declared from a mesh take precedence over + all other temperature definition (cell, material, and global). trace : tuple or list Show detailed information about a single particle, indicated by three integers: the batch number, generation number, and particle number @@ -1686,7 +1691,9 @@ def _create_temperature_field_subelement(self, root, mesh_memo=None): # Add temperature values subelement = ET.SubElement(element, "values") - subelement.text = " ".join([str(i) for i in self.temperature_field.values]) + subelement.text = "\n ".join( + textwrap.wrap(" ".join( + [str(i) for i in self.temperature_field.values]), 80)) def _create_trigger_subelement(self, root): if self._trigger_active is not None: From d20e3ea7793bb21e1960ae9117cdc7d380a9af47 Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Fri, 19 Dec 2025 16:41:11 -0600 Subject: [PATCH 10/18] Regression test --- .../temperature_field/inputs_true.dat | 48 ++++++++++++++ .../temperature_field/results_true.dat | 2 + .../temperature_field/test.py | 66 +++++++++++++++++++ 3 files changed, 116 insertions(+) create mode 100644 tests/regression_tests/temperature_field/inputs_true.dat create mode 100644 tests/regression_tests/temperature_field/results_true.dat create mode 100755 tests/regression_tests/temperature_field/test.py diff --git a/tests/regression_tests/temperature_field/inputs_true.dat b/tests/regression_tests/temperature_field/inputs_true.dat new file mode 100644 index 00000000000..d4fec68fb12 --- /dev/null +++ b/tests/regression_tests/temperature_field/inputs_true.dat @@ -0,0 +1,48 @@ + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 200 + 20 + + + 0.0 0.0 0.0 5.0 5.0 5.0 + + + true + + + + 1 + 294.0 394.0 494.0 594.0 694.0 794.0 894.0 994.0 + + + 2 2 2 + 0.0 0.0 0.0 + 5.0 5.0 5.0 + + true + 1000 + + diff --git a/tests/regression_tests/temperature_field/results_true.dat b/tests/regression_tests/temperature_field/results_true.dat new file mode 100644 index 00000000000..744dbf08866 --- /dev/null +++ b/tests/regression_tests/temperature_field/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.507597E+00 1.811642E-02 diff --git a/tests/regression_tests/temperature_field/test.py b/tests/regression_tests/temperature_field/test.py new file mode 100755 index 00000000000..66f2aa4afce --- /dev/null +++ b/tests/regression_tests/temperature_field/test.py @@ -0,0 +1,66 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def temperature_field_model(): + """Create a temperature field from a regular mesh over a box with + different temperature for each cell. + + """ + model = openmc.Model() + + # Material + mat = openmc.Material() + mat.add_nuclide("U235", 0.2) + mat.add_nuclide("U238", 0.8) + mat.add_element("O", 2.0) + mat.add_element("H", 4.0) + mat.set_density("g/cm3", 5.0) + mat.add_s_alpha_beta('c_H_in_H2O') + model.materials = openmc.Materials([mat]) + + # Create mesh + dim = 2 + lower_left = (0., 0., 0.) + upper_right = (5.0, 5.0, 5.0) + mesh = openmc.RegularMesh() + mesh.lower_left = lower_left + mesh.upper_right = upper_right + mesh.dimension = (dim, dim, dim) + + # Temperature values + temperature_values = [294.0 + i * 100 for i in range(dim**3)] + + # Create the temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Create geometry + box = openmc.model.RectangularParallelepiped( + xmin=lower_left[0], xmax=upper_right[0], + ymin=lower_left[1], ymax=upper_right[1], + zmin=lower_left[2], zmax=upper_right[2], + boundary_type="reflective" + ) + cell = openmc.Cell(fill=mat, region=-box) + model.geometry = openmc.Geometry([cell]) + + # Settings + settings = openmc.Settings() + settings.batches = 20 + settings.particles = 200 + settings.temperature_field = temperature_field + spatial_dist = openmc.stats.Box(lower_left, upper_right) + settings.source = openmc.IndependentSource( + space=spatial_dist, constraints={"fissionable": True}) + settings.temperature = {'tolerance': 1000, 'multipole': True} + model.settings = settings + + return model + + +def test_temperature_field(temperature_field_model): + harness = PyAPITestHarness('statepoint.20.h5', temperature_field_model) + harness.main() From 6c9b91843ccbfbf73458a1867fd14ed315589e0c Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Fri, 19 Dec 2025 20:37:54 -0600 Subject: [PATCH 11/18] Update documentation --- docs/source/pythonapi/base.rst | 11 ++++++++++ docs/source/usersguide/materials.rst | 31 ++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+) diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index ce2f6f0f857..3770bc628bb 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -168,6 +168,17 @@ Meshes openmc.SphericalMesh openmc.UnstructuredMesh +Fields +------ + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclassinherit.rst + + openmc.ScalarField + openmc.TemperatureField + Geometry Plotting ----------------- diff --git a/docs/source/usersguide/materials.rst b/docs/source/usersguide/materials.rst index 83af5580574..b33b1421f40 100644 --- a/docs/source/usersguide/materials.rst +++ b/docs/source/usersguide/materials.rst @@ -191,6 +191,37 @@ attribute, e.g., :attr:`Material.temperature` or :attr:`Cell.temperature` attributes, respectively. +Alternatively, temperatures can be specified using a temperature field composed +of a geometric mesh and a map that associates a temperature value to each mesh +cell. During the simulation, temperatures associated with particles are also +updated every time a temperature field cell surface is crossed. While a particle +is contained inside the temperature mesh, temperatures from the temperature field +take precedence over temperature declared for a cell, a material or globally. + +The following example shows how to specify temperatures using a temperarature +field based on a regular mesh: + +.. code-block:: python + + # Define a mesh (regular mesh) + dim = 5 + mesh = openmc.RegularMesh() + mesh.lower_left = (0., 0., 0.) + mesh.upper_right = (10.0, 10.0, 10.0) + mesh.dimension = (dim, dim, dim) + + # Define temperature values for each cell + temperature_values = [273.0 + i * 10 for i in range(dim**3)] + + # Create a temperature field + temperature_field = openmc.TemperatureField(mesh, temperature_values) + + # Register the temperature field in the settings + settings = openmc.Settings() + settings.temperature_field = temperature_field + +.. note:: Temperature fields are currently limited to structured meshes only. + ----------------- Material Mixtures ----------------- From 4b660d6224448977001f037e7a11fff6c2029b53 Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Fri, 19 Dec 2025 21:22:50 -0600 Subject: [PATCH 12/18] Create a main C++ ScalarField class --- include/openmc/field.h | 49 ++++++++++++++++++++++++++---------------- src/field.cpp | 14 ++++++------ 2 files changed, 38 insertions(+), 25 deletions(-) diff --git a/include/openmc/field.h b/include/openmc/field.h index 842886b3fd0..ced871146d9 100644 --- a/include/openmc/field.h +++ b/include/openmc/field.h @@ -6,13 +6,12 @@ namespace openmc { -class TemperatureField { - +class ScalarField { public: //---------------------------------------------------------------------------- // Constructors - TemperatureField() {}; - TemperatureField(Mesh* mesh_ptr, vector values); + ScalarField() = default; + ScalarField(Mesh* mesh_ptr, vector values); //---------------------------------------------------------------------------- // Methods @@ -26,6 +25,34 @@ class TemperatureField { //! \return The distance in cm to the next mesh boundary double distance_to_next_boundary(const Position& r, const Direction& d); + //---------------------------------------------------------------------------- + // Accessors + + // Mesh pointer + Mesh* mesh_ptr() const { return mesh_ptr_; } + + // Values + double& value(int i) { return values_[i]; } + const double& value(int i) const { return values_[i]; } + const vector& values() const { return values_; } + +private: + //---------------------------------------------------------------------------- + // Data members + Mesh* mesh_ptr_; //!< Pointer to the geometric mesh + vector values_; //!< Values associated with each mesh cell +}; + +class TemperatureField : public ScalarField { +public: + //---------------------------------------------------------------------------- + // Constructors + TemperatureField() = default; + TemperatureField(Mesh* mesh_ptr, vector values) : ScalarField(mesh_ptr, values) {}; + + //---------------------------------------------------------------------------- + // Methods + //! Returns the temperature in Kelvin corresponding to the given position. // //! \param[in] r Position of the particle @@ -46,20 +73,6 @@ class TemperatureField { // //! \param[inout] p Particle void update_particle_temperature(Particle& p); - - //---------------------------------------------------------------------------- - // Accessors - - // Temperature values - double& value(int i) { return values_[i]; } - const double& value(int i) const { return values_[i]; } - const vector& values() const { return values_; } - -private: - //---------------------------------------------------------------------------- - // Data members - Mesh* mesh_ptr_; //!< Pointer to the geometric mesh - vector values_; //!< Temperature values }; } // namespace openmc diff --git a/src/field.cpp b/src/field.cpp index dab214efcf9..55d9083a7d2 100644 --- a/src/field.cpp +++ b/src/field.cpp @@ -6,7 +6,7 @@ namespace openmc { -TemperatureField::TemperatureField(Mesh* mesh_ptr, vector values) +ScalarField::ScalarField(Mesh* mesh_ptr, vector values) { this->mesh_ptr_ = mesh_ptr; for (double v : values) { @@ -14,26 +14,26 @@ TemperatureField::TemperatureField(Mesh* mesh_ptr, vector values) } } -double TemperatureField::distance_to_next_boundary( +double ScalarField::distance_to_next_boundary( const Position& r, const Direction& u) { if (this->mesh_ptr_ != nullptr) { return this->mesh_ptr_->distance_to_next_boundary(r, u); } else { - fatal_error("No mesh found for the temperature field!"); + fatal_error("No mesh found for the scalar field!"); } } double TemperatureField::get_temperature(const Position& r) { - if (this->mesh_ptr_ != nullptr) { + if (this->mesh_ptr() != nullptr) { // Get bin from position - int i = this->mesh_ptr_->get_bin(r); + int i = this->mesh_ptr()->get_bin(r); // If we have a bin, we use it to locate the value - if (i >= 0 && i < this->values_.size()) { - return this->values_[i]; + if (i >= 0 && i < this->values().size()) { + return this->value(i); } // No values were found (outside the mesh) From 3b3b32f97faeb68aad749983864e46116a1b1fe1 Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Fri, 19 Dec 2025 23:26:32 -0600 Subject: [PATCH 13/18] Check that size of mesh and size of values are consistent in the field - C++ side --- include/openmc/field.h | 22 +++++++++++++++----- src/field.cpp | 47 +++++++++++++++++++++++------------------- 2 files changed, 43 insertions(+), 26 deletions(-) diff --git a/include/openmc/field.h b/include/openmc/field.h index ced871146d9..2bfb71fa63b 100644 --- a/include/openmc/field.h +++ b/include/openmc/field.h @@ -11,7 +11,8 @@ class ScalarField { //---------------------------------------------------------------------------- // Constructors ScalarField() = default; - ScalarField(Mesh* mesh_ptr, vector values); + ScalarField(Mesh* mesh_ptr, vector values, const std::string& field_name); + ScalarField(Mesh* mesh_ptr, vector values) : ScalarField(mesh_ptr, values, "ScalarField") {}; //---------------------------------------------------------------------------- // Methods @@ -28,8 +29,17 @@ class ScalarField { //---------------------------------------------------------------------------- // Accessors + // Field type + const std::string& field_type() const { return this->field_type_; } + // Mesh pointer - Mesh* mesh_ptr() const { return mesh_ptr_; } + Mesh* mesh_ptr() const { + if (this->mesh_ptr_ == nullptr) { + fatal_error(fmt::format("No mesh found for {}!", this->field_type_)); + } else { + return this->mesh_ptr_; + } + } // Values double& value(int i) { return values_[i]; } @@ -39,8 +49,9 @@ class ScalarField { private: //---------------------------------------------------------------------------- // Data members - Mesh* mesh_ptr_; //!< Pointer to the geometric mesh - vector values_; //!< Values associated with each mesh cell + std::string field_type_; //! Name of field type + Mesh* mesh_ptr_; //!< Pointer to the geometric mesh + vector values_; //!< Values associated with each mesh cell }; class TemperatureField : public ScalarField { @@ -48,7 +59,8 @@ class TemperatureField : public ScalarField { //---------------------------------------------------------------------------- // Constructors TemperatureField() = default; - TemperatureField(Mesh* mesh_ptr, vector values) : ScalarField(mesh_ptr, values) {}; + TemperatureField(Mesh* mesh_ptr, vector values) + : ScalarField(mesh_ptr, values, "TemperatureField") {}; //---------------------------------------------------------------------------- // Methods diff --git a/src/field.cpp b/src/field.cpp index 55d9083a7d2..2926f4b7de6 100644 --- a/src/field.cpp +++ b/src/field.cpp @@ -6,9 +6,24 @@ namespace openmc { -ScalarField::ScalarField(Mesh* mesh_ptr, vector values) +ScalarField::ScalarField( + Mesh* mesh_ptr, vector values, const std::string& field_type) { - this->mesh_ptr_ = mesh_ptr; + this->field_type_ = field_type; + + if (mesh_ptr != nullptr) { + this->mesh_ptr_ = mesh_ptr; + } else { + fatal_error(fmt::format("No mesh found for {}!", field_type)); + } + + if (this->mesh_ptr_->n_bins() != values.size()) { + fatal_error( + fmt::format("The number of bins in the mesh is not consistent with the " + "number of values declared in {}!", + field_type)); + } + for (double v : values) { this->values_.push_back(v); } @@ -17,31 +32,21 @@ ScalarField::ScalarField(Mesh* mesh_ptr, vector values) double ScalarField::distance_to_next_boundary( const Position& r, const Direction& u) { - if (this->mesh_ptr_ != nullptr) { - return this->mesh_ptr_->distance_to_next_boundary(r, u); - } else { - fatal_error("No mesh found for the scalar field!"); - } + return this->mesh_ptr()->distance_to_next_boundary(r, u); } double TemperatureField::get_temperature(const Position& r) { - if (this->mesh_ptr() != nullptr) { - - // Get bin from position - int i = this->mesh_ptr()->get_bin(r); - - // If we have a bin, we use it to locate the value - if (i >= 0 && i < this->values().size()) { - return this->value(i); - } - - // No values were found (outside the mesh) - return -1.0; + // Get bin from position + int i = this->mesh_ptr()->get_bin(r); - } else { - fatal_error("No mesh found for the temperature field!"); + // If we have a bin, we use it to locate the value + if (i >= 0 && i < this->values().size()) { + return this->value(i); } + + // No values were found (outside the mesh) + return -1.0; } double TemperatureField::get_sqrtkT(const Position& r) From 6dd684da8b712c56f42997e83e0e9e36fb0f1f8f Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Fri, 19 Dec 2025 23:48:52 -0600 Subject: [PATCH 14/18] Formatting --- include/openmc/field.h | 9 ++++++--- include/openmc/particle_data.h | 1 - include/openmc/settings.h | 6 +++--- src/finalize.cpp | 2 +- src/geometry_aux.cpp | 6 +++--- src/particle.cpp | 4 ++-- src/simulation.cpp | 34 +++++++++++++++++----------------- 7 files changed, 32 insertions(+), 30 deletions(-) diff --git a/include/openmc/field.h b/include/openmc/field.h index 2bfb71fa63b..58b2441a52e 100644 --- a/include/openmc/field.h +++ b/include/openmc/field.h @@ -11,8 +11,10 @@ class ScalarField { //---------------------------------------------------------------------------- // Constructors ScalarField() = default; - ScalarField(Mesh* mesh_ptr, vector values, const std::string& field_name); - ScalarField(Mesh* mesh_ptr, vector values) : ScalarField(mesh_ptr, values, "ScalarField") {}; + ScalarField( + Mesh* mesh_ptr, vector values, const std::string& field_name); + ScalarField(Mesh* mesh_ptr, vector values) + : ScalarField(mesh_ptr, values, "ScalarField") {}; //---------------------------------------------------------------------------- // Methods @@ -33,7 +35,8 @@ class ScalarField { const std::string& field_type() const { return this->field_type_; } // Mesh pointer - Mesh* mesh_ptr() const { + Mesh* mesh_ptr() const + { if (this->mesh_ptr_ == nullptr) { fatal_error(fmt::format("No mesh found for {}!", this->field_type_)); } else { diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 829d503d0fa..6fb1be09830 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -774,7 +774,6 @@ class ParticleData : public GeometryState { // Next event identifier for history-based transport int next_event() const { return next_event_; } int& next_event() { return next_event_; } - }; } // namespace openmc diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 6f9979d9e42..e16ee305680 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -69,9 +69,9 @@ extern bool delayed_photon_scaling; //!< Scale fission photon yield to include delayed extern "C" bool entropy_on; //!< calculate Shannon entropy? extern "C" bool - event_based; //!< use event-based mode (instead of history-based) -extern bool temperature_field_on ; //!< Is there a temperature field defined? -extern bool ifp_on; //!< Use IFP for kinetics parameters? + event_based; //!< use event-based mode (instead of history-based) +extern bool temperature_field_on; //!< Is there a temperature field defined? +extern bool ifp_on; //!< Use IFP for kinetics parameters? extern bool legendre_to_tabular; //!< convert Legendre distributions to tabular? extern bool material_cell_offsets; //!< create material cells offsets? extern "C" bool output_summary; //!< write summary.h5? diff --git a/src/finalize.cpp b/src/finalize.cpp index b9c000528d7..b2369745e18 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -159,7 +159,7 @@ int openmc_finalize() simulation::entropy_mesh = nullptr; simulation::ufs_mesh = nullptr; - simulation::temperature_field; + simulation::temperature_field = TemperatureField(); data::energy_max = {INFTY, INFTY, INFTY, INFTY}; data::energy_min = {0.0, 0.0, 0.0, 0.0}; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index bdb340b963f..1e81851bdbd 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -269,15 +269,15 @@ void get_temperatures( // could limit the declarations of temperature from the temperature field to // impacted nuclides only. if (settings::temperature_field_on) { - for (auto t: simulation::temperature_field.values()) { + for (auto t : simulation::temperature_field.values()) { // Nuclide temperatures - for (size_t i = 0; i < nuc_temps.size() ; i++){ + for (size_t i = 0; i < nuc_temps.size(); i++) { if (!contains(nuc_temps[i], t)) { nuc_temps[i].push_back(t); } } // Thermal scattering temperatures - for (size_t i = 0; i < thermal_temps.size() ; i++){ + for (size_t i = 0; i < thermal_temps.size(); i++) { if (!contains(thermal_temps[i], t)) { thermal_temps[i].push_back(t); } diff --git a/src/particle.cpp b/src/particle.cpp index 0f68dba912e..ad469166d17 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -262,8 +262,8 @@ void Particle::event_advance() (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY; // Select smaller distance - double distance = - std::min({boundary().distance(), collision_distance(), distance_cutoff, distance_tmesh}); + double distance = std::min({boundary().distance(), collision_distance(), + distance_cutoff, distance_tmesh}); // Prepare the stop condition if (distance == distance_cutoff) { diff --git a/src/simulation.cpp b/src/simulation.cpp index 2255522ba2d..f0e5a547c50 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -814,23 +814,23 @@ void transport_history_based_single_particle(Particle& p) } if (p.alive()) { switch (p.next_event()) { - case EVENT_CROSS_SURFACE: - p.event_cross_surface(); - break; - case EVENT_COLLIDE: - p.event_collide(); - break; - case EVENT_CROSS_TEMPERATURE_MESH: - p.event_cross_temperature_mesh(); - break; - case EVENT_TIME_CUTOFF: - p.wgt() = 0.0; - break; - default: - fatal_error(fmt::format( - "Unknown event '{}' in history-based transport!", p.next_event())); - break; - } + case EVENT_CROSS_SURFACE: + p.event_cross_surface(); + break; + case EVENT_COLLIDE: + p.event_collide(); + break; + case EVENT_CROSS_TEMPERATURE_MESH: + p.event_cross_temperature_mesh(); + break; + case EVENT_TIME_CUTOFF: + p.wgt() = 0.0; + break; + default: + fatal_error(fmt::format( + "Unknown event '{}' in history-based transport!", p.next_event())); + break; + } } p.event_revive_from_secondary(); } From 76cdf20b838acb54b9069dee88cfa3e553e9c63a Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Sat, 20 Dec 2025 00:25:27 -0600 Subject: [PATCH 15/18] Manage incompatibility with event-based mode --- src/settings.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/settings.cpp b/src/settings.cpp index 1d4fb8d3041..60104d86a31 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -1228,6 +1228,10 @@ void read_settings_xml(pugi::xml_node root) // Check whether to use event-based parallelism if (check_for_node(root, "event_based")) { event_based = get_node_value_bool(root, "event_based"); + if (temperature_field_on && event_based) { + fatal_error( + "Event-based transport is not yet compatible with temperature fields."); + } } // Check whether material cell offsets should be generated From 49bfb76ea0f9ddeb1edea93fa87905122ac301d4 Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Sat, 20 Dec 2025 01:11:09 -0600 Subject: [PATCH 16/18] Unit test for the distance to next boundary method --- include/openmc/mesh.h | 6 ++--- src/mesh.cpp | 6 ++--- tests/cpp_unit_tests/test_mesh.cpp | 36 ++++++++++++++++++++++++++++++ 3 files changed, 42 insertions(+), 6 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 43988d0049a..1dccf91ea90 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -181,7 +181,7 @@ class Mesh { //! \param[in] r Position of the particle //! \param[in] u Particle direction //! \return Distance to the next boundary - virtual double distance_to_next_boundary(Position r, Position u) const = 0; + virtual double distance_to_next_boundary(Position r, Direction u) const = 0; //! Get bin at a given position in space // @@ -312,7 +312,7 @@ class StructuredMesh : public Mesh { void surface_bins_crossed(Position r0, Position r1, const Direction& u, vector& bins) const override; - double distance_to_next_boundary(Position r, Position u) const override; + double distance_to_next_boundary(Position r, Direction u) const override; //! Determine which cell or surface bins were crossed by a particle // @@ -696,7 +696,7 @@ class UnstructuredMesh : public Mesh { UnstructuredMesh(pugi::xml_node node); UnstructuredMesh(hid_t group); - double distance_to_next_boundary(Position r, Position u) const override; + double distance_to_next_boundary(Position r, Direction u) const override; static const std::string mesh_type; virtual std::string get_mesh_type() const override; diff --git a/src/mesh.cpp b/src/mesh.cpp index 0c696eed220..afa027d21d6 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -720,7 +720,7 @@ UnstructuredMesh::UnstructuredMesh(hid_t group) : Mesh(group) } } -double UnstructuredMesh::distance_to_next_boundary(Position r, Position u) const +double UnstructuredMesh::distance_to_next_boundary(Position r, Direction u) const { fatal_error("Not implemented"); return -1.0; @@ -1181,12 +1181,12 @@ void StructuredMesh::surface_bins_crossed( raytrace_mesh(r0, r1, u, SurfaceAggregator(this, bins)); } -double StructuredMesh::distance_to_next_boundary(Position r, Position u) const +double StructuredMesh::distance_to_next_boundary(Position r, Direction u) const { Position global_r = r; Position local_r = local_coords(r); - double distance = INFTY; + double distance = 0.0; const int n = n_dimension_; bool in_mesh; diff --git a/tests/cpp_unit_tests/test_mesh.cpp b/tests/cpp_unit_tests/test_mesh.cpp index 24c4f77373f..1881e6d729e 100644 --- a/tests/cpp_unit_tests/test_mesh.cpp +++ b/tests/cpp_unit_tests/test_mesh.cpp @@ -255,3 +255,39 @@ TEST_CASE("Test multiple meshes HDF5 roundtrip - spherical") REQUIRE(regular_mesh_hdf5->lower_left() == regular_mesh_xml->lower_left()); REQUIRE(regular_mesh_hdf5->upper_right() == regular_mesh_xml->upper_right()); } + +TEST_CASE("Test distance_to_next_boundary() method for regular mesh") +{ + // The XML data as a string + std::string xml_string = R"( + + 2 2 2 + -1 -1 -1 + 1 1 1 + + )"; + + // Create the mesh from a file + pugi::xml_document doc; + pugi::xml_parse_result result = doc.load_string(xml_string.c_str()); + pugi::xml_node root = doc.child("mesh"); + auto mesh = RegularMesh(root); + + Position r; + Position u; + + // Test inside the mesh + r = Position(0.0, 0.0, 0.0); + u = Position(1.0, 0.0, 0.0); + REQUIRE(mesh.distance_to_next_boundary(r, u) == 1.0); + + // Test outside the mesh, going toward the mesh + r = Position(-1.5, 0.0, 0.0); + u = Position(1.0, 0.0, 0.0); + REQUIRE(mesh.distance_to_next_boundary(r, u) == 0.5); + + // Test outside the mesh, not going toward the mesh + r = Position(-2.0, 0.0, 0.0); + u = Position(-1.0, 0.0, 0.0); + REQUIRE(mesh.distance_to_next_boundary(r, u) == INFTY); +} \ No newline at end of file From 5e79b7b6f0848939e02fc225c0868b1ee6885fb2 Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Thu, 15 Jan 2026 12:48:40 -0600 Subject: [PATCH 17/18] Missing init file in new regression test --- tests/regression_tests/temperature_field/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 tests/regression_tests/temperature_field/__init__.py diff --git a/tests/regression_tests/temperature_field/__init__.py b/tests/regression_tests/temperature_field/__init__.py new file mode 100644 index 00000000000..e69de29bb2d From 0c670167203e5795d0d1234be22c307a9f273a79 Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Thu, 15 Jan 2026 13:24:44 -0600 Subject: [PATCH 18/18] Formatting --- src/geometry_aux.cpp | 2 +- src/mesh.cpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 1e81851bdbd..b7d954eaad9 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -17,8 +17,8 @@ #include "openmc/lattice.h" #include "openmc/material.h" #include "openmc/settings.h" -#include "openmc/surface.h" #include "openmc/simulation.h" +#include "openmc/surface.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/filter_cell_instance.h" #include "openmc/tallies/filter_distribcell.h" diff --git a/src/mesh.cpp b/src/mesh.cpp index afa027d21d6..809128da93f 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -720,7 +720,8 @@ UnstructuredMesh::UnstructuredMesh(hid_t group) : Mesh(group) } } -double UnstructuredMesh::distance_to_next_boundary(Position r, Direction u) const +double UnstructuredMesh::distance_to_next_boundary( + Position r, Direction u) const { fatal_error("Not implemented"); return -1.0;