Address #1023 comments

docs/source/methods/photon_physics.rst

@@ -167,8 +167,8 @@ the two authors who discovered it:
     :label: klein-nishina
 
     \frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{\alpha'}{\alpha} \right
-    ) \left [ \frac{\alpha'}{\alpha} + \frac{\alpha}{\alpha'} + \mu^2 - 1 \right
-    ]
+    )^2 \left [ \frac{\alpha'}{\alpha} + \frac{\alpha}{\alpha'} + \mu^2 - 1
+    \right ]
 
 where :math:`\alpha` and :math:`\alpha'` are the ratios of the incoming and
 exiting photon energies to the electron rest mass energy equivalent (0.511 MeV),

docs/source/usersguide/beginners.rst

@@ -8,19 +8,19 @@ A Beginner's Guide to OpenMC
 What does OpenMC do?
 --------------------
 
-In a nutshell, OpenMC simulates neutral particles (presently only neutrons)
-moving stochastically through an arbitrarily defined model that represents an
-real-world experimental setup. The experiment could be as simple as a sphere of
-metal or as complicated as a full-scale `nuclear reactor`_. This is what's known
-as `Monte Carlo`_ simulation. In the case of a nuclear reactor model, neutrons
-are especially important because they are the particles that induce `fission`_
-in isotopes of uranium and other elements. Knowing the behavior of neutrons
-allows one to determine how often and where fission occurs. The amount of energy
-released is then directly proportional to the fission reaction rate since most
-heat is produced by fission. By simulating many neutrons (millions or billions),
-it is possible to determine the average behavior of these neutrons (or the
-behavior of the energy produced, or any other quantity one is interested in)
-very accurately.
+In a nutshell, OpenMC simulates neutral particles (presently neutrons and
+photons) moving stochastically through an arbitrarily defined model that
+represents an real-world experimental setup. The experiment could be as simple
+as a sphere of metal or as complicated as a full-scale `nuclear reactor`_. This
+is what's known as `Monte Carlo`_ simulation. In the case of a nuclear reactor
+model, neutrons are especially important because they are the particles that
+induce `fission`_ in isotopes of uranium and other elements. Knowing the
+behavior of neutrons allows one to determine how often and where fission
+occurs. The amount of energy released is then directly proportional to the
+fission reaction rate since most heat is produced by fission. By simulating
+many neutrons (millions or billions), it is possible to determine the average
+behavior of these neutrons (or the behavior of the energy produced, or any
+other quantity one is interested in) very accurately.
 
 Using Monte Carlo methods to determine the average behavior of various physical
 quantities in a system is quite different from other means of solving the same

docs/source/usersguide/scripts.rst

@@ -181,10 +181,10 @@ use with OpenMC. This script has the following optional arguments:
 ``openmc-get-photon-data``
 --------------------------
 
-This script downloads `ENDF/B-VII.1 <http://www.nndc.bnl.gov/endf/b7.1/zips/>`_
-ENDF data from NNDC for photo-atomic and atomic relaxation sublibraries and
-converts it to an HDF5 library for use with photon transport in OpenMC. This
-script has the following optional arguments:
+This script downloads `ENDF data <http://www.nndc.bnl.gov/endf/b7.1/zips/>`_
+from NNDC for photo-atomic and atomic relaxation sublibraries and converts it
+to an HDF5 library for use with photon transport in OpenMC. This script has the
+following optional arguments:
 
 -b, --batch
     Suppress standard in

openmc/data/library.py

@@ -74,9 +74,6 @@ class DataLibrary(EqualityMixin):
         ----------
         path : str
             Path to file to write. Defaults to 'cross_sections.xml'.
