Expose Material::depletable in the CAPI (openmc-dev#2843)

Co-authored-by: Paul Romano <[email protected]>

include/openmc/capi.h

@@ -93,6 +93,8 @@ int openmc_material_set_id(int32_t index, int32_t id);
 int openmc_material_get_name(int32_t index, const char** name);
 int openmc_material_set_name(int32_t index, const char* name);
 int openmc_material_set_volume(int32_t index, double volume);
+int openmc_material_get_depletable(int32_t index, bool* depletable);
+int openmc_material_set_depletable(int32_t index, bool depletable);
 int openmc_material_filter_get_bins(
   int32_t index, const int32_t** bins, size_t* n);
 int openmc_material_filter_set_bins(

include/openmc/material.h

@@ -94,7 +94,7 @@ class Material {
   //  which will get auto-assigned to the next available ID. After creating
   //  the new material, it is added to openmc::model::materials.
   //! \return reference to the cloned material
-  Material & clone();
+  Material& clone();
 
   //----------------------------------------------------------------------------
   // Accessors
@@ -149,6 +149,7 @@ class Material {
   //! Get whether material is fissionable
   //! \return Whether material is fissionable
   bool fissionable() const { return fissionable_; }
+  bool& fissionable() { return fissionable_; }
 
   //! Get volume of material
   //! \return Volume in [cm^3]
@@ -158,6 +159,10 @@ class Material {
   //! \return Temperature in [K]
   double temperature() const;
 
+  //! Whether or not the material is depletable
+  bool depletable() const { return depletable_; }
+  bool& depletable() { return depletable_; }
+
   //! Get pointer to NCrystal material object
   //! \return Pointer to NCrystal material object
   const NCrystalMat& ncrystal_mat() const { return ncrystal_mat_; };
@@ -173,11 +178,8 @@ class Material {
   double density_;                      //!< Total atom density in [atom/b-cm]
   double density_gpcc_;                 //!< Total atom density in [g/cm^3]
   double volume_ {-1.0};                //!< Volume in [cm^3]
-  bool fissionable_ {
-    false};                 //!< Does this material contain fissionable nuclides
-  bool depletable_ {false}; //!< Is the material depletable?
-  vector<bool> p0_;         //!< Indicate which nuclides are to be treated with
-                            //!< iso-in-lab scattering
+  vector<bool> p0_; //!< Indicate which nuclides are to be treated with
+                    //!< iso-in-lab scattering
 
   // To improve performance of tallying, we store an array (direct address
   // table) that indicates for each nuclide in data::nuclides the index of the
@@ -210,6 +212,9 @@ class Material {
   // Private data members
   gsl::index index_;
 
+  bool depletable_ {false}; //!< Is the material depletable?
+  bool fissionable_ {
+    false}; //!< Does this material contain fissionable nuclides
   //! \brief Default temperature for cells containing this material.
   //!
   //! A negative value indicates no default temperature was specified.

openmc/lib/material.py

@@ -1,5 +1,5 @@
 from collections.abc import Mapping
-from ctypes import c_int, c_int32, c_double, c_char_p, POINTER, c_size_t
+from ctypes import c_bool, c_int, c_int32, c_double, c_char_p, POINTER, c_size_t
 from weakref import WeakValueDictionary
 
 import numpy as np
@@ -60,6 +60,12 @@
 _dll.openmc_material_set_volume.argtypes = [c_int32, c_double]
 _dll.openmc_material_set_volume.restype = c_int
 _dll.openmc_material_set_volume.errcheck = _error_handler
+_dll.openmc_material_get_depletable.argtypes = [c_int32, POINTER(c_bool)]
+_dll.openmc_material_get_depletable.restype = c_int
+_dll.openmc_material_get_depletable.errcheck = _error_handler
+_dll.openmc_material_set_depletable.argtypes = [c_int32, c_bool]
+_dll.openmc_material_set_depletable.restype = c_int
+_dll.openmc_material_set_depletable.errcheck = _error_handler
 _dll.n_materials.argtypes = []
 _dll.n_materials.restype = c_size_t
 
@@ -89,6 +95,8 @@ class Material(_FortranObjectWithID):
         List of nuclides in the material
     densities : numpy.ndarray
         Array of densities in atom/b-cm
+    depletable : bool
+        Whether this material is marked as depletable
     name : str
         Name of the material
     temperature : float
@@ -169,6 +177,16 @@ def volume(self):
     def volume(self, volume):
         _dll.openmc_material_set_volume(self._index, volume)
 
+    @property
+    def depletable(self):
+        depletable = c_bool()
+        _dll.openmc_material_get_depletable(self._index, depletable)
+        return depletable.value
+
+    @depletable.setter
+    def depletable(self, depletable):
+        _dll.openmc_material_set_depletable(self._index, depletable)
+
     @property
     def nuclides(self):
         return self._get_densities()[0]

src/event.cpp

@@ -51,7 +51,7 @@ void dispatch_xs_event(int64_t buffer_idx)
 {
   Particle& p = simulation::particles[buffer_idx];
   if (p.material() == MATERIAL_VOID ||
-      !model::materials[p.material()]->fissionable_) {
+      !model::materials[p.material()]->fissionable()) {
     simulation::calculate_nonfuel_xs_queue.thread_safe_append({p, buffer_idx});
   } else {
     simulation::calculate_fuel_xs_queue.thread_safe_append({p, buffer_idx});

src/material.cpp

@@ -372,8 +372,8 @@ Material& Material::clone()
   mat->density_ = density_;
   mat->density_gpcc_ = density_gpcc_;
   mat->volume_ = volume_;
-  mat->fissionable_ = fissionable_;
-  mat->depletable_ = depletable_;
+  mat->fissionable() = fissionable_;
+  mat->depletable() = depletable_;
   mat->p0_ = p0_;
   mat->mat_nuclide_index_ = mat_nuclide_index_;
   mat->thermal_tables_ = thermal_tables_;
@@ -1068,7 +1068,7 @@ void Material::to_hdf5(hid_t group) const
 {
   hid_t material_group = create_group(group, "material " + std::to_string(id_));
 
-  write_attribute(material_group, "depletable", static_cast<int>(depletable_));
+  write_attribute(material_group, "depletable", static_cast<int>(depletable()));
   if (volume_ > 0.0) {
     write_attribute(material_group, "volume", volume_);
   }
@@ -1550,6 +1550,30 @@ extern "C" int openmc_material_set_volume(int32_t index, double volume)
   }
 }
 
+extern "C" int openmc_material_get_depletable(int32_t index, bool* depletable)
+{
+  if (index < 0 || index >= model::materials.size()) {
+    set_errmsg("Index in materials array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  *depletable = model::materials[index]->depletable();
+
+  return 0;
+}
+
+extern "C" int openmc_material_set_depletable(int32_t index, bool depletable)
+{
+  if (index < 0 || index >= model::materials.size()) {
+    set_errmsg("Index in materials array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  model::materials[index]->depletable() = depletable;
+
+  return 0;
+}
+
 extern "C" int openmc_extend_materials(
   int32_t n, int32_t* index_start, int32_t* index_end)
 {

src/mgxs_interface.cpp

@@ -284,7 +284,7 @@ void mark_fissionable_mgxs_materials()
   for (const auto& mat : model::materials) {
     for (int i_nuc : mat->nuclide_) {
       if (data::mg.nuclides_[i_nuc].fissionable) {
-        mat->fissionable_ = true;
+        mat->fissionable() = true;
       }
     }
   }

src/physics_mg.cpp

@@ -43,7 +43,7 @@ void sample_reaction(Particle& p)
   // change when sampling fission sites. The following block handles all
   // absorption (including fission)
 
-  if (model::materials[p.material()]->fissionable_) {
+  if (model::materials[p.material()]->fissionable()) {
     if (settings::run_mode == RunMode::EIGENVALUE ||
         (settings::run_mode == RunMode::FIXED_SOURCE &&
           settings::create_fission_neutrons)) {

src/source.cpp

@@ -201,7 +201,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
           if (mat_index == MATERIAL_VOID) {
             found = false;
           } else {
-            found = model::materials[mat_index]->fissionable_;
+            found = model::materials[mat_index]->fissionable();
           }
         }
       }

tests/unit_tests/test_lib.py

@@ -206,6 +206,10 @@ def test_material(lib_init):
     m.name = "Not hot borated water"
     assert m.name == "Not hot borated water"
 
+    assert m.depletable == False
+    m.depletable = True
+    assert m.depletable == True
+
 
 def test_properties_density(lib_init):
     m = openmc.lib.materials[1]