Merge pull request #1193 from paulromano/various-patches

Hotfix for heating data, and other patches

CMakeLists.txt

@@ -255,13 +255,12 @@ set_target_properties(libopenmc PROPERTIES
   OUTPUT_NAME openmc)
 
 target_include_directories(libopenmc
-  PUBLIC include
-  PRIVATE ${HDF5_INCLUDE_DIRS})
+  PUBLIC include ${HDF5_INCLUDE_DIRS})
 
 # Set compile flags
 target_compile_options(libopenmc PRIVATE ${cxxflags})
 
-target_compile_definitions(libopenmc PRIVATE -DMAX_COORD=${maxcoord})
+target_compile_definitions(libopenmc PUBLIC -DMAX_COORD=${maxcoord})
 if (HDF5_IS_PARALLEL)
   target_compile_definitions(libopenmc PRIVATE -DPHDF5)
 endif()

docs/source/publications.rst

@@ -210,6 +210,11 @@ Multigroup Cross Section Generation
   transport simulations <https://doi.org/10.1016/j.anucene.2018.11.017>`_,"
   *Ann. Nucl. Energy*, **125**, 261-271 (2019).
 
+- Kun Zhuang, Xiaobin Tang, and Liangzhi Cao, "`Development and verification of
+  a model for generation of MSFR few-group homogenized cross-sections based on a
+  Monte Carlo code OpenMC <https://doi.org/10.1016/j.anucene.2018.09.037>`_,"
+  *Ann. Nucl. Energy*, **124**, 187-197 (2019).
+
 - Changho Lee and Yeon Sang Jung, "Verification of the Cross Section Library
   Generated Using OpenMC and MC\ :sup:`2`-3 for PROTEUS," *Proc. PHYSOR*, Cancun,
   Mexico, Apr. 22-26 (2018).

include/openmc/lattice.h

@@ -152,8 +152,7 @@ public:
   int indx_;  //!< An index to a Lattice universes or offsets array.
 
   LatticeIter(Lattice &lat, int indx)
-    : lat_(lat),
-      indx_(indx)
+    : indx_(indx), lat_(lat)
   {}
 
   bool operator==(const LatticeIter &rhs) {return (indx_ == rhs.indx_);}

openmc/capi/core.py

@@ -17,7 +17,8 @@ class _Bank(Structure):
                 ('u', c_double*3),
                 ('E', c_double),
                 ('wgt', c_double),
-                ('delayed_group', c_int)]
+                ('delayed_group', c_int),
+                ('particle', c_int)]
 
 
 # Define input type for numpy arrays that will be passed into C++ functions

openmc/data/neutron.py

@@ -622,7 +622,7 @@ class IncidentNeutron(EqualityMixin):
         data.energy[strT] = energy
         total_xs = ace.xss[i + n_energy : i + 2*n_energy]
         absorption_xs = ace.xss[i + 2*n_energy : i + 3*n_energy]
-        heating_number = ace.xss[i + 3*n_energy : i + 4*n_energy]*EV_PER_MEV
+        heating_number = ace.xss[i + 4*n_energy : i + 5*n_energy]*EV_PER_MEV
 
         # Create redundant reactions (total, absorption, and heating)
         total = Reaction(1)
@@ -637,7 +637,7 @@ class IncidentNeutron(EqualityMixin):
             data.reactions[101] = absorption
 
         heating = Reaction(301)
-        heating.xs[strT] = Tabulated1D(energy, heating_number)
+        heating.xs[strT] = Tabulated1D(energy, heating_number*total_xs)
         heating.redundant = True
         data.reactions[301] = heating
 

openmc/geometry.py

@@ -97,9 +97,14 @@ class Geometry(object):
         # Clean the indentation in the file to be user-readable
         xml.clean_indentation(root_element)
 
+        # Check if path is a directory
+        p = Path(path)
+        if p.is_dir():
+            p /= 'geometry.xml'
+
         # Write the XML Tree to the geometry.xml file
         tree = ET.ElementTree(root_element)
-        tree.write(path, xml_declaration=True, encoding='utf-8')
+        tree.write(str(p), xml_declaration=True, encoding='utf-8')
 
     @classmethod
     def from_xml(cls, path='geometry.xml', materials=None):

openmc/material.py

@@ -1,6 +1,7 @@
 from collections import OrderedDict
 from copy import deepcopy
 from numbers import Real, Integral
+from pathlib import Path
 import warnings
 from xml.etree import ElementTree as ET
 