-        append : bool
-            Whether to append to an existing file, if it exists.
-            Defaults to False.
 
         """
         root = ET.Element('cross_sections')
@@ -87,10 +84,9 @@ class DataLibrary(EqualityMixin):
         if common_dir == '':
             common_dir = '.'
 
-        directory = os.path.relpath(common_dir, os.path.dirname(path))
-        if directory != '.':
+        if os.path.relpath(common_dir, os.path.dirname(path)) != '.':
             dir_element = ET.SubElement(root, "directory")
-            dir_element.text = directory
+            dir_element.text = os.path.realpath(common_dir)
 
         for library in self.libraries:
             lib_element = ET.SubElement(root, "library")

openmc/particle_restart.py

@@ -27,6 +27,8 @@ class Particle(object):
         Type of simulation (criticality or fixed source)
     id : long
         Identifier of the particle
+    type : int
+        Particle type (1 = neutron, 2 = photon, 3 = electron, 4 = positron)
     weight : float
         Weight of the particle
     energy : float
@@ -65,6 +67,10 @@ class Particle(object):
     def id(self):
         return self._f['id'].value
 
+    @property
+    def type(self):
+        return self._f['type'].value
+
     @property
     def n_particles(self):
         return self._f['n_particles'].value

openmc/statepoint.py

@@ -83,6 +83,8 @@ class StatePoint(object):
         Number of tally realizations
     path : str
         Working directory for simulation
+    photon_transport : bool
+        Indicate whether photon transport is active
     run_mode : str
         Simulation run mode, e.g. 'eigenvalue'
     runtime : dict
@@ -322,6 +324,10 @@ class StatePoint(object):
     def path(self):
         return self._f.attrs['path'].decode()
 
+    @property
+    def photon_transport(self):
+        return self._f.attrs['photon_transport'] > 0
+
     @property
     def run_mode(self):
         return self._f['run_mode'].value.decode()

src/nuclide_header.F90

@@ -994,6 +994,10 @@ contains
           micro_xs % fission         = ZERO
           micro_xs % nu_fission      = ZERO
         end if
+
+        ! Calculate microscopic nuclide photon production cross section
+        micro_xs % photon_prod = (ONE - f) * xs % value(XS_PHOTON_PROD,i_grid) &
+             + f * xs % value(XS_PHOTON_PROD,i_grid + 1)
       end associate
 
       ! Depletion-related reactions

src/particle_header.F90

@@ -337,6 +337,7 @@ contains
         call write_dataset(file_id, 'run_mode', 'particle restart')
       end select
       call write_dataset(file_id, 'id', this % id)
+      call write_dataset(file_id, 'type', this % type)
       call write_dataset(file_id, 'weight', src % wgt)
       call write_dataset(file_id, 'energy', src % E)
       call write_dataset(file_id, 'xyz', src % xyz)

src/particle_restart.F90

@@ -100,6 +100,7 @@ contains
       previous_run_mode = MODE_FIXEDSOURCE
     end select
     call read_dataset(p % id, file_id, 'id')
+    call read_dataset(p % type, file_id, 'type')
     call read_dataset(p % wgt, file_id, 'weight')
     call read_dataset(p % E, file_id, 'energy')
     call read_dataset(p % coord(1) % xyz, file_id, 'xyz')

src/physics.F90

@@ -567,7 +567,7 @@ contains
     integer, intent(in)  :: i_nuclide  ! index in nuclides array
     real(8), intent(in)  :: E          ! energy of neutron
     integer, intent(out) :: i_reaction ! index in nuc % reactions array
-    integer, intent(out) :: i_product  ! index in nuc % reactions array
+    integer, intent(out) :: i_product  ! index in reaction % products array
 
     integer :: i_grid
     integer :: i_temp
@@ -582,7 +582,7 @@ contains
     ! Get pointer to nuclide
     associate (nuc => nuclides(i_nuclide))
 
-      ! Get grid index and interpolation factor and sample proton production cdf
+      ! Get grid index and interpolation factor and sample photon production cdf
       i_temp = micro_xs(i_nuclide) % index_temp
       i_grid = micro_xs(i_nuclide) % index_grid
       f      = micro_xs(i_nuclide) % interp_factor
@@ -592,14 +592,13 @@ contains
       ! Loop through each reaction type
       REACTION_LOOP: do i_reaction = 1, size(nuc % reactions)
         associate (rx => nuc % reactions(i_reaction))
+          threshold = rx % xs(i_temp) % threshold
+
+          ! if energy is below threshold for this reaction, skip it
+          if (i_grid < threshold) cycle
+
           do i_product = 1, size(rx % products)
             if (rx % products(i_product) % particle == PHOTON) then
-
-              threshold = rx % xs(i_temp) % threshold
-
-              ! if energy is below threshold for this reaction, skip it
-              if (i_grid < threshold) cycle
-
               ! add to cumulative probability
               yield = rx % products(i_product) % yield % evaluate(E)
               prob = prob + ((ONE - f) * rx % xs(i_temp) % value(i_grid - threshold + 1) &
@@ -1725,7 +1724,8 @@ contains
     integer :: i
 
     ! Sample the number of photons produced
-    nu_t = micro_xs(i_nuclide) % photon_prod / micro_xs(i_nuclide) % total
+    nu_t =  p % wgt / keff * micro_xs(i_nuclide) % photon_prod / &
+         micro_xs(i_nuclide) % total
     if (prn() > nu_t - int(nu_t)) then
       nu = int(nu_t)
     else

src/state_point.F90

@@ -126,6 +126,11 @@ contains
       case (MODE_EIGENVALUE)
         call write_dataset(file_id, "run_mode", "eigenvalue")
       end select
+      if (photon_transport) then
+        call write_attribute(file_id, "photon_transport", 1)
+      else
+        call write_attribute(file_id, "photon_transport", 0)
+      end if
       call write_dataset(file_id, "n_particles", n_particles)
       call write_dataset(file_id, "n_batches", n_batches)
 
@@ -678,6 +683,12 @@ contains
     case ('eigenvalue')
       run_mode = MODE_EIGENVALUE
     end select
+    call read_attribute(int_array(1), file_id, "photon_transport")
+    if (int_array(1) == 1) then
+      photon_transport = .true.
+    else
+      photon_transport = .false.
+    end if
     call read_dataset(n_particles, file_id, "n_particles")
     call read_dataset(int_array(1), file_id, "n_batches")