DNP drift (regular mesh only)

CMakeLists.txt

@@ -387,10 +387,12 @@ list(APPEND libopenmc_SOURCES
   src/distribution_energy.cpp
   src/distribution_multi.cpp
   src/distribution_spatial.cpp
+  src/dnp_drift.cpp
   src/eigenvalue.cpp
   src/endf.cpp
   src/error.cpp
   src/event.cpp
+  src/field.cpp
   src/file_utils.cpp
   src/finalize.cpp
   src/geometry.cpp
@@ -442,6 +444,7 @@ list(APPEND libopenmc_SOURCES
   src/simulation.cpp
   src/source.cpp
   src/state_point.cpp
+  src/streamline_integrator.cpp
   src/string_utils.cpp
   src/summary.cpp
   src/surface.cpp

docs/source/capi/index.rst

@@ -1007,3 +1007,9 @@ Functions
    :type scores: const int*
    :return: Return status (negative if an error occurred)
    :rtype: int
+
+.. doxygenfunction:: openmc_temperature_field_set_temperature
+
+.. doxygenfunction:: openmc_temperature_field_size
+
+.. doxygenfunction:: openmc_temperature_field_get_value

docs/source/io_formats/settings.rst

@@ -1469,6 +1469,27 @@ a material default temperature.
 
   *Default*: 293.6 K
 
+.. _temperature_field:
+
+-------------------------------
+``<temperature_field>`` Element
+-------------------------------
+
+The ``<temperature_field>`` element specifies temperatures in a model in Kelvin
+using a temperature field based on a geometric mesh. Temperatures declared using
+a temperature field take precedence over all other temperature definition (cell,
+material, and global). This element has the following sub-elements:
+
+  :mesh:
+    The ID of the geometric mesh on which the temperature field is defined.
+
+  :values:
+    A whitespace-separated list of temperatures in Kelvin, with one value for
+    each mesh element. The number of values must match the number of mesh
+    elements.
+
+  *Default*: None
+
 .. _temperature_method:
 
 --------------------------------

docs/source/pythonapi/base.rst

@@ -172,6 +172,17 @@ Meshes
    openmc.SphericalMesh
    openmc.UnstructuredMesh
 
+Fields
+------
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+   :template: myclassinherit.rst
+
+   openmc.ScalarField
+   openmc.TemperatureField
+
 Geometry Plotting
 -----------------
 

docs/source/usersguide/materials.rst

@@ -191,6 +191,38 @@ attribute, e.g.,
              :attr:`Material.temperature` or :attr:`Cell.temperature`
              attributes, respectively.
 
+Alternatively, temperatures can be specified using a temperature field composed
+of a mesh and a map that associates a temperature value to each mesh element.
+During the simulation, temperatures associated with particles are also updated
+every time a temperature field mesh surface is crossed. While a particle is
+contained inside the temperature mesh, temperatures from the temperature field
+take precedence over any temperature declared for a cell, a material, or
+globally.
+
+The following example shows how to specify temperatures using a temperature
+field based on a regular mesh:
+
+.. code-block:: python
+
+   # Define a 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 mesh element
+   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
 -----------------

include/openmc/capi.h

@@ -256,6 +256,26 @@ int openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides);
 int openmc_tally_set_scores(int32_t index, int n, const char** scores);
 int openmc_tally_set_type(int32_t index, const char* type);
 int openmc_tally_set_writable(int32_t index, bool writable);
+
+//! Set a value at a given index in the temperature field
+//
+//! \param[in] index Index in the field
+//! \param[in] temperature Temperature value to apply
+//! \return Status (negative if an error occurred)
+int openmc_temperature_field_set_temperature(int32_t index, double temperature);
+
+//! Return the temperature field size (number of cells)
+//
+//! \return Number of cell in the temperature field
+size_t openmc_temperature_field_size();
+
+//! Get the value at a given index in the temperature field
+//
+//! \param[in] index Index in the field
+//! \param[out] temperature Temperature value at the given index
+//! \return Status (negative if an error occurred)
+int openmc_temperature_field_get_value(int32_t index, double* temperature);
+
 int openmc_get_weight_windows_index(int32_t id, int32_t* idx);
 int openmc_weight_windows_get_id(int32_t index, int32_t* id);
 int openmc_weight_windows_set_id(int32_t index, int32_t id);

include/openmc/constants.h

@@ -392,6 +392,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_UNDEFINED = 0;
+const int EVENT_CROSS_SURFACE = 1;
+const int EVENT_COLLIDE = 2;
+const int EVENT_TIME_CUTOFF = 3;
+
 } // namespace openmc
 
 #endif // OPENMC_CONSTANTS_H

