Address comments on initial PR, all changed other than name of class and the use of python API to combine track outputs

include/openmc/distribution_spatial.h

@@ -112,11 +112,10 @@ private:
   UnstructuredMesh* umesh_ptr_;
   int32_t mesh_map_idx_;
   std::string sample_scheme_;
-  double total_weight_;
+  double total_strength_;
   std::vector<double> mesh_CDF_;
-  std::vector<double> mesh_weights_;
+  std::vector<double> mesh_strengths_;
   int64_t tot_bins_;
-  
 };
 
 //==============================================================================

include/openmc/mesh.h

@@ -468,9 +468,6 @@ public:
   static const std::string mesh_type;
   virtual std::string get_mesh_type() const override;
 
-  //! Sample a tetrahedron for an unstructured mesh
-  virtual Position sample_tet(std::array<Position, 4> coords, uint64_t* seed) const;
-
   // Overridden Methods
   // TODO Position sample(uint64_t* seed) const=0;
 
@@ -545,6 +542,9 @@ protected:
     1.0}; //!< Constant multiplication factor to apply to mesh coordinates
   bool specified_length_multiplier_ {false};
 
+  //! Sample a tetrahedron for an unstructured mesh
+  Position sample_tet(std::array<Position, 4> coords, uint64_t* seed) const;
+
 private:
   //! Setup method for the mesh. Builds data structures,
   //! sets up element mapping, creates bounding boxes, etc.

openmc/settings.py

@@ -978,7 +978,9 @@ class Settings:
         for source in self.source:
             root.append(source.to_xml_element())
             if isinstance(source.space, MeshIndependent):
-                root.append(source.space.mesh.to_xml_element())
+                path = f"./mesh[@id='{source.space.mesh.id}']" 
+                if root.find(path) is None: 
+                    root.append(source.space.mesh.to_xml_element())
 
     def _create_volume_calcs_subelement(self, root):
         for calc in self.volume_calculations:

openmc/stats/multivariate.py

@@ -3,7 +3,6 @@ from collections.abc import Iterable
 from math import pi, cos
 from numbers import Real
 from xml.etree import ElementTree as ET
-from xmlrpc.client import boolean
 
 import numpy as np
 
@@ -629,29 +628,25 @@ class MeshIndependent(Spatial):
 
     Parameters
     ----------
-    elem_weight_scheme : str, optional
-        The scheme for weighting and sampling elements from the mesh. Options are 'file' and 'volume' based weights.
     mesh : openmc.MeshBase
         The mesh instance used for sampling, mesh is written into settings.xml, mesh.id is written into the source distribution
-    weights_from_file : Iterable of Real, optional
+    strengths : Iterable of Real, optional
         A list of values which represent the weights of each element
 
 
     Attributes
     ----------
-    elem_weight_scheme : str, optional
-        Weighting scheme for sampling element from mesh, defaults to volume
     mesh : openmc.MeshBase
         The mesh instance used for sampling, mesh is written into settings.xml, mesh.id is written into the source distribution
