Merge branch 'mcpl_input' of github.com:ebknudsen/openmc into mcpl_input

include/openmc/settings.h

@@ -47,6 +47,7 @@ extern "C" bool run_CE;            //!< run with continuous-energy data?
 extern bool source_latest;         //!< write latest source at each batch?
 extern bool source_separate;       //!< write source to separate file?
 extern bool source_write;          //!< write source in HDF5 files?
+extern bool source_mcpl_write;     //!< write source in mcpl files?
 extern bool surf_source_write;     //!< write surface source file?
 extern bool surf_mcpl_write;       //!< write surface mcpl file?
 extern bool surf_source_read;      //!< read surface source file?

openmc/data/decay.py

@@ -1,10 +1,9 @@
 from collections.abc import Iterable
 from io import StringIO
 from math import log
-from pathlib import Path
 import re
+from typing import Optional
 from warnings import warn
-from xml.etree import ElementTree as ET
 
 import numpy as np
 from uncertainties import ufloat, UFloat
@@ -579,7 +578,7 @@ class Decay(EqualityMixin):
 _DECAY_PHOTON_ENERGY = {}
 
 
-def decay_photon_energy(nuclide: str) -> Univariate:
+def decay_photon_energy(nuclide: str) -> Optional[Univariate]:
     """Get photon energy distribution resulting from the decay of a nuclide
 
     This function relies on data stored in a depletion chain. Before calling it
@@ -595,9 +594,10 @@ def decay_photon_energy(nuclide: str) -> Univariate:
 
     Returns
     -------
-    openmc.stats.Univariate
-        Distribution of energies in [eV] of photons emitted from decay. Note
-        that the probabilities represent intensities, given as [decay/sec].
+    openmc.stats.Univariate or None
+        Distribution of energies in [eV] of photons emitted from decay, or None
+        if no photon source exists. Note that the probabilities represent
+        intensities, given as [decay/sec].
     """
     if not _DECAY_PHOTON_ENERGY:
         chain_file = openmc.config.get('chain_file')

openmc/material.py

@@ -92,10 +92,11 @@ class Material(IDManagerMixin):
     fissionable_mass : float
         Mass of fissionable nuclides in the material in [g]. Requires that the
         :attr:`volume` attribute is set.
-    decay_photon_energy : openmc.stats.Univariate
+    decay_photon_energy : openmc.stats.Univariate or None
         Energy distribution of photons emitted from decay of unstable nuclides
-        within the material. The integral of this distribution is the total
-        intensity of the photon source in [decay/sec].
+        within the material, or None if no photon source exists. The integral of
+        this distribution is the total intensity of the photon source in
+        [decay/sec].
 
         .. versionadded:: 0.14.0
 
@@ -264,7 +265,7 @@ class Material(IDManagerMixin):
         return density*self.volume
 
     @property
-    def decay_photon_energy(self) -> Univariate:
+    def decay_photon_energy(self) -> Optional[Univariate]:
         atoms = self.get_nuclide_atoms()
         dists = []
         probs = []
@@ -273,7 +274,7 @@ class Material(IDManagerMixin):
             if source_per_atom is not None:
                 dists.append(source_per_atom)
                 probs.append(num_atoms)
-        return openmc.data.combine_distributions(dists, probs)
+        return openmc.data.combine_distributions(dists, probs) if dists else None
 
     @classmethod
     def from_hdf5(cls, group: h5py.Group):

setup.py

@@ -5,11 +5,7 @@ import sys
 import numpy as np
 
 from setuptools import setup, find_packages
-try:
-    from Cython.Build import cythonize
-    have_cython = True
-except ImportError:
-    have_cython = False
+from Cython.Build import cythonize
 
 
 # Determine shared library suffix
@@ -76,13 +72,9 @@ kwargs = {
         'test': ['pytest', 'pytest-cov', 'colorama'],
         'vtk': ['vtk'],
     },
+    # Cython is used to add resonance reconstruction and fast float_endf
+    'ext_modules': cythonize('openmc/data/*.pyx'),
+    'include_dirs': [np.get_include()]
 }
 
-# If Cython is present, add resonance reconstruction and fast float_endf
-if have_cython:
-    kwargs.update({
-        'ext_modules': cythonize('openmc/data/*.pyx'),
-        'include_dirs': [np.get_include()]
-    })
-
 setup(**kwargs)

src/settings.cpp

