Merge pull request #1 from pyne/develop

update my fork

Author: nathanb01

docs/install/vb.rst

@@ -17,3 +17,5 @@ the VirtualBox follow these steps:
      /home/pyne-user/pyne.
   #. Email any questions to [email protected].
 
+The default password for the user "pyne-user" is "pyne".
+

news/speed_up_source_sampling.rst

@@ -0,0 +1,16 @@
+**Added:** None
+
+**Changed:**
+
+* Pass cell_list back to Fortran, to speed up source sampling.
+* Update the test functions according to the cell_list change.
+
+**Deprecated:** None
+
+**Removed:** None
+
+**Fixed:**
+
+* Error in voxel R2S.
+
+**Security:** None

news/unstructured_r2s_step2_fix.rst

@@ -0,0 +1,14 @@
+**Added:** None
+
+**Changed:**
+
+* Fix a problem of creating mesh from reading h5m files in unstructued R2S
+* Update test function in test_r2s, allowing step2 run with the blank_mesh.h5m from step1
+
+**Deprecated:** None
+
+**Removed:** None
+
+**Fixed:** None
+
+**Security:** None

pyne/cpp_source_sampling.pxd

@@ -36,7 +36,7 @@ cdef extern from "source_sampling.h" namespace "pyne":
     cdef cppclass SourceParticle:
         # constructors
         SourceParticle() except +
-        SourceParticle(double, double, double, double, double, int) except +
+        SourceParticle(double, double, double, double, double, cpp_vector[int]) except +
 
         # attributes
 
@@ -46,7 +46,7 @@ cdef extern from "source_sampling.h" namespace "pyne":
         double get_z() except +
         double get_e() except +
         double get_w() except +
-        int get_c() except +
+        cpp_vector[int] get_cell_list() except +
 
 
 cdef extern from "source_sampling.h" namespace "pyne":

pyne/material.pyx

@@ -2314,6 +2314,8 @@ cdef class _MaterialLibrary(object):
             matsmetadata = f.get_node(datapath + '_metadata').read()
         for i in range(len(matstable)):
             row = matstable[i]
+            if len(nucs) == 0:
+                row = (row[0], row[1], row[2], [])
             comp = dict((<int> k, v) for k, v in zip(nucs, row[3]) if v != 0.0)
             mat = Material(comp, mass=row[0], density=row[1],
                                     atoms_per_molecule=row[2])

pyne/source_sampling.pyx

@@ -500,7 +500,7 @@ cdef class Sampler:
                 (<cpp_source_sampling.Sampler *> self._inst)\
                 .particle_birth(rands_proxy)
         return SourceParticle(c_src.get_x(), c_src.get_y(), c_src.get_z(), \
-                c_src.get_e(), c_src.get_w(), c_src.get_c())
+                c_src.get_e(), c_src.get_w(), c_src.get_cell_list())
 
 
 cdef class SourceParticle:
@@ -518,8 +518,8 @@ cdef class SourceParticle:
         The energy of the source particle
     w : double
         The weight of the source particle
-    c : int
-        The cell number of the source particle
+    cell_list : vector[int]
+        The cell list for available cells of the source particle
     
     Notes
     -----
