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/boundary_condition.h

@@ -21,12 +21,38 @@ class BoundaryCondition {
 public:
   virtual ~BoundaryCondition() = default;
 
+  //! Apply the geometric transformation.
+  //! \param[in] p Particle
+  //! \param[in] surf Current surface
+  //! \param[out] new_r Transformed position
+  //! \param[out] new_u Transformed direction
+  //! \param[out] new_surface New surface index
+  virtual void transform(Particle& p, const Surface& surf, Position& new_r,
+    Direction& new_u, int& new_surface) const = 0;
+
+  //! Handle the crossing of the boundary condition by the particle
+  //! \param[inout] p Particle
+  //! \param[in] surf Current surface
+  //! \param[in] new_r Transformed position
+  //! \param[in] new_u Transformed direction
+  //! \param[in] new_surface New surface index
+  virtual void cross_bc(Particle& p, const Surface& surf, const Position& new_r,
+    const Direction& new_u, int new_surface) const = 0;
+
   //! Perform tracking operations for a particle that strikes the boundary.
   //! \param p The particle that struck the boundary.  This class is not meant
   //!   to directly modify anything about the particle, but it will do so
   //!   indirectly by calling the particle's appropriate cross_*_bc function.
   //! \param surf The specific surface on the boundary the particle struck.
-  virtual void handle_particle(Particle& p, const Surface& surf) const = 0;
+  virtual void handle_particle(Particle& p, const Surface& surf) const
+  {
+    Position new_r;
+    Direction new_u;
+    int new_surface;
+    transform(p, surf, new_r, new_u, new_surface);
+    BoundaryCondition::handle_albedo(p, surf);
+    cross_bc(p, surf, new_r, new_u, new_surface);
+  }
 
   //! Modify the incident particle's weight according to the boundary's albedo.
   //! \param p The particle that struck the boundary.  This function calculates
@@ -64,6 +90,9 @@ class BoundaryCondition {
   //! Return if this BC has an albedo.
   bool has_albedo() const { return (albedo_ > 0.0); }
 
+  //! Return albedo value.
+  double albedo() const { return albedo_; }
+
 private:
   double albedo_ = -1.0;
 };
@@ -74,6 +103,16 @@ class BoundaryCondition {
 
 class VacuumBC : public BoundaryCondition {
 public:
+  // Not used - handle_particle is fully overriden
+  void transform(Particle& p, const Surface& surf, Position& new_r,
+    Direction& new_u, int& new_surface) const override
+  {}
+
+  // Not used - handle_particle is fully overriden
+  void cross_bc(Particle& p, const Surface& surf, const Position& new_r,
+    const Direction& new_u, int new_surface) const override
+  {}
+
   void handle_particle(Particle& p, const Surface& surf) const override;
 
   std::string type() const override { return "vacuum"; }
@@ -85,7 +124,11 @@ class VacuumBC : public BoundaryCondition {
 
 class ReflectiveBC : public BoundaryCondition {
 public:
-  void handle_particle(Particle& p, const Surface& surf) const override;
+  void transform(Particle& p, const Surface& surf, Position& new_r,
+    Direction& new_u, int& new_surface) const override;
+
+  void cross_bc(Particle& p, const Surface& surf, const Position& new_r,
+    const Direction& new_u, int new_surface) const override;
 
   std::string type() const override { return "reflective"; }
 };
@@ -96,7 +139,11 @@ class ReflectiveBC : public BoundaryCondition {
 
 class WhiteBC : public BoundaryCondition {
 public:
-  void handle_particle(Particle& p, const Surface& surf) const override;
+  void transform(Particle& p, const Surface& surf, Position& new_r,
+    Direction& new_u, int& new_surface) const override;
+
+  void cross_bc(Particle& p, const Surface& surf, const Position& new_r,
+    const Direction& new_u, int new_surface) const override;
 
   std::string type() const override { return "white"; }
 };
@@ -128,7 +175,11 @@ class TranslationalPeriodicBC : public PeriodicBC {
 public:
   TranslationalPeriodicBC(int i_surf, int j_surf);
 
-  void handle_particle(Particle& p, const Surface& surf) const override;
+  void transform(Particle& p, const Surface& surf, Position& new_r,
+    Direction& new_u, int& new_surface) const override;
+
+  void cross_bc(Particle& p, const Surface& surf, const Position& new_r,
+    const Direction& new_u, int new_surface) const override;
 
 protected:
   //! Vector along which incident particles will be moved
@@ -145,9 +196,12 @@ class RotationalPeriodicBC : public PeriodicBC {
 public:
   enum PeriodicAxis { x, y, z };
   RotationalPeriodicBC(int i_surf, int j_surf, PeriodicAxis axis);
-  double compute_periodic_rotation(
-    double rise_1, double run_1, double rise_2, double run_2) const;
-  void handle_particle(Particle& p, const Surface& surf) const override;
+
+  void transform(Particle& p, const Surface& surf, Position& new_r,
+    Direction& new_u, int& new_surface) const override;
+
+  void cross_bc(Particle& p, const Surface& surf, const Position& new_r,
+    const Direction& new_u, int new_surface) const override;
 
 protected:
   //! Angle about the axis by which particle coordinates will be rotated

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,75 @@
+#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[in,out] 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[in,out] 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[in,out] 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 with the model geometry
+//! after mesh-based drift.
+//!
+//! After transporting a delayed neutron precursor on a mesh, the emission site
+//! may land outside the model or on a geometry boundary. This function
+//! identifies if the site is valid and applies boundary condition corrections
+//! as needed. If the function returns true, the site is guaranteed to be usable
+//! for the next fission generation.
+//!
+//! \param[in,out] site  Fission site corresponding to the DNP
+//! \param[in] particle_id Particle ID used to initialize the RNG seed for the
+//!                        boundary condition sampling (e.g., white BC).
+//!\return True if the site is usable for the next generation, false if the site
+//!        should be discarded.
+bool reconcile_precursor_drift(SourceSite& site, int64_t particle_id);
+
+} // 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