New Feature Development: The Implementation of chord length sampling (CLS) method for stochastic media( #3286 )

CMakeLists.txt

@@ -404,6 +404,7 @@ list(APPEND libopenmc_SOURCES
   src/simulation.cpp
   src/source.cpp
   src/state_point.cpp
+  src/stochastic_media.cpp
   src/string_utils.cpp
   src/summary.cpp
   src/surface.cpp

include/openmc/cell.h

@@ -26,7 +26,7 @@ namespace openmc {
 // Constants
 //==============================================================================
 
-enum class Fill { MATERIAL, UNIVERSE, LATTICE };
+enum class Fill { MATERIAL, UNIVERSE, LATTICE, STOCHASTIC_MEDIA };
 
 constexpr int32_t OP_LEFT_PAREN {std::numeric_limits<int32_t>::max()};
 constexpr int32_t OP_RIGHT_PAREN {std::numeric_limits<int32_t>::max() - 1};

include/openmc/constants.h

@@ -355,6 +355,12 @@ enum class GeometryType { CSG, DAG };
 // representations. This value represents no surface.
 constexpr int32_t SURFACE_NONE {0};
 
+enum class ParticleStatus {
+  OUTSIDE,   //  Particle is outside the cell containing stochastic media
+  IN_MATRIX, //  Particle is inside the matrix but not in the stochastic media
+  IN_STOCHASTIC_MEDIA, //  Particle is inside the stochastic media
+};
+
 } // namespace openmc
 
 #endif // OPENMC_CONSTANTS_H

include/openmc/particle.h

@@ -79,6 +79,8 @@ class Particle : public ParticleData {
   //! Cross a surface and handle boundary conditions
   void cross_surface(const Surface& surf);
 
+  void cross_surface_in_stochmedia();
+
   //! Cross a vacuum boundary condition.
   //
   //! \param surf The surface (with the vacuum boundary condition) that the

include/openmc/particle_data.h

@@ -1,3 +1,3 @@
 #ifndef OPENMC_PARTICLE_DATA_H
 #define OPENMC_PARTICLE_DATA_H
 
@@ -192,7 +192,8 @@
     SURFACE_NONE}; //!< surface token, non-zero if boundary is surface
   int coord_level; //!< coordinate level after crossing boundary
   array<int, 3>
-    lattice_translation {}; //!< which way lattice indices will change
+    lattice_translation {};           //!< which way lattice indices will change
+  bool if_stochastic_surface {false}; //!< is the surface in stochastic media?
 
   void reset()
   {
@@ -330,14 +331,18 @@
   // Boundary information
   BoundaryInfo& boundary() { return boundary_; }
 
+  // Particle status
+  ParticleStatus& status() { return status_; }
+  const ParticleStatus& status() const { return status_; }
+
 #ifdef DAGMC
   // DagMC state variables
   moab::DagMC::RayHistory& history() { return history_; }
   Direction& last_dir() { return last_dir_; }
 #endif
 
   // material of current and last cell
-  int& material() { return material_; }
+  int& material() { return material_; } 
   const int& material() const { return material_; }
   int& material_last() { return material_last_; }
   const int& material_last() const { return material_last_; }
@@ -369,6 +374,9 @@
 
   BoundaryInfo boundary_; //!< Info about the next intersection
 
+  ParticleStatus status_ {
+    ParticleStatus::OUTSIDE}; // The status of the particle in stochastic media
+
   int material_ {-1};      //!< index for current material
   int material_last_ {-1}; //!< index for last material
 

include/openmc/stochastic_media.h

@@ -0,0 +1,104 @@
+#include <string>
+#include <unordered_map>
+
+#include "pugixml.hpp"
+#include "xtensor/xtensor.hpp"
+
+#include "openmc/constants.h"
+#include "openmc/memory.h" // for unique_ptr
+#include "openmc/vector.h"
+
+namespace openmc {
+//==============================================================================
+// Global variables
+//==============================================================================
+class Stochastic_Media;
+
+namespace model {
+extern std::unordered_map<int32_t, int32_t> stochastic_media_map;
+extern vector<unique_ptr<Stochastic_Media>> stochastic_media;
+} // namespace model
+//==============================================================================
+class Stochastic_Media {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors, factory functions
+
+  explicit Stochastic_Media(pugi::xml_node cell_node);
+  Stochastic_Media() {};
+  virtual ~Stochastic_Media();
+
+  // Accessors
+  //! Get name
+  //! \return Stochastic Media name
+  const std::string& name() const { return name_; }
+
+  //! Set name
+  void set_name(const std::string& name) { name_ = name; }
+
+  //! Get ID of stochastic media
+  //! \return ID of stochastic media
+  int32_t id() const { return id_; }
+
+  //! Assign a unique ID to the stochastic media
+  //! \param[in] Unique ID to assign. A value of -1 indicates that an ID
+  //!   should be automatically assigned.
+  void set_id(int32_t id_);
+
+  //! Get the particle radius
+  double radius() const { return radius_; }
+
+  //! Get the packing fraction of the particle material
+  double pf() const { return pf_; }
+
+  uint64_t* current_seed() { return &seed_; }
+
+  //! Get the material of the particle
+  int32_t particle_mat(int32_t i) const { return particle_mat_[i]; }
+  vector<int32_t> particle_mat() const { return particle_mat_; }
+
+  //! Get the material of the matrix
+  int32_t matrix_mat() const { return matrix_mat_; }
+
+  //----------------------------------------------------------------------------
+  // Data
+  int32_t id_ {C_NONE}; //!< Unique ID
+  std::string name_;    //!< Name of stochastic_media
+  // Parameters: the  particle radius, now this method only support sphere
+  double radius_;
+  // The packing fraction of the particle material
+  double pf_;
+  // The material of the particle
+  vector<int32_t> particle_mat_;
+  // The material of the matrix
+  int32_t matrix_mat_;
+  //----------------------------------------------------------------------------
+
+  virtual void sample_material(GeometryState& p) = 0;
+  virtual void adjust_indices();
+
+protected:
+  // Protected data members
+  int64_t index_;
+  uint64_t seed_ {init_seed(id_, STREAM_SOURCE)};
+};
+
+class CLS_Media : public Stochastic_Media {
+public:
+  explicit CLS_Media(pugi::xml_node cell_node);
+  CLS_Media() {};
+
+  void sample_material(GeometryState& p) override;
+};
+//==============================================================================
+// Non-member functions
+//==============================================================================
+//! Read stochastic media data XML node
+//! \param[in] root node of stochastic_media XML element
+void read_stochastic_media(pugi::xml_node root);
+
+double distance_to_stochamedia(Particle& p);
+
+void free_memory_stochastic_media();
+
+} // namespace openmc

src/finalize.cpp

@@ -20,6 +20,7 @@
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
 #include "openmc/source.h"
+#include "openmc/stochastic_media.h"
 #include "openmc/surface.h"
 #include "openmc/tallies/tally.h"
 #include "openmc/thermal.h"
@@ -36,6 +37,7 @@ void free_memory()
   free_memory_geometry();
   free_memory_surfaces();
   free_memory_material();
+  free_memory_stochastic_media();
   free_memory_volume();
   free_memory_simulation();
   free_memory_photon();

src/geometry.cpp

@@ -12,6 +12,7 @@
 #include "openmc/simulation.h"
 #include "openmc/string_utils.h"
 #include "openmc/surface.h"
+#include <openmc/stochastic_media.h>
 
 namespace openmc {
 
@@ -177,9 +178,16 @@ bool find_cell_inner(
       p.material() = c.material(p.cell_instance());
       p.sqrtkT_last() = p.sqrtkT();
       p.sqrtkT() = c.sqrtkT(p.cell_instance());
+      p.status() = ParticleStatus::OUTSIDE;
 
       return true;
 
+    } else if (c.type_ == Fill::STOCHASTIC_MEDIA) {
+      //========================================================================
+      //! Found stochastic media cell, means this is the lowest coord level.
+      model::stochastic_media[c.fill_]->sample_material(p);
+      return true;
+
     } else if (c.type_ == Fill::UNIVERSE) {
       //========================================================================
       //! Found a lower universe, update this coord level then search the next.

src/geometry_aux.cpp

@@ -17,6 +17,7 @@
 #include "openmc/lattice.h"
 #include "openmc/material.h"
 #include "openmc/settings.h"
+#include "openmc/stochastic_media.h"
 #include "openmc/surface.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/filter_cell_instance.h"
@@ -66,10 +67,11 @@ void read_geometry_xml()
 
 void read_geometry_xml(pugi::xml_node root)
 {
-  // Read surfaces, cells, lattice
+  // Read surfaces, cells, lattice, and stochastic media
   read_surfaces(root);
   read_cells(root);
   read_lattices(root);
+  read_stochastic_media(root);
 
   // Check to make sure a boundary condition was applied to at least one
   // surface
@@ -103,15 +105,19 @@ void adjust_indices()
       int32_t id = c->fill_;
       auto search_univ = model::universe_map.find(id);
       auto search_lat = model::lattice_map.find(id);
+      auto search_media = model::stochastic_media_map.find(id);
       if (search_univ != model::universe_map.end()) {
         c->type_ = Fill::UNIVERSE;
         c->fill_ = search_univ->second;
       } else if (search_lat != model::lattice_map.end()) {
         c->type_ = Fill::LATTICE;
         c->fill_ = search_lat->second;
+      } else if (search_media != model::stochastic_media_map.end()) {
+        c->type_ = Fill::STOCHASTIC_MEDIA;
+        c->fill_ = search_media->second;
       } else {
-        fatal_error(fmt::format("Specified fill {} on cell {} is neither a "
-                                "universe nor a lattice.",
+        fatal_error(fmt::format("Specified fill {} on cell {} must be one of  "
+                                "universe  lattice and stochastic media.",
           id, c->id_));
       }
     } else {
@@ -146,6 +152,10 @@ void adjust_indices()
   for (auto& l : model::lattices) {
     l->adjust_indices();
   }
+
+  for (auto& media : model::stochastic_media) {
+    media->adjust_indices();
+  }
 }
 
 //==============================================================================
@@ -187,8 +197,16 @@ void assign_temperatures()
 {
   for (auto& c : model::cells) {
     // Ignore non-material cells and cells with defined temperature.
-    if (c->material_.size() == 0)
-      continue;
+    if (c->material_.size() == 0) {
+      if (c->type_ != Fill::STOCHASTIC_MEDIA)
+        continue;
+      c->sqrtkT_.reserve(1);
+
+      auto fill_idx = model::stochastic_media_map[c->fill_];
+      auto mat_id {model::stochastic_media[fill_idx]->particle_mat()[0]};
+      const auto& mat {model::materials[mat_id]};
+      c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * mat->temperature()));
+    }
     if (c->sqrtkT_.size() > 0)
       continue;
 
@@ -639,6 +657,9 @@ void free_memory_geometry()
   model::lattices.clear();
   model::lattice_map.clear();
 
+  model::stochastic_media.clear();
+  model::stochastic_media_map.clear();
+
   model::overlap_check_count.clear();
 }
 

src/particle.cpp

@@ -25,6 +25,7 @@
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
 #include "openmc/source.h"
+#include "openmc/stochastic_media.h"
 #include "openmc/surface.h"
 #include "openmc/tallies/derivative.h"
 #include "openmc/tallies/tally.h"
@@ -229,6 +230,14 @@ void Particle::event_advance()
   // Select smaller of the two distances
   double distance = std::min(boundary().distance, collision_distance());
 
+  if (this->status() != ParticleStatus::OUTSIDE) {
+    double stocha_media_distance = distance_to_stochamedia(*this);
+    if (stocha_media_distance < distance) {
+      distance = stocha_media_distance;
+      boundary().if_stochastic_surface = true;
+    }
+  }
+
   // Advance particle in space and time
   // Short-term solution until the surface source is revised and we can use
   // this->move_distance(distance)
@@ -607,6 +616,50 @@ void Particle::cross_surface(const Surface& surf)
   }
 }
 
+void Particle::cross_surface_in_stochmedia()
+{
+  surface() = SURFACE_NONE;
+  bool verbose = settings::verbosity >= 10 || trace();
+
+  double i_cell = this->lowest_coord().cell;
+  Cell& c {*model::cells[i_cell]};
+
+  this->cell_instance() = 0;
+  // Find the distribcell instance number.
+  if (c.distribcell_index_ >= 0) {
+    this->cell_instance() = cell_instance_at_level(*this, this->n_coord() - 1);
+  }
+
+  // If the Particle is in a cell containing a stochastic medium, update the
+  // status of particle based on the current stochastic medium.
+  Stochastic_Media& media {*model::stochastic_media[c.fill_]};
+  int32_t material;
+  if (this->status() == ParticleStatus::IN_STOCHASTIC_MEDIA) {
+
+    material = media.matrix_mat();
+    this->status() = ParticleStatus::IN_MATRIX;
+    if (verbose) {
+      write_message(1,
+        "    Crossing surface in stochasitic media {} from particle to matrix",
+        media.id_);
+    }
+
+  } else if (this->status() == ParticleStatus::IN_MATRIX) {
+    // short time implementation, only one particle type in matrix
+    material = media.particle_mat(0);
+    this->status() = ParticleStatus::IN_STOCHASTIC_MEDIA;
+    if (verbose) {
+      write_message(1,
+        "    Crossing surface in stochasitic media {} from matrix to particle",
+        media.id_);
+    }
+  }
+  this->material_last() = this->material();
+  this->material() = material;
+  this->sqrtkT_last() = this->sqrtkT();
+  this->sqrtkT() = c.sqrtkT(this->cell_instance());
+}
+
 void Particle::cross_vacuum_bc(const Surface& surf)
 {
   // Score any surface current tallies -- note that the particle is moved