DNP drift (regular mesh only)

CMakeLists.txt

@@ -370,10 +370,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
@@ -425,6 +427,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

@@ -1030,3 +1030,12 @@ Functions
    :type scores: const int*
    :return: Return status (negative if an error occurred)
    :rtype: int
+
+.. c:function:: int openmc_temperature_field_set_temperature(int32_t index, double temperature)
+
+   Set the temperature value of a given cell in the temperature field
+
+   :param int32_t index: Index in the tempererature mesh
+   :param double temperature: Temperature in Kelvin
+   :return: Return status (negative if an error occurred)
+   :rtype: int

docs/source/pythonapi/base.rst

@@ -170,6 +170,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,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
 -----------------

include/openmc/constants.h

@@ -378,6 +378,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,57 @@
+#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 { PLACE_AT_INLET, BLOCK_AT_OUTLET, BLOCK_AT_LOCATION };
+
+//! 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
+bool transport_dnp(SourceSite& site, double decay_time, uint64_t* seed);
+
+} // namespace openmc
+
+#endif // OPENMC_DNP_DRIFT_H

include/openmc/event.h

@@ -107,6 +107,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