@@ -1065,9 +1066,14 @@ class Materials(cv.CheckedList):
         # Clean the indentation in the file to be user-readable
         clean_indentation(root_element)
 
+        # Check if path is a directory
+        p = Path(path)
+        if p.is_dir():
+            p /= 'materials.xml'
+
         # Write the XML Tree to the materials.xml file
         tree = ET.ElementTree(root_element)
-        tree.write(path, xml_declaration=True, encoding='utf-8')
+        tree.write(str(p), xml_declaration=True, encoding='utf-8')
 
     @classmethod
     def from_xml(cls, path='materials.xml'):

openmc/model/model.py

@@ -1,4 +1,5 @@
 from collections.abc import Iterable
+from pathlib import Path
 
 import openmc
 from openmc.checkvalue import check_type, check_value
@@ -154,27 +155,39 @@ class Model(object):
                                        'si_celi', 'si_leqi', 'celi', 'leqi'))
         getattr(dep.integrator, method)(op, timesteps, **kwargs)
 
-    def export_to_xml(self):
-        """Export model to XML files."""
+    def export_to_xml(self, directory='.'):
+        """Export model to XML files.
 
-        self.settings.export_to_xml()
+        Parameters
+        ----------
+        directory : str
+            Directory to write XML files to. If it doesn't exist already, it
+            will be created.
+
+        """
+        # Create directory if
+        d = Path(directory)
+        if not d.is_dir():
+            d.mkdir(parents=True)
+
+        self.settings.export_to_xml(d)
         if not self.settings.dagmc:
-            self.geometry.export_to_xml()
+            self.geometry.export_to_xml(d)
 
         # If a materials collection was specified, export it. Otherwise, look
         # for all materials in the geometry and use that to automatically build
         # a collection.
         if self.materials:
-            self.materials.export_to_xml()
+            self.materials.export_to_xml(d)
         else:
             materials = openmc.Materials(self.geometry.get_all_materials()
                                          .values())
-            materials.export_to_xml()
+            materials.export_to_xml(d)
 
         if self.tallies:
-            self.tallies.export_to_xml()
+            self.tallies.export_to_xml(d)
         if self.plots:
