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/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 ----------------- 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/include/openmc/field.h b/include/openmc/field.h new file mode 100644 index 00000000000..58b2441a52e --- /dev/null +++ b/include/openmc/field.h @@ -0,0 +1,94 @@ +#ifndef OPENMC_FIELD_H +#define OPENMC_FIELD_H + +#include "openmc/mesh.h" +#include "openmc/vector.h" + +namespace openmc { + +class ScalarField { +public: + //---------------------------------------------------------------------------- + // 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") {}; + + //---------------------------------------------------------------------------- + // 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); + + //---------------------------------------------------------------------------- + // Accessors + + // Field type + const std::string& field_type() const { return this->field_type_; } + + // Mesh pointer + 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]; } + const double& value(int i) const { return values_[i]; } + const vector& values() const { return values_; } + +private: + //---------------------------------------------------------------------------- + // Data members + 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 { +public: + //---------------------------------------------------------------------------- + // Constructors + TemperatureField() = default; + TemperatureField(Mesh* mesh_ptr, vector values) + : ScalarField(mesh_ptr, values, "TemperatureField") {}; + + //---------------------------------------------------------------------------- + // Methods + + //! 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); + + //! 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); +}; + +} // namespace openmc +#endif // OPENMC_FIELD_H diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 5c9272e93b9..1dccf91ea90 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -173,6 +173,16 @@ 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. + // + //! \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, Direction u) const = 0; + //! Get bin at a given position in space // //! \param[in] r Position to get bin for @@ -302,6 +312,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, Direction u) const override; + //! Determine which cell or surface bins were crossed by a particle // //! \param[in] r0 Previous position of the particle @@ -684,6 +696,8 @@ class UnstructuredMesh : public Mesh { UnstructuredMesh(pugi::xml_node node); UnstructuredMesh(hid_t group); + 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/include/openmc/particle.h b/include/openmc/particle.h index 0f37719b94f..8daf8290b12 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -68,6 +68,7 @@ class Particle : public ParticleData { void event_calculate_xs(); void event_advance(); void event_cross_surface(); + void event_cross_temperature_mesh(); void event_collide(); void event_revive_from_secondary(); void event_death(); diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index fdacfa765b3..6fb1be09830 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,10 @@ 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/include/openmc/settings.h b/include/openmc/settings.h index b369c99fef8..e16ee305680 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -69,8 +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 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/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..463d997e602 --- /dev/null +++ b/openmc/field.py @@ -0,0 +1,67 @@ +from abc import ABC, abstractmethod + +import openmc + + +class ScalarField(ABC): + """Scalar field defined on a geometric mesh. + + Attributes + ---------- + mesh : Mesh + Spatial mesh associated with the field + values : iterable of float + List of values associated with each mesh cell + + """ + 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_exodus_file(cls, filepath): + """Construct a ScalarField from an Exodus mesh file + + Parameters + ---------- + filepath : path-like or str + Path to the Exodus mesh file + + """ + #TODO + raise NotImplementedError("Constructor not yet implemented.") + + +class TemperatureField(ScalarField): + """Temperature field. + + Attributes + ---------- + mesh : Mesh + Spatial mesh associated with the field + values : iterable of float + List of values associated with each mesh cell + + """ + 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 + + """ + super().__init__(mesh, values) diff --git a/openmc/settings.py b/openmc/settings.py index 43c1fe06982..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 @@ -14,6 +15,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 @@ -308,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 @@ -400,6 +406,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 +721,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 +1668,33 @@ 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 = "\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: trigger_element = ET.SubElement(root, "trigger") @@ -2135,6 +2180,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 +2467,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 +2582,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..2926f4b7de6 --- /dev/null +++ b/src/field.cpp @@ -0,0 +1,80 @@ +#include "openmc/field.h" +#include "openmc/cell.h" +#include "openmc/constants.h" +#include "openmc/mesh.h" +#include "openmc/vector.h" + +namespace openmc { + +ScalarField::ScalarField( + Mesh* mesh_ptr, vector values, const std::string& field_type) +{ + 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); + } +} + +double ScalarField::distance_to_next_boundary( + const Position& r, const Direction& u) +{ + return this->mesh_ptr()->distance_to_next_boundary(r, u); +} + +double TemperatureField::get_temperature(const Position& r) +{ + // 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; +} + +double TemperatureField::get_sqrtkT(const Position& r) +{ + double temperature = this->get_temperature(r); + if (temperature >= 0) { + return sqrt(temperature * K_BOLTZMANN); + } + return -1.0; +} + +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/finalize.cpp b/src/finalize.cpp index 9ee94340997..b2369745e18 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 = TemperatureField(); + 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..b7d954eaad9 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -17,6 +17,7 @@ #include "openmc/lattice.h" #include "openmc/material.h" #include "openmc/settings.h" +#include "openmc/simulation.h" #include "openmc/surface.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/filter_cell_instance.h" @@ -261,6 +262,28 @@ void get_temperatures( } } } + + // 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); + } + } + } + } } //============================================================================== diff --git a/src/mesh.cpp b/src/mesh.cpp index 610d057cf08..809128da93f 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -720,6 +720,13 @@ UnstructuredMesh::UnstructuredMesh(hid_t group) : Mesh(group) } } +double UnstructuredMesh::distance_to_next_boundary( + Position r, Direction u) const +{ + fatal_error("Not implemented"); + return -1.0; +} + void UnstructuredMesh::determine_bounds() { double xmin = INFTY; @@ -1175,6 +1182,49 @@ void StructuredMesh::surface_bins_crossed( raytrace_mesh(r0, r1, u, SurfaceAggregator(this, bins)); } +double StructuredMesh::distance_to_next_boundary(Position r, Direction u) const +{ + Position global_r = r; + Position local_r = local_coords(r); + + double distance = 0.0; + const int n = n_dimension_; + 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..ad469166d17 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -193,6 +193,11 @@ 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) { + simulation::temperature_field.update_particle_temperature(*this); + } } // Write particle track. @@ -231,7 +236,7 @@ void Particle::event_calculate_xs() void Particle::event_advance() { - // Find the distance to the nearest boundary + // Find the distance to the nearest geometry boundary boundary() = distance_to_boundary(*this); // Sample a distance to collision @@ -243,14 +248,33 @@ 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_boundary(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 - double distance = - std::min({boundary().distance(), collision_distance(), distance_cutoff}); + // Select smaller distance + double distance = std::min({boundary().distance(), collision_distance(), + distance_cutoff, distance_tmesh}); + + // Prepare the stop condition + if (distance == distance_cutoff) { + next_event() = EVENT_TIME_CUTOFF; + } else if (distance == boundary().distance()) { + next_event() = EVENT_CROSS_SURFACE; + } else if (distance == distance_tmesh) { + next_event() = EVENT_CROSS_TEMPERATURE_MESH; + } else if (distance == collision_distance()) { + next_event() = EVENT_COLLIDE; + } // Advance particle in space and time this->move_distance(distance); @@ -278,11 +302,12 @@ void Particle::event_advance() if (!model::active_tallies.empty()) { score_track_derivative(*this, distance); } +} - // Set particle weight to zero if it hit the time boundary - if (distance == distance_cutoff) { - wgt() = 0.0; - } +void Particle::event_cross_temperature_mesh() +{ + // Update temperature of the particle + simulation::temperature_field.update_particle_temperature(*this); } void Particle::event_cross_surface() @@ -322,6 +347,12 @@ void Particle::event_cross_surface() } event() = TallyEvent::SURFACE; } + + // Update temperature of the particle if temperature field + if (settings::temperature_field_on) { + simulation::temperature_field.update_particle_temperature(*this); + } + // 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..60104d86a31 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 + 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")); @@ -1189,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 diff --git a/src/simulation.cpp b/src/simulation.cpp index b536ae5881f..f0e5a547c50 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -3,10 +3,12 @@ #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" #include "openmc/event.h" +#include "openmc/field.h" #include "openmc/geometry_aux.h" #include "openmc/ifp.h" #include "openmc/material.h" @@ -320,6 +322,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; @@ -809,10 +813,23 @@ void transport_history_based_single_particle(Particle& p) p.event_advance(); } if (p.alive()) { - if (p.collision_distance() > p.boundary().distance()) { + switch (p.next_event()) { + case EVENT_CROSS_SURFACE: p.event_cross_surface(); - } else if (p.alive()) { + 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(); 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 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 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()