Allow mesh material volume calculations outside model geometry

include/openmc/mesh.h

@@ -252,6 +252,9 @@ class Mesh {
   virtual std::string get_mesh_type() const = 0;
 
   //! Determine volume of materials within each mesh element
+  //!
+  //! Portions of mesh elements outside the model geometry are treated as void.
+  //! Universe fills within the model must still define all enclosed space.
   //
   //! \param[in] nx Number of samples in x direction
   //! \param[in] ny Number of samples in y direction
@@ -264,6 +267,9 @@ class Mesh {
     int32_t* materials, double* volumes) const;
 
   //! Determine volume and bounding boxes of materials within each mesh element
+  //!
+  //! Portions of mesh elements outside the model geometry are treated as void.
+  //! Universe fills within the model must still define all enclosed space.
   //
   //! \param[in] nx Number of samples in x direction
   //! \param[in] ny Number of samples in y direction

openmc/lib/mesh.py

@@ -196,7 +196,10 @@ def material_volumes(
         This method works by raytracing repeatedly through the mesh to count the
         estimated volume of each material in all mesh elements. Three sets of
         rays are used: one set parallel to the x-axis, one parallel to the
-        y-axis, and one parallel to the z-axis.
+        y-axis, and one parallel to the z-axis. Regions of the mesh that are
+        outside the model geometry are treated as void, equivalent to a cell
+        with no material. Universe fills within the model must still define all
+        enclosed space.
 
         .. versionadded:: 0.15.0
 

openmc/mesh.py

@@ -451,7 +451,10 @@ def material_volumes(
         This method works by raytracing repeatedly through the mesh to count the
         estimated volume of each material in all mesh elements. Three sets of
         rays are used: one set parallel to the x-axis, one parallel to the
-        y-axis, and one parallel to the z-axis.
+        y-axis, and one parallel to the z-axis. Regions of the mesh that are
+        outside the model geometry are treated as void, equivalent to a cell
+        with no material. Universe fills within the model must still define all
+        enclosed space.
 
         .. versionadded:: 0.15.1
 

src/dagmc.cpp

@@ -857,12 +857,15 @@ std::pair<double, int32_t> DAGCell::distance(
       dag_univ->surf_idx_offset_ + dagmc_ptr_->index_by_handle(hit_surf);
   } else if (!dagmc_ptr_->is_implicit_complement(vol) ||
              is_root_universe(dag_univ->id_)) {
-    // surface boundary conditions are ignored for projection plotting, meaning
+    // Surface boundary conditions are ignored for projection plotting, meaning
     // that the particle may move through the graveyard (bounding) volume and
-    // into the implicit complement on the other side where no intersection will
-    // be found. Treating this as a lost particle is problematic when plotting.
-    // Instead, the infinite distance and invalid surface index are returned.
-    if (settings::run_mode == RunMode::PLOTTING)
+    // into the implicit complement on the other side where no intersection
+    // will be found. A no-hit result is also expected when querying root cells
+    // for the next boundary from undefined space, when no containing cell is
+    // assigned. In both cases, return an infinite distance and invalid surface
+    // index rather than marking a particle as lost.
+    if (settings::run_mode == RunMode::PLOTTING ||
+        p->lowest_coord().cell() == C_NONE)
       return {INFTY, -1};
 
     // the particle should be marked as lost immediately if an intersection

src/geometry.cpp

@@ -281,6 +281,16 @@ bool find_cell_inner(
 bool neighbor_list_find_cell(GeometryState& p, bool verbose)
 {
 
+#ifdef OPENMC_DAGMC_ENABLED
+  // A CSG crossing can move the particle into another instance of the same
+  // DAGMC universe, where the previous facet history is no longer valid.
+  if (p.surface() != SURFACE_NONE) {
+    const auto& surf = model::surfaces[p.surface_index()];
+    if (surf->geom_type() == GeometryType::CSG)
+      p.history().reset();
+  }
+#endif
+
   // Reset all the deeper coordinate levels.
   for (int i = p.n_coord(); i < model::n_coord_levels; i++) {
     p.coord(i).reset();

src/mesh.cpp

@@ -482,9 +482,6 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
   width.y = (ny > 0) ? width.y / ny : 0.0;
   width.z = (nz > 0) ? width.z / nz : 0.0;
 
-  // Set flag for mesh being contained within model
-  bool out_of_model = false;
-
 #pragma omp parallel
   {
     // Preallocate vector for mesh indices and length fractions and particle
@@ -496,6 +493,32 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
     site.E = 1.0;
     site.particle = ParticleType::neutron();
 
+    bool verbose = settings::verbosity >= 10;
+
+    // Save the cells occupied immediately before a boundary crossing.
+    auto save_cell_state = [&p]() {
+      for (int j = 0; j < p.n_coord(); ++j) {
+        p.cell_last(j) = p.coord(j).cell();
+      }
+      p.n_coord_last() = p.n_coord();
+    };
+
+    // Initialize cell history after locating a ray inside the model.
+    auto initialize_cell_state = [&p, &save_cell_state]() {
+      if (p.cell_born() == C_NONE)
+        p.cell_born() = p.lowest_coord().cell();
+
+      save_cell_state();
+    };
+
+    // Reset a failed coordinate search while preserving position and direction.
+    auto reset_geometry_state = [&p]() {
+      Position r = p.r();
+      Direction u = p.u();
+      p.init_from_r_u(r, u);
+      p.coord(0).universe() = model::root_universe;
+    };
+
     for (int axis = 0; axis < 3; ++axis) {
       // Set starting position and direction
       site.r = {0.0, 0.0, 0.0};
@@ -524,6 +547,50 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
       int i1_start = mpi::rank * min_work + std::min(mpi::rank, remainder);
       int i1_end = i1_start + n1_local;
 
+      // Add the contribution from a ray segment. The positions used here are
+      // kept separate from the particle position because the latter is moved a
+      // tiny distance across each surface for robust geometry searches.
+      auto add_segment = [&](const Position& r0, const Position& r1,
+                           int i_material) {
+        double distance = r1[axis] - r0[axis];
+        if (distance <= 0.0)
+          return;
+
+        bins.clear();
+        length_fractions.clear();
+        this->bins_crossed(r0, r1, site.u, bins, length_fractions);
+
+        double cumulative_frac = 0.0;
+        for (int i_bin = 0; i_bin < bins.size(); i_bin++) {
+          int mesh_index = bins[i_bin];
+          double length = distance * length_fractions[i_bin];
+          double volume = length * d1 * d2;
+
+          if (compute_bboxes) {
+            double axis_start = r0[axis] + distance * cumulative_frac;
+            double axis_end = axis_start + length;
+            cumulative_frac += length_fractions[i_bin];
+
+            Position contrib_min = site.r;
+            Position contrib_max = site.r;
+
+            contrib_min[ax1] = site.r[ax1] - 0.5 * d1;
+            contrib_max[ax1] = site.r[ax1] + 0.5 * d1;
+            contrib_min[ax2] = site.r[ax2] - 0.5 * d2;
+            contrib_max[ax2] = site.r[ax2] + 0.5 * d2;
+            contrib_min[axis] = std::min(axis_start, axis_end);
+            contrib_max[axis] = std::max(axis_start, axis_end);
+
+            BoundingBox contrib_bbox {contrib_min, contrib_max};
+            contrib_bbox &= bbox;
+
+            result.add_volume(mesh_index, i_material, volume, &contrib_bbox);
+          } else {
+            result.add_volume(mesh_index, i_material, volume);
+          }
+        }
+      };
+
       // Loop over rays on face of bounding box
 #pragma omp for collapse(2)
       for (int i1 = i1_start; i1 < i1_end; ++i1) {
@@ -533,98 +600,115 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
 
           p.from_source(&site);
 
+          // Set the physical endpoint of this ray at the far mesh face.
+          Position r_mesh_end = site.r;
+          r_mesh_end[axis] = bbox.max[axis];
+
           // Determine particle's location
-          if (!exhaustive_find_cell(p)) {
-            out_of_model = true;
-            continue;
+          bool inside_model = exhaustive_find_cell(p, verbose);
+
+          if (inside_model) {
+            initialize_cell_state();
+          } else {
+            // Clear any partial descent into nested universes before searching
+            // for the first root-universe boundary from undefined space.
+            reset_geometry_state();
           }
 
-          // Set birth cell attribute
-          if (p.cell_born() == C_NONE)
-            p.cell_born() = p.lowest_coord().cell();
+          // Physical position through which volume has been accumulated. This
+          // differs by TINY_BIT from p.r() after crossing a surface.
+          Position r_scored = site.r;
+
+          while (r_scored[axis] < r_mesh_end[axis]) {
+            if (!inside_model) {
+              // The ray is outside the model. Advance to the next surface of
+              // any cell in the root universe, as is done for ray-traced
+              // plots. Undefined space traversed along the way is void.
+              Position r0 = p.r();
+              p.advance_to_boundary_from_void();
+
+              // If no model surface lies before the mesh edge, score the
+              // remaining exterior interval as void and finish the ray.
+              double distance_to_mesh_end = r_mesh_end[axis] - r0[axis];
+              if (p.boundary().surface() == SURFACE_NONE ||
+                  p.boundary().distance() >= distance_to_mesh_end) {
+                add_segment(r_scored, r_mesh_end, MATERIAL_VOID);
+                break;
+              }
 
-          // Initialize last cells from current cell
-          for (int j = 0; j < p.n_coord(); ++j) {
-            p.cell_last(j) = p.coord(j).cell();
-          }
-          p.n_coord_last() = p.n_coord();
+              // Determine the physical position of the model boundary.
+              Position r_boundary = r0 + p.boundary().distance() * p.u();
 
-          while (true) {
-            // Ray trace from r_start to r_end
-            Position r0 = p.r();
-            double max_distance = bbox.max[axis] - r0[axis];
+              // Score the exterior interval and record its physical endpoint.
+              add_segment(r_scored, r_boundary, MATERIAL_VOID);
+              r_scored = r_boundary;
+
+              // Check whether advancing through the surface entered the model.
+              inside_model = exhaustive_find_cell(p, verbose);
+              if (inside_model) {
+                initialize_cell_state();
+              } else {
+                // Clear any partial coordinate search before looking for the
+                // next surface from undefined space.
+                reset_geometry_state();
+              }
+              continue;
+            }
 
             // Find the distance to the nearest boundary
             BoundaryInfo boundary = distance_to_boundary(p);
 
-            // Advance particle forward
-            double distance = std::min(boundary.distance(), max_distance);
-            p.move_distance(distance);
-
-            // Determine what mesh elements were crossed by particle
-            bins.clear();
-            length_fractions.clear();
-            this->bins_crossed(r0, p.r(), p.u(), bins, length_fractions);
-
-            // Add volumes to any mesh elements that were crossed
+            // Convert the material index to a user-facing ID
             int i_material = p.material();
             if (i_material != C_NONE) {
               i_material = model::materials[i_material]->id();
             }
-            double cumulative_frac = 0.0;
-            for (int i_bin = 0; i_bin < bins.size(); i_bin++) {
-              int mesh_index = bins[i_bin];
-              double length = distance * length_fractions[i_bin];
-              double volume = length * d1 * d2;
 
-              if (compute_bboxes) {
-                double axis_start = r0[axis] + distance * cumulative_frac;
-                double axis_end = axis_start + length;
-                cumulative_frac += length_fractions[i_bin];
+            // If no model boundary lies before the mesh edge, score the
+            // remaining material interval and finish the ray.
+            double distance_to_mesh_end = r_mesh_end[axis] - p.r()[axis];
+            if (boundary.distance() >= distance_to_mesh_end) {
+              add_segment(r_scored, r_mesh_end, i_material);
+              break;
+            }
 
-                Position contrib_min = site.r;
-                Position contrib_max = site.r;
+            // Determine the physical position of the model boundary.
+            Position r_boundary = p.r() + boundary.distance() * p.u();
 
-                contrib_min[ax1] = site.r[ax1] - 0.5 * d1;
-                contrib_max[ax1] = site.r[ax1] + 0.5 * d1;
-                contrib_min[ax2] = site.r[ax2] - 0.5 * d2;
-                contrib_max[ax2] = site.r[ax2] + 0.5 * d2;
-                contrib_min[axis] = std::min(axis_start, axis_end);
-                contrib_max[axis] = std::max(axis_start, axis_end);
+            // Score the material interval and record its physical endpoint.
+            add_segment(r_scored, r_boundary, i_material);
+            r_scored = r_boundary;
 
-                BoundingBox contrib_bbox {contrib_min, contrib_max};
-                contrib_bbox &= bbox;
+            // Cross the next geometric surface. The small forward movement
+            // and neighbor-list search mirror Ray::trace, allowing a failed
+            // search to mean that the ray has left the model rather than that
+            // a transport particle has been lost.
+            save_cell_state();
 
-                result.add_volume(
-                  mesh_index, i_material, volume, &contrib_bbox);
-              } else {
-                // Add volume to result
-                result.add_volume(mesh_index, i_material, volume);
-              }
-            }
-
-            if (distance == max_distance)
-              break;
-
-            // cross next geometric surface
-            for (int j = 0; j < p.n_coord(); ++j) {
-              p.cell_last(j) = p.coord(j).cell();
-            }
-            p.n_coord_last() = p.n_coord();
+            // Move just beyond the surface to make the next search robust.
+            p.move_distance(boundary.distance() + TINY_BIT);
 
             // Set surface that particle is on and adjust coordinate levels
             p.surface() = boundary.surface();
             p.n_coord() = boundary.coord_level();
 
+            // Update the geometry state according to the boundary type.
             if (boundary.lattice_translation()[0] != 0 ||
                 boundary.lattice_translation()[1] != 0 ||
                 boundary.lattice_translation()[2] != 0) {
               // Particle crosses lattice boundary
-              cross_lattice(p, boundary);
+              cross_lattice(p, boundary, verbose);
+              inside_model = true;
             } else {
-              // Particle crosses surface
-              const auto& surf {model::surfaces[p.surface_index()].get()};
-              p.cross_surface(*surf);
+              // Search for the cell on the opposite side of a surface.
+              inside_model = neighbor_list_find_cell(p, verbose);
+            }
+
+            // Treat a failed cell search as a transition to exterior void.
+            if (!inside_model) {
+              // Reset the geometry state so the next iteration can search for
+              // another disjoint portion of the model.
+              reset_geometry_state();
             }
           }
         }
@@ -633,9 +717,7 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
   }
 
   // Check for errors
-  if (out_of_model) {
-    throw std::runtime_error("Mesh not fully contained in geometry.");
-  } else if (result.table_full()) {
+  if (result.table_full()) {
     throw std::runtime_error("Maximum number of materials for mesh material "
                              "volume calculation insufficient.");
   }

src/particle.cpp

@@ -663,12 +663,6 @@ void Particle::cross_surface(const Surface& surf)
     write_message(1, "    Crossing surface {}", surf.id_);
   }
 
-// if we're crossing a CSG surface, make sure the DAG history is reset
-#ifdef OPENMC_DAGMC_ENABLED
-  if (surf.geom_type() == GeometryType::CSG)
-    history().reset();
-#endif
-
   // Handle any applicable boundary conditions.
   if (surf.bc_ && settings::run_mode != RunMode::PLOTTING &&
       settings::run_mode != RunMode::VOLUME) {