@@ -60,6 +60,7 @@ bool run_CE {true};
 bool source_latest {false};
 bool source_separate {false};
 bool source_write {true};
+bool source_mcpl_write {true};
 bool surf_source_write {false};
 bool surf_mcpl_write {false};
 bool surf_source_read {false};
@@ -653,6 +654,9 @@ void read_settings_xml()
     if (check_for_node(node_sp, "write")) {
       source_write = get_node_value_bool(node_sp, "write");
     }
+    if (check_for_node(node_sp, "mcpl")) {
+      source_write = get_node_value_bool(node_sp, "mcpl");
+    }
     if (check_for_node(node_sp, "overwrite_latest")) {
       source_latest = get_node_value_bool(node_sp, "overwrite_latest");
       source_separate = source_latest;

src/simulation.cpp

@@ -389,16 +389,33 @@ void finalize_batch()
 
   if (settings::run_mode == RunMode::EIGENVALUE) {
     // Write out a separate source point if it's been specified for this batch
-    if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
-        settings::source_write && settings::source_separate) {
-      write_source_point(nullptr);
-    }
+#ifdef OPENMC_MCPL
+    if(! settings::source_mcpl_write) {
+#endif
+      if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
+          settings::source_write && settings::source_separate) {
+        write_source_point(nullptr);
+      }
 
-    // Write a continously-overwritten source point if requested.
-    if (settings::source_latest) {
-      auto filename = settings::path_output + "source.h5";
-      write_source_point(filename.c_str());
+      // Write a continously-overwritten source point if requested.
+      if (settings::source_latest) {
+        auto filename = settings::path_output + "source.h5";
+        write_source_point(filename.c_str());
+      }
+#ifdef OPENMC_MCPL
+    } else {
+      if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
+            settings::source_mcpl_write && settings::source_separate) {
+          write_mcpl_source_point(nullptr);
+      }
+
+      // Write a continously-overwritten source point if requested.
+      if (settings::source_latest && settings::source_mcpl_write) {
+        auto filename = settings::path_output + "source.mcpl";
+        write_mcpl_source_point(filename.c_str());
+      }
     }
+#endif
   }
 
   // Write out surface source if requested.

src/state_point.cpp

@@ -637,7 +637,7 @@ void write_mcpl_source_point(const char *filename, bool surf_source_bank)
 
   if (mpi::master) {
     //change this - this is h5 specific
-    mcpl_closeandgzip_outfile(file_id);
+    mcpl_close_outfile(file_id);
   }
 
 }

tests/regression_tests/source_mcpl_file/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-3.017557E-01 3.398770E-03
+3.039964E-01 3.654869E-04

tests/regression_tests/source_mcpl_file/settings.xml

@@ -1,11 +1,15 @@
 <?xml version="1.0"?>
 <settings>
+  <state_point batches="10" />
+  <source_point mcpl="true" separate="true" />
   <eigenvalue>
     <batches>10</batches>
     <inactive>5</inactive>
     <particles>1000</particles>
   </eigenvalue>
   <source>
-    <mcpl>source.10.mcpl</mcpl>
+    <space type="box">
+      <parameters>-4 -4 -4  4  4  4</parameters>
+    </space>
   </source>
 </settings>

tests/regression_tests/source_mcpl_file/test.py

@@ -1,65 +1,101 @@
-import os
-import subprocess
+#!/usr/bin/env python
+
 import glob
+import os
 
-from tests.testing_harness import TestHarness
+from tests.testing_harness import *
 
 
-class SourceMCPLFileTestHarness(TestHarness):
-  def execute_test(self):
-    """Run OpenMC with the appropriate arguments and check the outputs."""
-    try:
-      self._create_input()
-      self._test_input_created()
-      self._run_openmc()
-      self._test_output_created()
-      results = self._get_results()
-      self._write_results(results)
-      self._compare_results()
-    finally:
-      self._cleanup()
+settings1="""<?xml version="1.0"?>
+<settings>
+  <state_point batches="10" />
+  <source_point mcpl="true" separate="true" />
+  <eigenvalue>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <particles>1000</particles>
+  </eigenvalue>
+  <source>
+    <space type="box">
+      <parameters>-4 -4 -4  4  4  4</parameters>
+    </space>
+  </source>
+</settings>
+"""
 
