Specify temperature from a field (structured mesh only)

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

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
 -----------------
 

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

@@ -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

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<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 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<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 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

include/openmc/mesh.h

@@ -173,6 +173,16 @@ 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] 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<int>& 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;
 

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();

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

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?

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;
 

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 *

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)