Address PR openmc-dev#618 comments

docs/source/usersguide/input.rst

@@ -290,7 +290,7 @@ OpenMC can use it for on-the-fly Doppler-broadening of resolved resonance range
 cross sections. If this element is absent from the settings.xml file, the
 :envvar:`OPENMC_MULTIPOLE_LIBRARY` environment variable will be used.
 
-  .. note:: The <use_windowed_multipole> element must also be set to "True"
+  .. note:: The <use_windowed_multipole> element must also be set to "true"
     for windowed multipole functionality.
 
 ``<max_order>`` Element

openmc/universe.py

@@ -475,7 +475,7 @@ def create_surface_elements(node, element):
         if self.temperature is not None:
             if isinstance(self.temperature, Iterable):
                 element.set("temperature", ' '.join(
-                            [str(t) for t in self.temperature]))
+                            str(t) for t in self.temperature))
             else:
                 element.set("temperature", str(self.temperature))
 

src/ace.F90

@@ -56,6 +56,7 @@ subroutine read_ace_xs()
     integer :: temp_table   ! temporary value for sorting
     character(12)  :: name  ! name of isotope, e.g. 92235.03c
     character(12)  :: alias ! alias of nuclide, e.g. U-235.03c
+    logical :: mp_found     ! if windowed multipole libraries were found
     type(Material),   pointer :: mat
     type(NuclideCE), pointer :: nuc
     type(SAlphaBeta), pointer :: sab
@@ -239,6 +240,21 @@ subroutine read_ace_xs()
       end if
     end do
 