@@ -528,12 +528,12 @@ cdef class SourceParticle:
     The class is found in the "pyne" namespace"""
 
     # constuctors
-    def __cinit__(self, double x, double y, double z, double e, double w, int c):
-        self.c_src = cpp_source_sampling.SourceParticle(x, y, z, e, w, c)
+    def __cinit__(self, double x, double y, double z, double e, double w, cpp_vector[int] cell_list):
+        self.c_src = cpp_source_sampling.SourceParticle(x, y, z, e, w, cell_list)
 
     @property
-    def c(self):
-        return self.c_src.get_c()
+    def cell_list(self):
+        return self.c_src.get_cell_list()
 
     @property
     def x(self):

share/source.F90

@@ -12,36 +12,96 @@
 ! Full instructions on compiling and using MCNP5 with this subroutine are found
 ! in the PyNE user manual.
 
-function find_cell() result(icl_tmp)
-! This function to determines the current MCNP cell index location and exits if
-! no valid cell is found. This only works if there are no repeated geometries or
-! universes present in the model.
+function find_cell(cell_list, cell_list_size) result(icl_tmp)
+! This function determines the current MCNP cell index location.
+! Return a positive integer if a valid cell is found, otherwise it returns -1.
+! This only works if there are no repeated geometries or universes present in
+! the model.
+! Parameters:
+!     cell_list: array of integers. Contains cell numbers that the source
+!                particle possiblely located.
+!     cell_list_size: integer. Size of the cell_list.
+!
+! There are 3 types of not found icl_tmp (icl_tmp == -1):
+! - Type 1: Sorce particle is located in a cell that exist in neutron transport
+!             but removed in photon transport. Therefore, the cell number does
+!             not exist in the ncl list.
+!           This is not an error, happens with small frequency.
+!           Skip it and resample next particle without warning message.
+! - Type 2: Source particle is not located in the given cell_list for HEX
+!             mesh cases (both voxel and sub-voxel). Because cell_list from
+!             discretize_geom() missed some cells with very small volume
+!             fractions.
+!           This is caused by random error, happens with small frequency.
+!           Skip it andd resample next particle with warning message.
+! - Type 3: Source partilce with a specific coordinates can't be found in any
+!             cell from ncl(1) to ncl(mxa).
+!           This is an error. When this error happens, it means that there is
+!             something wrong in either the source.F90 file or the DAGMC
+!             geometry representation. It happens with low frenquency for some
+!             complex geometry. The results is suspicious under this condition.
+!           Skip it and resample next particle with error message.
+
 
     use mcnp_global
     use mcnp_debug
     ! xxx,yyy,zzz are global variables
     ! mxa is global
     integer :: i ! iterator variable
-    integer :: j ! tempory cell test
+    integer :: j ! temporary cell test
     integer :: icl_tmp ! temporary cell variable
+    integer, intent(in) :: cell_list_size
+    integer, dimension(cell_list_size), intent(in) :: cell_list
+    integer :: cidx ! cell index
+
     icl_tmp = -1
+    ! If the cell_list is given (for HEX mesh),
+    ! use it to find cell first
+    if (cell_list_size > 0) then
+        ! HEX mesh. VOXEL/SUBVOXEL
+        do i = 1, cell_list_size
+           if (cell_list(i) .le. 0) then
+               ! not a valid cell number (-1)
+               exit
+           endif
+           cidx = namchg(1, cell_list(i))
+           if (cidx .eq. 0) then
+               ! Type 1: cell index not found, skip and resampling
+               exit
+           endif
+           call chkcel(cidx, 0, j)
+           if (j .eq. 0) then
+              ! valid cell found
+              icl_tmp = cidx
+              exit
+           endif
+        enddo
+    endif
 
-    do i = 1, mxa
-      call chkcel(i, 0, j)
-      if (j .eq. 0) then
-         ! valid cel set
-         icl_tmp = i
-         ! icl now set to be valid cell
-         exit
+    ! If the icl_tmp is not found yet (for HEX mesh), type 2 or type 3 happens,
+    ! or the cell_list is not given (for TET mesh),
+    ! find it in the entire list of cells
+    if ((icl_tmp == -1) .or. (cell_list_size .eq. 0)) then
+        do i = 1, mxa
+           call chkcel(i, 0, j)
+           if (j .eq. 0) then
+              ! valid cell found
+              icl_tmp = i
+              if (cell_list_size > 0) then
+                 ! this is a type 2 problem, skip
+                 ! reset the icl_tmp to -1 because of the type 2 not found
+                 icl_tmp = -1
+              endif
+              exit
+           endif
+        enddo
+      ! icl now is -1, it is a type 3 error.
+      ! Skip and print error message
+      if(icl_tmp .le. 0) then
+        write(*,*) 'ERROR: history ', nps, 'at position ', xxx, yyy, zzz, ' not in any cell'
+        write(*,*) 'Skipping and resampling the source particle'
       endif
-    enddo
-
-    ! icl now is -1
-    if(icl_tmp .le. 0) then
-      write(*,*) 'history ', nps, 'at position ', xxx, yyy, zzz, ' not in any cell'
-      write(*,*) 'Skipping and resampling the source particle'
     endif
-
 end function find_cell
 
 subroutine source
@@ -52,16 +112,19 @@ subroutine source
     implicit real(dknd) (a-h,o-z)
     logical, save :: first_run = .true.
     real(dknd), dimension(6) :: rands
-    integer :: icl_tmp ! temporary cell variable
+    integer :: icl_tmp ! temporary cell index variable
     integer :: find_cell
     integer :: tries
-    integer :: cell_num
+    integer, save :: cell_list_size = 0
+    integer, dimension(:), allocatable, save :: cell_list
 
     if (first_run .eqv. .true.) then
-        call sampling_setup(idum(1))
+        ! set up, and return cell_list_size to create a cell_list
+        call sampling_setup(idum(1), cell_list_size)
+        allocate(cell_list(cell_list_size))
         first_run = .false.
     endif
- 
+
 100 continue 
    tries = 0
    rands(1) = rang() ! sample alias table
@@ -72,31 +135,31 @@ subroutine source
    rands(4) = rang() ! sample y
    rands(5) = rang() ! sample z
 
-   cell_num = -1
-   call particle_birth(rands, xxx, yyy, zzz, erg, wgt, cell_num)
-   icl_tmp = find_cell()
+   call particle_birth(rands, xxx, yyy, zzz, erg, wgt, cell_list)
+   ! Loop over cell_list to find icl_tmp
+   icl_tmp = find_cell(cell_list, cell_list_size)
 
-   ! update tries and resampling if icl_tmp == -1
-   if (icl_tmp .eq. -1) then
-      tries = tries + 1
-      goto 200
+   ! check whether this is a valid cell
+   if (icl_tmp .le. 0) then
+      goto 300
    endif
-
-   if (idum(1) > 2) then
-       if (cell_num .ne. ncl(icl_tmp) .and. tries < idum(2)) then
-           tries = tries + 1
-           goto 200
-       endif
+   
+   ! check whether the material of sampled cell is void
+   if (mat(icl_tmp).eq.0) then
+       goto 300
+   else
+       goto 400
    endif
-   if (mat(icl_tmp).eq.0 .and. tries < idum(2)) then
-       tries = tries + 1
-       goto 200
-   end if
 
-   if(tries.eq.idum(2)) then
+300 continue
+   tries = tries + 1
+   if(tries < idum(2)) then
+       goto 200
+   else
        goto 100
-   end if
+   endif
 
+400 continue
    icl = icl_tmp
    tme = 0.0
    ipt = idum(3)