-            self.plots.export_to_xml()
+            self.plots.export_to_xml(d)
 
     def run(self, **kwargs):
         """Creates the XML files, runs OpenMC, and returns k-effective

openmc/plots.py

@@ -1,9 +1,10 @@
 from collections.abc import Iterable, Mapping
 from numbers import Real, Integral
-from xml.etree import ElementTree as ET
+from pathlib import Path
 import subprocess
 import sys
 import warnings
+from xml.etree import ElementTree as ET
 
 import numpy as np
 
@@ -818,6 +819,11 @@ class Plots(cv.CheckedList):
         # Clean the indentation in the file to be user-readable
         clean_indentation(self._plots_file)
 
+        # Check if path is a directory
+        p = Path(path)
+        if p.is_dir():
+            p /= 'plots.xml'
+
         # Write the XML Tree to the plots.xml file
         tree = ET.ElementTree(self._plots_file)
-        tree.write(path, xml_declaration=True, encoding='utf-8', method="xml")
+        tree.write(str(p), xml_declaration=True, encoding='utf-8')

openmc/settings.py

@@ -1,4 +1,5 @@
 from collections.abc import Iterable, MutableSequence, Mapping
+from pathlib import Path
 from numbers import Real, Integral
 import warnings
 from xml.etree import ElementTree as ET
@@ -995,6 +996,11 @@ class Settings(object):
         # Clean the indentation in the file to be user-readable
         clean_indentation(root_element)
 
+        # Check if path is a directory
+        p = Path(path)
+        if p.is_dir():
+            p /= 'settings.xml'
+
         # Write the XML Tree to the settings.xml file
         tree = ET.ElementTree(root_element)
-        tree.write(path, xml_declaration=True, encoding='utf-8', method="xml")
+        tree.write(str(p), xml_declaration=True, encoding='utf-8')

openmc/tallies.py

@@ -5,6 +5,7 @@ from functools import partial, reduce
 from itertools import product
 from numbers import Integral, Real
 import operator
+from pathlib import Path
 import warnings
 from xml.etree import ElementTree as ET
 
@@ -3189,7 +3190,11 @@ class Tallies(cv.CheckedList):
         # Clean the indentation in the file to be user-readable
         clean_indentation(root_element)
 
+        # Check if path is a directory
+        p = Path(path)
+        if p.is_dir():
+            p /= 'tallies.xml'
+
         # Write the XML Tree to the tallies.xml file
         tree = ET.ElementTree(root_element)
-        tree.write(path, xml_declaration=True,
-                   encoding='utf-8', method="xml")
+        tree.write(str(p), xml_declaration=True, encoding='utf-8')

src/physics.cpp

@@ -267,7 +267,7 @@ void sample_photon_reaction(Particle* p)
     double phi = 2.0*PI*prn();
     double E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b;
     int electron = static_cast<int>(Particle::Type::electron);
-    if (E_electron >= settings::energy_cutoff[electron]) { 
+    if (E_electron >= settings::energy_cutoff[electron]) {
       double mu_electron = (alpha - alpha_out*mu)
         / std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu);
       Direction u = rotate_angle(p->u(), mu_electron, &phi);
@@ -504,6 +504,9 @@ Reaction* sample_fission(int i_nuclide, double E)
     // Create fission bank sites if fission occurs
     if (prob > cutoff) return rx;
   }
+
+  // If we reached here, no reaction was sampled
+  throw std::runtime_error{"No fission reaction was sampled for " + nuc->name_};
 }
 
 void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product)
@@ -901,8 +904,11 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
         }
       }
     }
-  } // case RVS, DBRC
+    } // case RVS, DBRC
   } // switch (sampling_method)
+
+  throw std::runtime_error{"Unable to sample target velocity for "
+    "elastic scattering."};
 }
 
 Direction

src/state_point.cpp

@@ -480,6 +480,7 @@ hid_t h5banktype() {
   H5Tinsert(banktype, "E", HOFFSET(Particle::Bank, E), H5T_NATIVE_DOUBLE);
   H5Tinsert(banktype, "wgt", HOFFSET(Particle::Bank, wgt), H5T_NATIVE_DOUBLE);
   H5Tinsert(banktype, "delayed_group", HOFFSET(Particle::Bank, delayed_group), H5T_NATIVE_INT);
+  H5Tinsert(banktype, "particle", HOFFSET(Particle::Bank, particle), H5T_NATIVE_INT);
 
   H5Tclose(postype);
   return banktype;

src/tallies/filter_sph_harm.cpp

@@ -80,9 +80,10 @@ SphericalHarmonicsFilter::text_label(int bin) const
     if (bin < (n + 1) * (n + 1)) {
       int m = (bin - n*n) - n;
       out << "Spherical harmonic expansion, Y" << n << "," << m;
-      return out.str();
+      break;
     }
   }
+  return out.str();
 }
 
 //==============================================================================

src/tallies/filter_zernike.cpp

@@ -65,9 +65,10 @@ ZernikeFilter::text_label(int bin) const
       int first = last - (n + 1);
       int m = -n + (bin - first) * 2;
       out << "Zernike expansion, Z" << n << "," << m;
-      return out.str();
+      break;
     }
   }
+  return out.str();
 }
 
 void

src/tallies/trigger.cpp

@@ -83,6 +83,9 @@ check_tally_triggers(double& ratio, int& tally_id, int& score)
             case TriggerMetric::relative_error:
               uncertainty = rel_err;
               break;
+            case TriggerMetric::not_active:
+              uncertainty = 0.0;
+              break;
           }
 
           // Compute the uncertainty / threshold ratio.
@@ -109,7 +112,6 @@ double
 check_keff_trigger()
 {
   if (settings::run_mode != RUN_MODE_EIGENVALUE) return 0.;
-  if (settings::keff_trigger.metric == TriggerMetric::not_active) return 0.;
 
   double k_combined[2];
   int err = openmc_get_keff(k_combined);
@@ -124,6 +126,9 @@ check_keff_trigger()
     break;
   case TriggerMetric::relative_error:
     uncertainty = k_combined[1] / k_combined[0];
+    break;
+  case TriggerMetric::not_active:
+    break;
   }
 
   double ratio = uncertainty / settings::keff_trigger.threshold;

tests/regression_tests/tallies/results_true.dat

@@ -1 +1 @@
-de773a799f84348241ebd65bf7beae8114861d8674b652bfd443d578c146a7d37ca38f2efc5d05124fdbb1cf12829907768351e6b2d6b5f4bd2c121417757507
+c37e0468f1684f7810429332721762884f3cdf0429d38caa99500bbde04ad84cf1e1639adc46b6db44b28b6f31028ed029f85920c8164d77463b6d89336043a8