Merge branch 'openmc-dev:develop' into Sensitivity-Analysis

Author: Ebiwonjumi

include/openmc/constants.h

@@ -63,6 +63,11 @@ constexpr int MAX_SAMPLE {100000};
 // source region in the random ray solver
 constexpr double MIN_HITS_PER_BATCH {1.5};
 
+// The minimum flux value to be considered non-zero when computing adjoint
+// sources. Positive values below this cutoff will be treated as zero, so as to
+// prevent extremely large adjoint source terms from being generated.
+constexpr double ZERO_FLUX_CUTOFF {1e-22};
+
 // ============================================================================
 // MATH AND PHYSICAL CONSTANTS
 

openmc/dagmc.py

@@ -13,6 +13,7 @@
 from .surface import _BOUNDARY_TYPES
 from .bounding_box import BoundingBox
 from .utility_funcs import input_path
+from .plots import add_plot_params
 
 
 class DAGMCUniverse(openmc.UniverseBase):
@@ -566,6 +567,12 @@ def sync_dagmc_cells(self, mats: Iterable[openmc.Material]):
                 fill = mats_per_id[dag_cell.fill.id] if dag_cell.fill else None
             self.add_cell(openmc.DAGMCCell(cell_id=dag_cell_id, fill=fill))
 
+    @add_plot_params
+    def plot(self, *args, **kwargs):
+        """Display a slice plot of the DAGMCUniverse.
+        """
+        return openmc.Geometry(self).plot(*args, **kwargs)
+
 
 class DAGMCCell(openmc.Cell):
     """A cell class for DAGMC-based geometries.

openmc/data/data.py

@@ -94,7 +94,7 @@
     'Yb174': 0.32025, 'Yb176': 0.12995, 'Lu175': 0.97401,
     'Lu176': 0.02599, 'Hf174': 0.0016, 'Hf176': 0.0526,
     'Hf177': 0.186, 'Hf178': 0.2728, 'Hf179': 0.1362,
-    'Hf180': 0.3508, 'Ta180': 0.0001201, 'Ta181': 0.9998799,
+    'Hf180': 0.3508, 'Ta180_m1': 0.0001201, 'Ta181': 0.9998799,
     'W180': 0.0012, 'W182': 0.265, 'W183': 0.1431,
     'W184': 0.3064, 'W186': 0.2843, 'Re185': 0.374,
     'Re187': 0.626, 'Os184': 0.0002, 'Os186': 0.0159,

openmc/element.py

@@ -144,8 +144,7 @@ def expand(self, percent, percent_type, enrichment=None,
             root = tree.getroot()
             for child in root.findall('library'):
                 nuclide = child.attrib['materials']
-                if re.match(r'{}\d+'.format(self), nuclide) and \
-                   '_m' not in nuclide:
+                if re.match(r'{}\d+'.format(self), nuclide):
                     library_nuclides.add(nuclide)
 
             # Get a set of the mutual and absent nuclides. Convert to lists
@@ -185,8 +184,10 @@ def expand(self, percent, percent_type, enrichment=None,
                 for nuclide in absent_nuclides:
                     if nuclide in ['O17', 'O18'] and 'O16' in mutual_nuclides:
                         abundances['O16'] += NATURAL_ABUNDANCE[nuclide]
-                    elif nuclide == 'Ta180' and 'Ta181' in mutual_nuclides:
-                        abundances['Ta181'] += NATURAL_ABUNDANCE[nuclide]
+                    elif nuclide == 'Ta180_m1' and 'Ta180' in library_nuclides:
+                            abundances['Ta180'] = NATURAL_ABUNDANCE[nuclide]
+                    elif nuclide == 'Ta180_m1' and 'Ta181' in mutual_nuclides:
+                            abundances['Ta181'] += NATURAL_ABUNDANCE[nuclide]
                     elif nuclide == 'W180' and 'W182' in mutual_nuclides:
                         abundances['W182'] += NATURAL_ABUNDANCE[nuclide]
                     else:

openmc/geometry.py

@@ -754,4 +754,17 @@ def plot(self, *args, **kwargs):
 
         .. versionadded:: 0.14.0
         """