-    weights_from_file : Iterable of Real, optional
+    strengths : Iterable of Real, optional
         A list of values which represent the weights of each element
 
     """
 
-    def __init__(self, mesh, weights_from_file=None, volume_weighting=False):
+    def __init__(self, mesh, strengths=None, volume_normalized=False):
         self.mesh = mesh        
-        self.weights_from_file = weights_from_file
-        self.volume_weighting = volume_weighting
+        self.strengths = strengths
+        self.volume_normalized = volume_normalized
 
     @property
     def mesh(self):
@@ -660,34 +655,35 @@ class MeshIndependent(Spatial):
     @mesh.setter
     def mesh(self, mesh):
         if mesh != None:
-            cv.check_type('Unstructured Mesh', mesh, MeshBase)
+            cv.check_type('Unstructured Mesh instance', mesh, MeshBase)
         self._mesh = mesh
 
     @property
-    def volume_weighting(self):
-        return self._volume_weighting
+    def volume_normalized(self):
+        return self._volume_normalized
 
-    @volume_weighting.setter
-    def volume_weighting(self, volume_weighting):
-        cv.check_type('Scheme for sampling an element from the mesh', volume_weighting, bool)
-        self._volume_weighting = volume_weighting
+    @volume_normalized.setter
+    def volume_normalized(self, volume_normalized):
+        cv.check_type('Normalize (multiply) strengths by volume', volume_normalized, bool)
+        self._volume_normalized = volume_normalized
 
     @property
-    def weights_from_file(self):
-        return self._weights_from_file
+    def strengths(self):
+        return self._strengths
 
-    @weights_from_file.setter
-    def weights_from_file(self, given_weights):
-        if given_weights is not None:
-            self._weights_from_file = (np.array(given_weights, dtype=float)).flatten()
+    @strengths.setter
+    def strengths(self, given_strengths):
+        if given_strengths is not None:
+            cv.check_type('strengths array passed in', given_strengths, Iterable, Real)
+            self._strengths = np.array(given_strengths, dtype=float).flatten()
         else:
-            self._weights_from_file = None
+            self._strengths = None
 
     @property
-    def num_weight_bins(self):
-        if self.weights_from_file == None:
-            raise ValueError('Weight bins are not set')
-        return self.weights_from_file.size
+    def num_strength_bins(self):
+        if self.strengths == None:
+            raise ValueError('Strengths are not set')
+        return self.strengths.size
                 
     def to_xml_element(self):
         """Return XML representation of the spatial distribution
@@ -701,11 +697,11 @@ class MeshIndependent(Spatial):
         element = ET.Element('space')
         element.set('type', 'mesh')
         element.set("mesh_id", str(self.mesh.id))
-        element.set("volume_weighting", str(self.volume_weighting))
+        element.set("volume_normalized", str(self.volume_normalized))
 
-        if self.weights_from_file is not None:
-            subelement = ET.SubElement(element, 'weights_from_file')
-            subelement.text = ' '.join(str(e) for e in self.weights_from_file)
+        if self.strengths is not None:
+            subelement = ET.SubElement(element, 'strengths')
+            subelement.text = ' '.join(str(e) for e in self.strengths)
 
         return element
 