+    ! If the user wants multipole, make sure we found a multipole library.
+    if (multipole_active) then
+      mp_found = .false.
+      do i = 1, n_nuclides_total
+        if (nuclides(i) % mp_present) then
+          mp_found = .true.
+          exit
+        end if
+      end do
+      if (.not. mp_found) call warning("Windowed multipole functionality is &
+           &turned on, but no multipole libraries were found.  Set the &
+           &<multipole_library> element in settings.xml or the &
+           &OPENMC_MULTIPOLE_LIBRARY environment variable.")
+    end if
+
   end subroutine read_ace_xs
 
 !===============================================================================
@@ -429,12 +445,12 @@ end subroutine read_ace_table
 
   subroutine read_multipole_data(i_table)
 
-    integer, intent(in) :: i_table   ! index in nuclides/sab_tables
+    integer, intent(in) :: i_table  ! index in nuclides/sab_tables
 
-    logical       :: file_exists   ! does multipole library exist?
-    character(7)  :: readable      ! is multipole library readable?
-    character(6)  :: zaid_string   ! String of the ZAID
-    character(MAX_FILE_LEN+9)  :: filename  ! path to multipole xs library
+    logical :: file_exists                 ! Does multipole library exist?
+    character(7) :: readable               ! Is multipole library readable?
+    character(6) :: zaid_string            ! String of the ZAID
+    character(MAX_FILE_LEN+9) :: filename  ! Path to multipole xs library
 
     ! For the time being, and I know this is a bit hacky, we just assume
     ! that the file will be zaid.h5.

src/input_xml.F90

@@ -167,8 +167,7 @@ subroutine read_settings_xml()
 
     ! Find the windowed multipole library
     if (run_mode /= MODE_PLOTTING) then
-      if (.not. check_for_node(doc, "multipole_library") .and. &
-           run_mode /= MODE_PLOTTING) then
+      if (.not. check_for_node(doc, "multipole_library")) then
         ! No library location specified in settings.xml, check
         ! environment variable
         call get_environment_variable("OPENMC_MULTIPOLE_LIBRARY", env_variable)
@@ -1109,9 +1108,9 @@ subroutine read_settings_xml()
     if (check_for_node(doc, "use_windowed_multipole")) then
       call get_node_value(doc, "use_windowed_multipole", temp_str)
       select case (to_lower(temp_str))
-      case ('true', 't', '1', 'y')
+      case ('true', '1')
         multipole_active = .true.
-      case ('false', 'f', '0', 'n')
+      case ('false', '0')
         multipole_active = .false.
       case default
         call fatal_error("Unrecognized value for <use_windowed_multipole> in &

src/multipole.F90

@@ -32,7 +32,6 @@ subroutine multipole_read(filename, multipole, i_table)
     integer :: i, j
     integer, allocatable :: MT(:)
     logical :: accumulated_fission
-    character(len=3) :: MT_string
     character(len=24) :: MT_n ! Takes the form '/nuclide/reactions/MT???'
     integer :: is_fissionable
 
@@ -91,13 +90,13 @@ subroutine multipole_read(filename, multipole, i_table)
       allocate(nuc % nu_fission(nuc % n_grid))
       allocate(nuc % absorption(nuc % n_grid))
 
-      nuc % total = ZERO
-      nuc % absorption = ZERO
-      nuc % fission = ZERO
+      nuc % total(:) = ZERO
+      nuc % absorption(:) = ZERO
+      nuc % fission(:) = ZERO
 
       ! Read in new energy axis (converting eV to MeV)
       call read_dataset(group_id, "energy_points", nuc % energy)
-      nuc % energy = nuc % energy / 1.0D6
+      nuc % energy = nuc % energy / 1.0e6_8
 
       ! Get count and list of MT tables
       call read_dataset(group_id, "MT_count", NMT)
@@ -111,20 +110,19 @@ subroutine multipole_read(filename, multipole, i_table)
 
       ! Loop over each MT entry and load it into a reaction.
       do i = 1, NMT
-        write(MT_string, '(I3.3)') MT(i)
-        MT_n = "/nuclide/reactions/MT" // MT_string
+        write(MT_n, '(A, I3.3)') '/nuclide/reactions/MT', MT(i)
 
         group_id = open_group(file_id, MT_n)
 
         ! Each MT needs to be treated slightly differently.
         select case (MT(i))
           case(ELASTIC)
             call read_dataset(group_id, "MT_sigma", nuc % elastic)
-            nuc % total = nuc % total + nuc % elastic
+            nuc % total(:) = nuc % total + nuc % elastic
           case(N_FISSION)
             call read_dataset(group_id, "MT_sigma", nuc % fission)
-            nuc % total = nuc % total + nuc % fission
-            nuc % absorption = nuc % absorption + nuc % fission
+            nuc % total(:) = nuc % total + nuc % fission
+            nuc % absorption(:) = nuc % absorption + nuc % fission
             accumulated_fission = .true.
           case default
             ! Search through all of our secondary reactions
@@ -136,36 +134,42 @@ subroutine multipole_read(filename, multipole, i_table)
                 ! fission cross section.
                 if ( (MT(i) == N_F .or. MT(i) == N_NF .or. MT(i) == N_2NF &
                       .or. MT(i) == N_3NF) .and. accumulated_fission) then
-                  nuc % total = nuc % total - nuc % fission
-                  nuc % absorption = nuc % absorption - nuc % fission
-                  nuc % fission = 0.0_8
+                  nuc % total(:) = nuc % total - nuc % fission
+                  nuc % absorption(:) = nuc % absorption - nuc % fission
+                  nuc % fission(:) = ZERO
                   accumulated_fission = .false.
                 end if
 
                 deallocate(nuc % reactions(j) % sigma)
                 allocate(nuc % reactions(j) % sigma(nuc % n_grid))
 
-                call read_dataset(group_id, "MT_sigma", nuc % reactions(j) % sigma)
-                call read_dataset(group_id, "Q_value", nuc % reactions(j) % Q_value)
-                call read_dataset(group_id, "threshold", nuc % reactions(j) % threshold)
+                call read_dataset(group_id, "MT_sigma", &
+                                  nuc % reactions(j) % sigma)
+                call read_dataset(group_id, "Q_value", &
+                                  nuc % reactions(j) % Q_value)
+                call read_dataset(group_id, "threshold", &
+                                  nuc % reactions(j) % threshold)
                 nuc % reactions(j) % threshold = 1 ! TODO: reconsider implications.
-                nuc % reactions(j) % Q_value = nuc % reactions(j) % Q_value / 1.0D6
+                nuc % reactions(j) % Q_value = nuc % reactions(j) % Q_value &
+                                               / 1.0e6_8
 
                 ! Accumulate total
                 if (MT(i) /= N_LEVEL .and. MT(i) <= N_DA) then
-                  nuc % total = nuc % total + nuc % reactions(j) % sigma
+                  nuc % total(:) = nuc % total + nuc % reactions(j) % sigma
                 end if
 
                 ! Accumulate absorption
                 if (MT(i) >= N_GAMMA .and. MT(i) <= N_DA) then
-                  nuc % absorption = nuc % absorption + nuc % reactions(j) % sigma
+                  nuc % absorption(:) = nuc % absorption &
+                                        + nuc % reactions(j) % sigma
                 end if
 
                 ! Accumulate fission (if needed)
                 if ( (MT(i) == N_F .or. MT(i) == N_NF .or. MT(i) == N_2NF &
                       .or. MT(i) == N_3NF) ) then
-                  nuc % fission = nuc % fission + nuc % reactions(j) % sigma
-                  nuc % absorption = nuc % absorption + nuc % reactions(j) % sigma
+                  nuc % fission(:) = nuc % fission + nuc % reactions(j) % sigma
+                  nuc % absorption(:) = nuc % absorption &
+                                        + nuc % reactions(j) % sigma
                 end if
               end if
             end do

tests/test_filter_distribcell/case-1/test.py

@@ -0,0 +1,7 @@
+import openmc
+
+
+su = openmc.Summary('summary.h5')
+sp = openmc.StatePoint('statepoint.1.h5')
+sp.link_with_summary(su)
+print(sp.tallies[1].get_pandas_dataframe(summary=su))

tests/test_multipole/test_multipole.py

@@ -1,5 +1,4 @@
 #!/usr/bin/env python
-
 import os
 import sys
 sys.path.insert(0, os.pardir)
@@ -9,7 +8,7 @@
 from openmc.source import Source
 
 
-class DistribmatTestHarness(PyAPITestHarness):
+class MultipoleTestHarness(PyAPITestHarness):
     def _build_inputs(self):
         ####################
         # Materials
@@ -116,8 +115,15 @@ def _build_inputs(self):
 
         plots_file.export_to_xml()
 
+    def execute_test(self):
+        if not 'OPENMC_MULTIPOLE_LIBRARY' in os.environ:
+            raise RuntimeError("The 'OPENMC_MULTIPOLE_LIBRARY' environment "
+                 "variable must be specified for this test.")
+        else:
+            super(MultipoleTestHarness, self).execute_test()
+
     def _get_results(self):
-        outstr = super(DistribmatTestHarness, self)._get_results()
+        outstr = super(MultipoleTestHarness, self)._get_results()
         su = openmc.Summary('summary.h5')
         outstr += str(su.get_cell_by_id(11))
         return outstr
@@ -126,9 +132,9 @@ def _cleanup(self):
         f = os.path.join(os.getcwd(), 'plots.xml')
         if os.path.exists(f):
             os.remove(f)
-        super(DistribmatTestHarness, self)._cleanup()
+        super(MultipoleTestHarness, self)._cleanup()
 
 
 if __name__ == '__main__':
-    harness = DistribmatTestHarness('statepoint.5.*')
+    harness = MultipoleTestHarness('statepoint.5.*')
     harness.main()