Specify temperature from a field (structured mesh only)

CMakeLists.txt

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

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 temperature 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,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/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/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

include/openmc/field.h

@@ -0,0 +1,96 @@
+#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<double> values, const std::string& field_name);
+  ScalarField(Mesh* mesh_ptr, vector<double> values)
+    : ScalarField(mesh_ptr, values, "ScalarField") {};
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  //! Returns the distance to the next mesh boundary given 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] current_bin Current bin number
+  //! \param[in] r Position of the particle
+  //! \param[in] u Direction of the particle
+  //! \param[out] bin_next Next bin number
+  //! \return The distance in cm to the next mesh boundary
+  double distance_to_next_boundary(
+    int current_bin, const Position& r, const Direction& d, int& bin_next);
+
+  //----------------------------------------------------------------------------
+  // 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<double>& values() const { return values_; }
+
+private:
+  //----------------------------------------------------------------------------
+  // Data members
+  std::string field_type_; //! Name of field type
+  Mesh* mesh_ptr_;         //!< Pointer to the geometric mesh
+  vector<double> values_;  //!< Values associated with each mesh cell
+};
+
+class TemperatureField : public ScalarField {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors
+  TemperatureField() = default;
+  TemperatureField(Mesh* mesh_ptr, vector<double> values)
+    : ScalarField(mesh_ptr, values, "TemperatureField") {};
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  //! Returns the temperature in Kelvin corresponding to a given bin number
+  //! relative to the mesh.
+  //
+  //! \param[in] r Position of the particle
+  //! \return Temperature in Kelvin
+  double get_temperature(int bin);
+
+  //! Returns the square root of the temperature multiplied by the Boltzmann
+  //! constant in eV for a given bin number relative to the mesh.
+  //
+  //! \param[in] bin Bin number
+  //! \return Sqrt(k_Boltzmann * temperature) in eV
+  double get_sqrtkT(int bin);
+
+  //! Returns the bin number corresponding to the location of the particle.
+  //
+  //! \param[in] r Position of the particle
+  //! \return Corresponding bin number or -1 if outside the mesh
+  int get_bin(const Position& r);
+};
+
+} //  namespace openmc
+#endif // OPENMC_FIELD_H

include/openmc/mesh.h