-  def _create_input(self):
-    compiled=subprocess.run(['gcc','-o','gen_dummy_mcpl.out','gen_dummy_mcpl.c','-lm','-lmcpl'])
-    assert compiled==0, 'Could not compile mcpl-file generator code'
-    subprocess.run(['./gen_dummy_mcpl.out'])
+settings2 = """<?xml version="1.0"?>
+<settings>
+  <eigenvalue>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <particles>1000</particles>
+  </eigenvalue>
+  <source>
+    <mcpl> source.10.{0} </mcpl>
+  </source>
+</settings>
+"""
 
-  def update_results(self):
-    """Update the results_true using the current version of OpenMC."""
-    try:
-      self._create_input()
-      self._test_input_created()
-      self._run_openmc()
-      self._test_output_created()
-      results = self._get_results()
-      self._write_results(results)
-      self._overwrite_results()
-    finally:
-      self._cleanup()
 
-  def _test_input_created(self):
-    """Check that the input mcpl.file was generated as it should"""
-    mcplfile=glob.glob(os.path.join(os.getcwd(),'source.10.mcpl'))
-    assert len(mcplfile) == 1, 'Either multiple or no mcpl files ' \
-      'exist.'
-    assert mcplfile[0].endswith('mcpl'), \
-      'output file does not end with mcpl.'
+class SourceFileTestHarness(TestHarness):
+    def execute_test(self):
+        """Run OpenMC with the appropriate arguments and check the outputs."""
+        try:
+            self._run_openmc()
+            self._test_output_created()
+            self._run_openmc_restart()
+            results = self._get_results()
+            self._write_results(results)
+            self._compare_results()
+        finally:
+            self._cleanup()
 
-  def _test_output_created(self):
-    """Check that the output files were created"""
-    statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))
-    assert len(statepoint) == 1, 'Either multiple or no statepoint files ' \
-      'exist.'
-    assert statepoint[0].endswith('h5'), \
-      'statepoint file does not end with h5.'
+    def update_results(self):
+        """Update the results_true using the current version of OpenMC."""
+        try:
+            self._run_openmc()
+            self._test_output_created()
+            self._run_openmc_restart()
+            results = self._get_results()
+            self._write_results(results)
+            self._overwrite_results()
+        finally:
+            self._cleanup()
 
-  def _cleanup(self):
-    super()._cleanup()
-    source_mcpl=glob.glob(os.path.join(os.getcwd(),'source*.mcpl'))
-    for f in source_mcpl:
-      if (os.path.exists(f)):
-        os.remove(f)
+    def _test_output_created(self):
+        """Make sure statepoint and source files have been created."""
+        statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))
+        assert len(statepoint) == 1, 'Either multiple or no statepoint files ' \
+             'exist.'
+        assert statepoint[0].endswith('h5'), \
+             'Statepoint file is not a HDF5 file.'
 
-def test_mcpl_source_file():
-  harness = SourceMCPLFileTestHarness('source.10.mcpl')
-  harness.main()
+        source = glob.glob(os.path.join(os.getcwd(), 'source.10.mcpl*'))
+        assert len(source) == 1, 'Either multiple or no source files exist.'
+        assert source[0].endswith('mcpl') or source[0].endswith('mcpl.gz'), \
+             'Source file is not a MCPL file.'
+
+    def _run_openmc_restart(self):
+        # Get the name of the source file.
+        source = glob.glob(os.path.join(os.getcwd(), 'source.10.*'))
+
+        # Write the new settings.xml file.
+        with open('settings.xml','w') as fh:
+            fh.write(settings2.format(source[0].split('.')[-1]))
+
+        # Run OpenMC.
+        self._run_openmc()
+
+    def _cleanup(self):
+        TestHarness._cleanup(self)
+        output = glob.glob(os.path.join(os.getcwd(), 'source.*'))
+        #for f in output:
+        #    if os.path.exists(f):
+        #        os.remove(f)
+        with open('settings.xml','w') as fh:
+            fh.write(settings1)
+
+
+def test_source_file():
+    harness = SourceFileTestHarness('statepoint.10.h5')
+    harness.main()

tests/unit_tests/test_material.py

@@ -573,3 +573,9 @@ def test_decay_photon_energy():
     assert src.integral() == pytest.approx(sum(
         intensity(decay_photon_energy(nuc)) for nuc in m.get_nuclides()
     ))
+
+    # A material with no unstable nuclides should have no decay photon source
+    stable = openmc.Material()
+    stable.add_nuclide('Gd156', 1.0)
+    stable.volume = 1.0
+    assert stable.decay_photon_energy is None