-        return self.root_universe.plot(*args, **kwargs)
+        model = openmc.Model()
+        model.geometry = self
+        model.materials = self.get_all_materials().values()
+
+        # Add placeholder materials for DAGMCUniverses
+        universes = self.get_all_universes()
+        for universe in universes.values():
+            if isinstance(universe, openmc.DAGMCUniverse):
+                for name in universe.material_names:
+                    mat_dag = openmc.Material(name=name)
+                    mat_dag.add_nuclide('H1', 1.0)
+                    model.materials.append(mat_dag)
+
+        return model.plot(*args, **kwargs)

openmc/model/model.py

@@ -986,7 +986,16 @@ def plot(
         y_min = (origin[y] - 0.5*width[1]) * axis_scaling_factor[axis_units]
         y_max = (origin[y] + 0.5*width[1]) * axis_scaling_factor[axis_units]
 
+        # Determine whether any materials contains macroscopic data and if so,
+        # set energy mode accordingly
+        _energy_mode = self.settings._energy_mode
+        for mat in self.geometry.get_all_materials().values():
+            if mat._macroscopic is not None:
+                self.settings.energy_mode = 'multi-group'
+                break
+
         with TemporaryDirectory() as tmpdir:
+            _plot_seed = self.settings.plot_seed
             if seed is not None:
                 self.settings.plot_seed = seed
 
@@ -1007,6 +1016,11 @@ def plot(
             # Run OpenMC in geometry plotting mode
             self.plot_geometry(False, cwd=tmpdir, openmc_exec=openmc_exec)
 
+            # Undo changes to model
+            self.plots.pop()
+            self.settings._plot_seed = _plot_seed
+            self.settings._energy_mode = _energy_mode
+
             # Read image from file
             img_path = Path(tmpdir) / f'plot_{plot.id}.png'
             if not img_path.is_file():

openmc/universe.py

@@ -337,14 +337,6 @@ def plot(self, *args, **kwargs):
         """
         model = openmc.Model()
         model.geometry = openmc.Geometry(self)
-
-        # Determine whether any materials contains macroscopic data and if
-        # so, set energy mode accordingly
-        for mat in self.get_all_materials().values():
-            if mat._macroscopic is not None:
-                model.settings.energy_mode = 'multi-group'
-                break
-
         return model.plot(*args, **kwargs)
 
     def get_nuclides(self):

src/random_ray/flat_source_domain.cpp

@@ -304,6 +304,13 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux()
           set_flux_to_source(sr, g);
         }
       }
+      // Halt if NaN implosion is detected
+      if (!std::isfinite(source_regions_.scalar_flux_new(sr, g))) {
+        fatal_error("A source region scalar flux is not finite. "
+                    "This indicates a numerical instability in the "
+                    "simulation. Consider increasing ray density or adjusting "
+                    "the source region mesh.");
+      }
     }
   }
 
@@ -1154,16 +1161,21 @@ void FlatSourceDomain::flatten_xs()
           m.get_xs(MgxsType::TOTAL, g_out, NULL, NULL, NULL, t, a);
         sigma_t_.push_back(sigma_t);
 
-        double nu_Sigma_f =
+        double nu_sigma_f =
           m.get_xs(MgxsType::NU_FISSION, g_out, NULL, NULL, NULL, t, a);
-        nu_sigma_f_.push_back(nu_Sigma_f);
+        nu_sigma_f_.push_back(nu_sigma_f);
 
         double sigma_f =
           m.get_xs(MgxsType::FISSION, g_out, NULL, NULL, NULL, t, a);
         sigma_f_.push_back(sigma_f);
 
         double chi =
           m.get_xs(MgxsType::CHI_PROMPT, g_out, &g_out, NULL, NULL, t, a);
+        if (!std::isfinite(chi)) {
+          // MGXS interface may return NaN in some cases, such as when material
+          // is fissionable but has very small sigma_f.
+          chi = 0.0;
+        }
         chi_.push_back(chi);
 
         for (int g_in = 0; g_in < negroups_; g_in++) {
@@ -1191,17 +1203,30 @@ void FlatSourceDomain::flatten_xs()
 
 void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
 {
-  // Set the external source to 1/forward_flux. If the forward flux is negative
-  // or zero, set the adjoint source to zero, as this is likely a very small
-  // source region that we don't need to bother trying to vector particles
-  // towards. Flux negativity in random ray is not related to the flux being
-  // small in magnitude, but rather due to the source region being physically
-  // small in volume and thus having a noisy flux estimate.
+  // Set the adjoint external source to 1/forward_flux. If the forward flux is
+  // negative, zero, or extremely close to zero, set the adjoint source to zero,
+  // as this is likely a very small source region that we don't need to bother
+  // trying to vector particles towards. In the case of flux "being extremely
+  // close to zero", we define this as being a fixed fraction of the maximum
+  // forward flux, below which we assume the flux would be physically
+  // undetectable.
+
+  // First, find the maximum forward flux value
+  double max_flux = 0.0;
+#pragma omp parallel for reduction(max : max_flux)
+  for (int64_t se = 0; se < n_source_elements(); se++) {
+    double flux = forward_flux[se];
+    if (flux > max_flux) {
+      max_flux = flux;
+    }
+  }
+
+  // Then, compute the adjoint source for each source region
 #pragma omp parallel for
   for (int64_t sr = 0; sr < n_source_regions(); sr++) {
     for (int g = 0; g < negroups_; g++) {
       double flux = forward_flux[sr * negroups_ + g];
-      if (flux <= 0.0) {
+      if (flux <= ZERO_FLUX_CUTOFF * max_flux) {
         source_regions_.external_source(sr, g) = 0.0;
       } else {
         source_regions_.external_source(sr, g) = 1.0 / flux;

tests/unit_tests/dagmc/test_plot.py

@@ -7,9 +7,9 @@
     not openmc.lib._dagmc_enabled(), reason="DAGMC CAD geometry is not enabled."
 )
 
-def test_plotting_dagmc_geometry(request):
-    """Test plotting a DAGMC geometry with OpenMC. This is different to CSG
-    geometry plotting as the path to the DAGMC file needs handling."""
+def test_plotting_dagmc_model(request):
+    """Test plotting a DAGMC model with OpenMC. This is different to CSG
+    model plotting as the path to the DAGMC file needs handling."""
 
     dag_universe = openmc.DAGMCUniverse(request.path.parent / 'dagmc.h5m')
     csg_with_dag_inside = dag_universe.bounded_universe()
@@ -30,3 +30,32 @@ def test_plotting_dagmc_geometry(request):
     model.settings.particles = 50
 
     model.plot()
+
+
+def test_plotting_dagmc_universe(request):
+    """Test plotting a DAGMCUniverse with OpenMC. This is different to plotting
+    UniverseBase as the materials are not defined withing the DAGMCUniverse."""
+
+    dag_universe = openmc.DAGMCUniverse(request.path.parent / 'dagmc.h5m')
+    dag_universe.plot()
+
+
+def test_plotting_geometry_filled_with_dagmc_universe(request):
+    """Test plotting a geometry with OpenMC. This is an edge case when plotting
+    geometry as often geometry objects don't include a DAGMCUniverse. The
+    inclusion of a DAGMCUniverse requires special handling for the materials."""
+
+    dag_universe = openmc.DAGMCUniverse(request.path.parent / 'dagmc.h5m', auto_geom_ids=True)
+
+    sphere1 = openmc.Sphere(r=50.0)
+    sphere2 = openmc.Sphere(r=60.0, boundary_type='vacuum')
+
+    # Adding a material to the CSG Universe to check all materials are accounted for
+    csg_material = openmc.Material(name='csg_material')
+    csg_material.add_nuclide("H1", 1.0)
+
+    cell1 = openmc.Cell(fill=dag_universe, region=-sphere1)
+    cell2 = openmc.Cell(fill=csg_material, region=+sphere1 & -sphere2)
+
+    geometry = openmc.Geometry([cell1, cell2])
+    geometry.plot()

tests/unit_tests/test_element.py

@@ -44,6 +44,13 @@ def test_expand_no_isotopes():
         element.expand(100.0, 'ao')
 
 
+def test_expand_ta():
+    ref = {'Ta180': 0.01201, 'Ta181': 99.98799}
+    element = openmc.Element('Ta')
+    for isotope in element.expand(100.0, 'ao'):
+        assert isotope[1] == approx(ref[isotope[0]])
+    
+
 def test_expand_exceptions():
     """ Test that correct exceptions are raised for invalid input """