@@ -726,13 +722,13 @@ class MeshIndependent(Spatial):
         """
 
         mesh_id = int(elem.get('mesh_id'))
-        volume_weighting = elem.get("volume_weighting")
-        volume_weighting = get_text(elem, 'volume_weighting').lower() == 'true'
-        if elem.get('weights_from_file') is not None:
-            weights_from_file = [float(b) for b in get_text(elem, 'weights_from_file').split()]
+        volume_normalized = elem.get("volume_normalized")
+        volume_normalized = get_text(elem, 'volume_normalized').lower() == 'true'
+        if elem.get('strengths') is not None:
+            strengths = [float(b) for b in get_text(elem, 'strengths').split()]
         else:
-            weights_from_file = None
-        return cls(volume_weighting, mesh_id, weights_from_file)
+            strengths = None
+        return cls(volume_normalized, mesh_id, strengths)
 
 
 class Box(Spatial):

src/distribution_spatial.cpp

@@ -6,9 +6,6 @@
 #include "openmc/mesh.h"
 #include "openmc/search.h"
 
-#include <iostream>
-#include <fstream>
-
 namespace openmc {
 
 //==============================================================================
@@ -202,48 +199,48 @@ MeshIndependent::MeshIndependent(pugi::xml_node node)
   umesh_ptr_ = dynamic_cast<UnstructuredMesh*>(mesh_ptr.get());
   if (!umesh_ptr_) {fatal_error("Mesh passed to spatial distribution is not an unstructured mesh object"); }
 
-  // Create CDF based on weighting scheme
+  // Initialize arrays for CDF creation
   tot_bins_ = umesh_ptr_->n_bins();
-  std::vector<double> weights = {};
-  weights.resize(tot_bins_);
-  double temp_total_weight = 0.0;
+  std::vector<double> strengths = {};
+  strengths.resize(tot_bins_);
+  double temp_total_strength = 0.0;
   mesh_CDF_.resize(tot_bins_+1);
   mesh_CDF_[0] = {0.0};
-  total_weight_ = 0.0;
-  mesh_weights_.resize(tot_bins_);
+  total_strength_ = 0.0;
+  mesh_strengths_.resize(tot_bins_);
 
   // Create cdfs for sampling for an element over a mesh
   // Volume scheme is weighted by the volume of each tet
   // File scheme is weighted by an array given in the xml file
 
-  if (check_for_node(node, "weights_from_file")) {
-    mesh_weights_ = get_node_array<double>(node, "weights_from_file");
-    if (mesh_weights_.size() != tot_bins_){
-      write_message("WARNING: The size of the weights array from the xml file does not equal the number of elements in the mesh.");
-    } if (get_node_value_bool(node, "volume_weighting")){
+  if (check_for_node(node, "strengths")) {
+    mesh_strengths_ = get_node_array<double>(node, "strengths");
+    if (mesh_strengths_.size() != tot_bins_){
+      fatal_error("The size of the strengths array from the xml file does not equal the number of elements in the mesh.");
+    } if (get_node_value_bool(node, "volume_normalized")){
       for (int i = 0; i < tot_bins_; i++){
-        weights[i] = mesh_weights_[i]*umesh_ptr_->volume(i);
+        strengths[i] = mesh_strengths_[i]*umesh_ptr_->volume(i);
       } 
-    } else if (!get_node_value_bool(node, "volume_weighting")){
+    } else if (!get_node_value_bool(node, "volume_normalized")){
       for (int i = 0; i < tot_bins_; i++){
-        weights[i] = mesh_weights_[i];  
+        strengths[i] = mesh_strengths_[i];  
       }
     }
-  } else if (get_node_value_bool(node, "volume_weighting")){
+  } else if (get_node_value_bool(node, "volume_normalized")){
     for (int i = 0; i<tot_bins_; i++){
-      weights[i] = umesh_ptr_->volume(i);
+      strengths[i] = umesh_ptr_->volume(i);
     }
   } else {
     for (int i = 0; i<tot_bins_; i++){
-      weights[i] = 1;
+      strengths[i] = 1;
     }
   }
   for (int i = 0; i<tot_bins_; i++){
-    temp_total_weight = temp_total_weight + weights[i];
+    temp_total_strength = temp_total_strength + strengths[i];
   }
-  total_weight_ = temp_total_weight;
+  total_strength_ = temp_total_strength;
   for (int i = 0; i<tot_bins_; i++){
-    mesh_CDF_[i+1] = mesh_CDF_[i] + weights[i]/total_weight_;
+    mesh_CDF_[i+1] = mesh_CDF_[i] + strengths[i]/total_strength_;
   }
 }
 

src/mesh.cpp

@@ -358,7 +358,7 @@ StructuredMesh::MeshIndex StructuredMesh::get_indices_from_bin(int bin) const
   return ijk;
 }
 
-Position StructuredMesh::sample(uint64_t* seed, int32_t tet_bin) const {
+Position StructuredMesh::sample(uint64_t* seed, int32_t bin) const {
   fatal_error("Position sampling on structured meshes is not yet implemented");
 }
 
@@ -2085,9 +2085,9 @@ std::string MOABMesh::library() const
 }
 
 // Sample position within a tet for MOAB type tets
-Position MOABMesh::sample(uint64_t* seed, int32_t tet_bin) const {
+Position MOABMesh::sample(uint64_t* seed, int32_t bin) const {
 
-  moab::EntityHandle tet_ent = get_ent_handle_from_bin(tet_bin);
+  moab::EntityHandle tet_ent = get_ent_handle_from_bin(bin);
 
   // Get vertex coordinates for MOAB tet
   vector<moab::EntityHandle> conn1;
@@ -2103,7 +2103,7 @@ Position MOABMesh::sample(uint64_t* seed, int32_t tet_bin) const {
 
   std::array<Position, 4> tet_verts;
   for (int i = 0; i < 4; i++) {
-    tet_verts[i] =  {p[i][0], p[i][1], p[i][2]};
+    tet_verts[i] = {p[i][0], p[i][1], p[i][2]};
   }
   // Samples position within tet using Barycentric stuff  
   Position r = this->sample_tet(tet_verts, seed);
@@ -2565,8 +2565,8 @@ void LibMesh::initialize()
 }
 
 // Sample position within a tet for LibMesh type tets
-Position LibMesh::sample(uint64_t* seed, int32_t tet_bin) const {
-  const auto& elem = get_element_from_bin(tet_bin);
+Position LibMesh::sample(uint64_t* seed, int32_t bin) const {
+  const auto& elem = get_element_from_bin(bin);
   // Get tet vertex coordinates from LibMesh
   std::array<Position, 4> tet_verts;
   for (int i = 0; i < elem.n_nodes(); i++) {

tests/regression_tests/unstructured_mesh/source_sampling/test.py

@@ -16,7 +16,6 @@ TETS_PER_VOXEL = 12
 
 class UnstructuredMeshSourceTest(PyAPITestHarness):
     def __init__(self, statepoint_name, model, inputs_true, schemes):
-
         super().__init__(statepoint_name, model, inputs_true)
         self.schemes = schemes
 
@@ -31,60 +30,53 @@ class UnstructuredMeshSourceTest(PyAPITestHarness):
         openmc.run(**kwargs)
 
     def _compare_results(self):
-
+        # There are 10000 particles and 1000 hexes, this leads to the average
+        # shown below
         average_in_hex = 10.0
 
-
         call(['../../../../scripts/openmc-track-combine', '-o', 'tracks.h5'] +
                 glob.glob('tracks_p*.h5'))
 
         # loop over the tracks and get data
         tracks = openmc.Tracks(filepath='tracks.h5')
         tracks_born = np.empty((len(tracks), 1))
-        decimals = 0
 
         instances = np.zeros(1000)
 
         for i in range(0, len(tracks)):
             tracks_born[i] = tracks[i].particle_tracks[0].states['cell_id'][0]
-            instances[int(tracks_born[i])-1] = instances[int(tracks_born[i])-1]+1
+            instances[int(tracks_born[i])-1] += 1
 
         if self.schemes == "file":
             assert(instances[0] > 0 and instances[1] > 0)
             assert(instances[0] > instances[1])
-            
-            for i in range(0, len(instances)):
-                if i != 0 and i != 1:
-                    assert(instances[i] == 0)
+
+            for i in range(2, len(instances)):
+                assert(instances[i] == 0)
 
         else:
             assert(np.average(instances) == average_in_hex)
-            assert(np.std(instances) < np.average(instances))  
+            assert(np.std(instances) < np.average(instances))
             assert(np.amax(instances) < np.average(instances)+6*np.std(instances))
-            
+
 
     def _cleanup(self):
         super()._cleanup()
-        output = glob.glob('tally*.vtk')
+        output = glob.glob('track*.h5')
         output += glob.glob('tally*.e')
         for f in output:
             if os.path.exists(f):
                 os.remove(f)
 
-    
-
-
 param_values = (['libmesh', 'moab'], # mesh libraries
                 ['volume', 'file']) # Element weighting schemes
 
 test_cases = []
-# for i, (lib, holes) in enumerate(product(*param_values)):
 for i, (lib, schemes) in enumerate(product(*param_values)):
     test_cases.append({'library' : lib,
                        'schemes' : schemes,
                        'inputs_true' : 'inputs_true{}.dat'.format(i)})
 
-
 @pytest.mark.parametrize("test_opts", test_cases)
 def test_unstructured_mesh(test_opts):
 
@@ -110,19 +102,20 @@ def test_unstructured_mesh(test_opts):
 
     ### Geometry ###
     dimen = 10
-    size_hex = 20.0/dimen
+    size_hex = 20.0 / dimen
 
     ### Geometry ###
     cell = np.empty((dimen, dimen, dimen), dtype=object)
-    surfaces = np.empty((dimen+1, 3), dtype=object)
+    surfaces = np.empty((dimen + 1, 3), dtype=object)
 
     geometry = openmc.Geometry()
     universe = openmc.Universe(universe_id=1, name="Contains all hexes")
 
     for i in range(0,dimen+1):
-        surfaces[i][0] = openmc.XPlane(-10.0+i*size_hex, name="X plane at "+str(-10.0+i*size_hex))
-        surfaces[i][1] = openmc.YPlane(-10.0+i*size_hex, name="Y plane at "+str(-10.0+i*size_hex))
-        surfaces[i][2] = openmc.ZPlane(-10.0+i*size_hex, name="Z plane at "+str(-10.0+i*size_hex))
+        coord = -10.0 + i * size_hex
+        surfaces[i][0] = openmc.XPlane(coord, name="X plane at "+str(coord))
+        surfaces[i][1] = openmc.YPlane(coord, name="Y plane at "+str(coord))
+        surfaces[i][2] = openmc.ZPlane(coord, name="Z plane at "+str(coord))
 
         surfaces[i][0].boundary_type = 'vacuum'
         surfaces[i][1].boundary_type = 'vacuum'
@@ -131,11 +124,11 @@ def test_unstructured_mesh(test_opts):
     for k in range(0,dimen):
         for j in range(0,dimen):
             for i in range(0,dimen):
-                cell[i][j][k] = openmc.Cell(name=("x " + str(i) +" y " + str(j) +" z " + str(k)))
+                cell[i][j][k] = openmc.Cell(name=("x = {}, y = {}, z = {}".format(i,j,k)))
                 cell[i][j][k].region = +surfaces[i][0] & -surfaces[i+1][0] & \
                                     +surfaces[j][1] & -surfaces[j+1][1] & \
                                     +surfaces[k][2] & -surfaces[k+1][2]
-                cell[i][j][k].fill = water_mat
+                cell[i][j][k].fill = None
                 universe.add_cell(cell[i][j][k])
 
     geometry = openmc.Geometry(universe)
@@ -161,19 +154,17 @@ def test_unstructured_mesh(test_opts):
     settings.particles = 5000
     settings.batches = 2
 
-    # settings.create_fission_neutrons = False
-    settings.tracks = [(1, 1, 1)]
     settings.max_tracks = 10000
 
     # source setup
     if test_opts['schemes'] == 'volume':
-        space = openmc.stats.MeshIndependent(volume_weighting=True, mesh=uscd_mesh)
+        space = openmc.stats.MeshIndependent(volume_normalized=True, mesh=uscd_mesh)
     elif test_opts['schemes'] == 'file':
         array = np.zeros(12000)
         for i in range(0, 12):
             array[i] = 10
             array[i+12] = 2
-        space = openmc.stats.MeshIndependent(volume_weighting=False, weights_from_file=array, mesh=uscd_mesh)
+        space = openmc.stats.MeshIndependent(volume_normalized=False, strengths=array, mesh=uscd_mesh)
 
     energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0])
     source = openmc.Source(space=space, energy=energy)
@@ -185,7 +176,7 @@ def test_unstructured_mesh(test_opts):
                                settings=settings)
 
     harness = UnstructuredMeshSourceTest('statepoint.2.h5',
-                                   model,
-                                   test_opts['inputs_true'],
-                                   test_opts['schemes'])
+                                         model,
+                                         test_opts['inputs_true'],
+                                         test_opts['schemes'])
     harness.main()