include/openmc/dnp_drift.h

@@ -0,0 +1,76 @@
+#ifndef OPENMC_DNP_DRIFT_H
+#define OPENMC_DNP_DRIFT_H
+
+#include "openmc/particle_data.h"
+#include "openmc/position.h"
+
+namespace openmc {
+
+const int DNP_DRIFT_TRANSPORT_MAX_ITER =
+  1000; //! Maximum number of iterations allowed in the DNP transport loop.
+const double DNP_DRIFT_DISTANCE_MIN =
+  1.0E-13; //! Minimum distance between two integration steps before stopping
+           //! the DNP transport loop.
+
+enum class Actions { REENTER, ESCAPE, DECAY_IN_PLACE };
+
+//! Adjust position linearly.
+//!
+//! This function is used to adjust the position of a DNP that reached its decay
+//! time between two integration steps (i.e., t is greater than decay_time). The
+//! final position of the DNP is adjusted linearly so that the time
+//! corresponding to this new position is equal to the sampled decay time.
+//!
+//! \param[inout] y_n Current position
+//! \param[in] y_n_minus_1 Previous position
+//! \param[in] t Time at current position
+//! \param[in] dt Time step
+//! \param[in] decay_time Decay time
+void _adjust_position(Position& y_n, const Position& y_n_minus_1, double t,
+  double dt, double decay_time);
+
+//! Adjust time linearly.
+//!
+//! This function is used to adjust the time associated with the position of a
+//! DNP that needs to be stopped between two integration steps. The position C,
+//! where the DNP is stopped, must be located between position A and B. Time is
+//! adjusted linearly.
+//!
+//! \param[inout] t In: time at position B, out: time at position C
+//! \param[in] ta Time at position A
+//! \param[in] pa Position A
+//! \param[in] pb Position B
+//! \param[in] pc Position C (between A and B)
+void _adjust_time(double& t, double ta, const Position& pa, const Position& pb,
+  const Position& pc);
+
+//! Explicitly transport a delayed neutron precursor using streamline
+//! integration.
+//!
+//! \param[inout] site Fission site corresponding to the DNP
+//! \param[in] decay_time Sampled decay time (in seconds)
+//! \param[in] seed Random number generator seed
+//! \return True if the site is considered inside the explicitly modeled part of
+//!         the system. False otherwise.
+bool transport_dnp(SourceSite& site, double decay_time, uint64_t* seed);
+
+//! Reconcile a delayed neutron precursor source site after mesh-based drift
+//! with the model geometry.
+//!
+//! After transporting a delayed neutron precursor on a mesh, the emission site
+//! may land on or near a geometry boundary. This function verifies that the
+//! site is valid within the geometry and applies boundary condition corrections
+//! if needed:
+//!   - Transmission: site is accepted as is.
+//!   - Vacuum: site is accepted if the direction points inward the current
+//!     cell; rejected otherwise.
+//!   - Reflective: direction is reflected if pointing outward.
+//!
+//! \param[inout] site  Fission site corresponding to the DNP
+//! \return True if the site is valid and usable for the next generation,
+//!         false if the site should be discarded.
+bool reconcile_precursor_drift(SourceSite& site);
+
+} // namespace openmc
+
+#endif // OPENMC_DNP_DRIFT_H

include/openmc/event.h

@@ -109,6 +109,10 @@ void process_surface_crossing_events();
 //! Execute the collision event for all particles in this event's buffer
 void process_collision_events();
 
+//! Execute the temperature mesh crossing event for all particles in this
+//! event's buffer
+void process_temperature_mesh_crossing_events();
+
 //! Execute the death event for all particles
 //
 //! \param n_particles The number of particles in the particle buffer