@@ -193,6 +193,19 @@ class Mesh {
   virtual void surface_bins_crossed(
     Position r0, Position r1, const Direction& u, vector<int>& 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] current_bin Current bin number
+  //! \param[in] r Position of the particle
+  //! \param[in] u Direction of the particle
+  //! \param[out] bin_next Next bin number
+  //! \return Distance to the next boundary
+  virtual double distance_to_next_boundary(
+    int current_bin, Position r, Direction u, int& bin_next) const = 0;
+
   //! Get bin at a given position in space
   //
   //! \param[in] r Position to get bin for
@@ -333,6 +346,9 @@ class StructuredMesh : public Mesh {
   void surface_bins_crossed(Position r0, Position r1, const Direction& u,
     vector<int>& bins) const override;
 
+  double distance_to_next_boundary(
+    int current_bin, Position r, Direction u, int& bin_next) const override;
+
   //! Determine which cell or surface bins were crossed by a particle
   //
   //! \param[in] r0 Previous position of the particle
@@ -416,6 +432,17 @@ class StructuredMesh : public Mesh {
   virtual MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
     const Position& r0, const Direction& u, double l) const = 0;
 
+  //! Find the closest distance from a point to the mesh boundaries that are
+  //! aligned with the k direction. The point has to be located outside the
+  //! mesh.
+  //!
+  //! \param[in] k direction index of grid surface
+  //! \param[in] r position, from where to calculate the distance
+  //! \param[in] u direction of flight
+  //! \return distance to the mesh boundary
+  virtual double distance_to_mesh_boundary_from_outside(
+    int k, const Position& r, const Direction& u) const = 0;
+
   //! Get a label for the mesh bin
   std::string bin_label(int bin) const override;
 
@@ -468,6 +495,13 @@ class PeriodicStructuredMesh : public StructuredMesh {
     return r - origin_;
   };
 
+  double distance_to_mesh_boundary_from_outside(
+    int k, const Position& r, const Direction& u) const override
+  {
+    fatal_error("Not implemented");
+    return -1.0;
+  }
+
   // Data members
   Position origin_ {0.0, 0.0, 0.0}; //!< Origin of the mesh
 };
@@ -493,6 +527,9 @@ class RegularMesh : public StructuredMesh {
   MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
     const Position& r0, const Direction& u, double l) const override;
 
+  double distance_to_mesh_boundary_from_outside(
+    int k, const Position& r, const Direction& u) const override;
+
   std::pair<vector<double>, vector<double>> plot(
     Position plot_ll, Position plot_ur) const override;
 
@@ -546,6 +583,9 @@ class RectilinearMesh : public StructuredMesh {
   MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
     const Position& r0, const Direction& u, double l) const override;
 
+  double distance_to_mesh_boundary_from_outside(
+    int k, const Position& r, const Direction& u) const override;
+
   std::pair<vector<double>, vector<double>> plot(
     Position plot_ll, Position plot_ur) const override;
 
@@ -715,6 +755,9 @@ class UnstructuredMesh : public Mesh {
   UnstructuredMesh(pugi::xml_node node);
   UnstructuredMesh(hid_t group);
 
+  double distance_to_next_boundary(
+    int current_bin, Position r, Direction u, int& bin_next) const override;
+
   static const std::string mesh_type;
   virtual std::string get_mesh_type() const override;
 

include/openmc/particle_data.h

@@ -265,6 +265,24 @@ class BoundaryInfo {
     0, 0, 0}; //!< which way lattice indices will change
 };
 
+struct TransportEvent {
+  int event_type = EVENT_UNDEFINED; //!< Type of transport event (crossing
+                                    //!< surface, collision, time cutoff)
+  bool cross_surface_geometry =
+    false; //!< if a surface of the geometry has been crossed during the event
+  bool cross_surface_temperature_field =
+    false; //!< if a surface of the temperature field has been crossed  during
+           //!< the event
+
+  // resets all information
+  void clear()
+  {
+    event_type = EVENT_UNDEFINED;
+    cross_surface_geometry = false;
+    cross_surface_temperature_field = false;
+  }
+};
+
 /*
  * Contains all geometry state information for a particle.
  */
@@ -412,6 +430,12 @@ class GeometryState {
   const double& density_mult() const { return density_mult_; }
   double& density_mult_last() { return density_mult_last_; }
 
+  // Temperature field related information
+  int& tf_bin() { return tf_bin_; }
+  const int& tf_bin() const { return tf_bin_; }
+  int& tf_bin_next() { return tf_bin_next_; }
+  const int& tf_bin_next() const { return tf_bin_next_; }
+
 private:
   int64_t id_ {-1}; //!< Unique ID
 
@@ -437,12 +461,15 @@ class GeometryState {
   int material_ {-1};      //!< index for current material
   int material_last_ {-1}; //!< index for last material
 
-  double sqrtkT_ {-1.0};     //!< sqrt(k_Boltzmann * temperature) in eV
-  double sqrtkT_last_ {0.0}; //!< last temperature
+  double sqrtkT_ {-1.0};      //!< sqrt(k_Boltzmann * temperature) in eV
+  double sqrtkT_last_ {-1.0}; //!< last temperature
 
   double density_mult_ {1.0};      //!< density multiplier
   double density_mult_last_ {1.0}; //!< last density multiplier
 
+  int tf_bin_ = C_NONE;      //!< Current temperature field bin
+  int tf_bin_next_ = C_NONE; //!< Next temperature field bin
+
 #ifdef OPENMC_DAGMC_ENABLED
   moab::DagMC::RayHistory history_;
   Direction last_dir_;
@@ -568,6 +595,8 @@ class ParticleData : public GeometryState {
 
   int64_t n_progeny_ {0};
 
+  TransportEvent next_event_ = TransportEvent();
+
 public:
   //----------------------------------------------------------------------------
   // Constructors
@@ -784,6 +813,10 @@ class ParticleData : public GeometryState {
       d = 0;
     }
   }
+
+  // Next event in transport
+  TransportEvent& next_event() { return next_event_; }
+  const TransportEvent& next_event() const { return next_event_; }
 };
 
 } // namespace openmc

include/openmc/settings.h

@@ -91,6 +91,7 @@ extern bool surf_mcpl_write;         //!< write surface mcpl file?
 extern bool surf_source_read;        //!< read surface source file?
 extern bool survival_biasing;        //!< use survival biasing?
 extern bool survival_normalization;  //!< use survival normalization?
+extern bool temperature_field_on;    //!< Is there a temperature field defined?
 extern bool temperature_multipole;   //!< use multipole data?
 extern "C" bool trigger_on;          //!< tally triggers enabled?
 extern bool trigger_predict;         //!< predict batches for triggers?

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<double> k_generation;
 extern vector<int64_t> work_index;