Initial open-source commit

Author: aaronjeffreyolson

.gitignore

@@ -0,0 +1,10 @@
+**/*pyc
+**/*txt
+**/*dat
+**/*h5
+**/*png
+**/*pdf
+**/*npy
+**/*csv
+**/*~
+

Core/1D/MonteCarloParticleSolverpy.py

@@ -0,0 +1,208 @@
+#!usr/bin/env python
+import sys
+sys.path.append('/../../Core/Tools')
+from RandomNumberspy import RandomNumbers
+from Particlepy import Particle
+from ClassToolspy import ClassTools
+import numpy as np
+import matplotlib.pyplot as plt
+
+## \brief Multi-D base geometry class - useless by itself, but other multi-D geometry classes inherit from it
+# \author Aaron Olson, [email protected], [email protected]
+#
+# Base class for multi-D geometries
+#
+class Geometry_Base(ClassTools):
+    def __init__(self):
+        super(Geometry_Base,self).__init__()
+        self.flshowplot = False; self.flsaveplot = False
+    
+    def __str__(self):
+        return str(self.__dict__)
+
+    ## \brief Define outer boundaries of parallelapided geometry and solve volume
+    #
+    # \param[in] xbounds, ybounds, zbounds lists of two floats, boundaries of problem domain
+    # \returns initializes self.xbounds, self.ybounds, self.zbounds, and calcs self.Volume
+    def defineGeometryBoundaries(self,xbounds=None,ybounds=None,zbounds=None):
+        assert isinstance(xbounds,list) and isinstance(ybounds,list) and isinstance(zbounds,list)
+        assert len(xbounds)==2 and len(ybounds)==2 and len(zbounds)==2
+        assert isinstance(xbounds[0],float) and isinstance(xbounds[1],float) and xbounds[0]<=xbounds[1]
+        if self.Part.numDims<2: assert xbounds[0]==-np.inf and xbounds[1]==np.inf
+        assert isinstance(ybounds[0],float) and isinstance(ybounds[1],float) and ybounds[0]<=ybounds[1]
+        if self.Part.numDims<3: assert ybounds[0]==-np.inf and ybounds[1]==np.inf
+        assert isinstance(zbounds[0],float) and isinstance(zbounds[1],float) and zbounds[0]<=zbounds[1]
+        if self.GeomType=='Markovian':
+            if not ( xbounds[0]==ybounds[0] and xbounds[0]==zbounds[0] and xbounds[1]==ybounds[1] and xbounds[1]==zbounds[1] and -xbounds[0]==xbounds[1] ):
+                raise Exception("PlaybookMC's current 3D Markovian geometry implementation requires a cubic domain centered at the origin")
+
+        self.xbounds = xbounds
+        self.ybounds = ybounds
+        self.zbounds = zbounds
+
+        self.Volume = self.zbounds[1]-zbounds[0]
+        if self.Part.numDims>=2: self.Volume *= (self.xbounds[1]-xbounds[0])
+        if self.Part.numDims==3: self.Volume *= (self.ybounds[1]-ybounds[0])
+
+
+    ## \brief Define boundary conditions
+    #
+    # \param[in] xBCs, yBCs, zBCs lists of two strings, 'vacuum' or 'reflective'
+    # \returns initializes self.xBCs, self.yBCs, self.zBCs
+    def defineBoundaryConditions(self,xBCs,yBCs,zBCs):
+        assert isinstance(xBCs,list) and isinstance(yBCs,list) and isinstance(zBCs,list)
+        assert len(xBCs)==2 and len(yBCs)==2 and len(zBCs)==2
+
+        if self.Part.numDims<2: assert xBCs[0]==None and xBCs[1]==None
+        else                  :
+            assert xBCs[0]=='vacuum' or xBCs[0]=='reflective'
+            assert xBCs[1]=='vacuum' or xBCs[1]=='reflective'
+
+        if self.Part.numDims<3: assert yBCs[0]==None and yBCs[1]==None
+        else                  :
+            assert yBCs[0]=='vacuum' or yBCs[0]=='reflective'
+            assert yBCs[1]=='vacuum' or yBCs[1]=='reflective'
+
+        assert zBCs[0]=='vacuum' or zBCs[0]=='reflective'
+        assert zBCs[1]=='vacuum' or zBCs[1]=='reflective'
+
+        self.xBCs = xBCs
+        self.yBCs = yBCs
+        self.zBCs = zBCs
+        
+
+    ## \brief Defines Cross Sections
+    #
+    # \param[in] totxs, list of floats, list of total cross sections
+    # \param scatxs, list of floats, list of scattering cross sections
+    # \returns sets cross sections
+    def defineCrossSections(self,totxs=None,scatxs=None):
+        assert isinstance(totxs,list)
+        self.nummats = len(totxs)
+        for i in range(0,self.nummats):
+            assert isinstance(totxs[i],float) and totxs[i]>=0.0
+        self.totxs = totxs
+        self.Majorant = max(totxs)
+        if not scatxs==None:
+            assert isinstance(scatxs,list) and len(scatxs)==self.nummats
+            for i in range(0,self.nummats):
+                assert isinstance(scatxs[i],float) and scatxs[i]>=0.0 and totxs[i]>=scatxs[i]
+            self.scatxs = scatxs
+            self.absxs = []
+            self.scatrat = []
+            for i in range(0,self.nummats):
+                self.absxs.append( self.totxs[i]-self.scatxs[i] )
+                self.scatrat.append( self.scatxs[i]/self.totxs[i] )
+
+    ## \brief Associates random number object - must be instance of RandomNumbers class
+    #
+    # \param[in] Rng RandomNumbers object
+    # \returns initializes self.Rng
+    def associateRng(self,Rng):
+        assert isinstance(Rng,RandomNumbers)
+        self.Rng = Rng
+
+    ## \brief Associates Particle object - must be instance of Particle class
+    #
+    # \param[in] Part Particle object
+    # \returns initializes self.Part
+    def associatePart(self,Part):
+        assert isinstance(Part,Particle)
+        self.Part = Part
+
+    ## \brief Tests for a leakage or boundary reflection event, flags the first, evaluates the second
+    #
+    # \returns str 'reflect', 'transmit', 'sideleak', or 'noleak'
+    def testForBoundaries(self):
+        if self.isclose(self.Part.z,self.zbounds[0]):
+            self.Part.z = self.zbounds[0]
+            if   self.zBCs[0]==    'vacuum': return 'reflect'
+            elif self.zBCs[0]=='reflective': self.Part.reflect_z()
+        if self.isclose(self.zbounds[1],self.Part.z):
+            self.Part.z = self.zbounds[1]
+            if   self.zBCs[1]==    'vacuum': return 'transmit'
+            elif self.zBCs[1]=='reflective': self.Part.reflect_z()
+
+        if self.isclose(self.Part.x,self.xbounds[0]):
+            self.Part.x = self.xbounds[0]
+            if   self.xBCs[0]==    'vacuum': return 'sideleak'
+            elif self.xBCs[0]=='reflective': self.Part.reflect_x()
+        if self.isclose(self.xbounds[1],self.Part.x):
+            self.Part.x = self.xbounds[1]
+            if   self.xBCs[1]==    'vacuum': return 'sideleak'
+            elif self.xBCs[1]=='reflective': self.Part.reflect_x()
+
+        if self.isclose(self.Part.y,self.ybounds[0]):
+            self.Part.y = self.ybounds[0]
+            if   self.yBCs[0]==    'vacuum': return 'sideleak'
+            elif self.yBCs[0]=='reflective': self.Part.reflect_y()
+        if self.isclose(self.ybounds[1],self.Part.y):
+            self.Part.y = self.ybounds[1]
+            if   self.yBCs[1]==    'vacuum': return 'sideleak'
+            elif self.yBCs[1]=='reflective': self.Part.reflect_y()
+
+        return 'noleak'
+
+
+    ## \brief Returns distance to geometry boundary, only boundary for WMC-based approach
+    #
+    # \returns positive float, distance along streaming path to boundary
+    def calculateDistanceToBoundary(self):
+        if   self.isclose(self.Part.u,0.0): dbx = np.inf
+        elif self.Part.u>0.0              : dbx = (self.xbounds[1]-self.Part.x)/self.Part.u
+        else                              : dbx = (self.xbounds[0]-self.Part.x)/self.Part.u
+
+        if   self.isclose(self.Part.v,0.0): dby = np.inf
+        elif self.Part.v>0.0              : dby = (self.ybounds[1]-self.Part.y)/self.Part.v
+        else                              : dby = (self.ybounds[0]-self.Part.y)/self.Part.v
+
+        if   self.isclose(self.Part.w,0.0): dbz = np.inf
+        elif self.Part.w>0.0              : dbz = (self.zbounds[1]-self.Part.z)/self.Part.w
+        else                              : dbz = (self.zbounds[0]-self.Part.z)/self.Part.w
+
+        return min(dbx,dby,dbz)
+
+
+    ## \brief Samples new point, chooses to accept or reject, and chooses to absorb or scatter
+    #
+    # \returns drives sampling of new points and returns str 'reject', 'absorb', or 'scatter'
+    # need to change to allow for standard tracking or Woodcock instead of just Woodcock
+    # pseudocode: if woodcock, sample whether or not it is accepted. If it's not woodcock,
+    # we just want to sample absorb or scatter. Will also do that if collision is accepted for woodcock.
+    def evaluateCollision(self):
+        if self.trackingType == "Woodcock": # is there a way to tell that we are using woodcock tracking at this point in the code??? If not, I think I need two functions
+            self.samplePoint() #define new point
+            if self.Rng.rand()>self.totxs[self.CurrentMatInd]/self.Majorant: return 'reject' #accept or reject potential collision
+        return 'absorb' if self.Rng.rand()>self.scatrat[self.CurrentMatInd] else 'scatter' #choose collision type
+
+
+    ## \brief samples random point within defined geometry 
+    #
+    # \returns list of floats that are x, y, z values
+    def _generateRandomPoint(self):
+        x = self.Rng.uniform(self.xbounds[0],self.xbounds[1])
+        y = self.Rng.uniform(self.ybounds[0],self.ybounds[1])
+        z = self.Rng.uniform(self.zbounds[0],self.zbounds[1])        
+
+        return [x,y,z]
+
+
+
+
+
+
+
+    ## \brief Define mixing parameters (correlation length and material probabilities)
+    # To be provided by non-base class
+    def defineMixingParams(self,laminf=None,prob=None):
+        assert 1==0
+
+    ## \brief Initializes list of material types and coordinate of location, initializes xyz boundary locations
+    # To be provided by non-base class
+    def initializeGeometryMemory(self):
+        assert 1==0
+
+    ## \brief Samples a new point according to the selected rules
+    # To be provided by non-base class
+    def samplePoint(self):
+        assert 1==0

Core/1D/OneDSlabpy.py

@@ -0,0 +1,150 @@
+#!usr/bin/env python
+from Geometry_Basepy import Geometry_Base
+import numpy as np
+import scipy.spatial
+import bisect as bisect
+import collections
+try:
+    import pandas as pd
+except:
+    print("Pandas not available, some functionality may run slower")
+    pass
+
+
+## \brief Multi-D Box-Poisson geometry class
+# \author Emily Vu, [email protected], [email protected]
+# \author Aaron Olson, [email protected], [email protected]
+#
+# Class for multi-D Box-Poisson geometry 
+#
+class Geometry_BoxPoisson(Geometry_Base):
+    def __init__(self):
+        super(Geometry_BoxPoisson,self).__init__()
+        self.flshowplot = False; self.flsaveplot = False
+        self.kdtrees = False
+        self.GeomType = 'BoxPoisson'
+
+    def __str__(self):
+        return str(self.__dict__)
+
+    
+    ## \brief Defines mixing parameters (correlation length and material probabilities)
+    # Note: laminf and rhoB notation and formula following LarmierJQSRT2017 different mixing types paper
+    #
+    # \param[in] laminf, float, correlation length of Markovian mixing
+    # \param[in] prob, list of floats, list of material abundances (probabilities)
+    # \returns sets corr length, mat. probs, and hyperplane density rhoB
+    def defineMixingParams(self,laminf=None,prob=None):
+        # Assert material chord lengths and slab length and store in object
+        assert isinstance(laminf,float) and laminf>0.0
+        assert isinstance(prob,list) and len(prob)==self.nummats
+        for i in range(0,self.nummats): assert isinstance(prob[i],float) and prob[i]>=0.0 and prob[i]<=1.0
+        assert self.isclose( np.sum(prob), 1.0 )
+        self.laminf = laminf
+        self.prob   = prob
+        self.rhoB   = 2.0 / 3.0 / self.laminf
+
+                
+    ## \brief Initializes list of material types and coordinate of location, initializes xyz boundary locations
+    #
+    # Note: Numpy arrays don't start blank and concatenate, where lists do, therefore plan to use lists, but convert to arrays using np.asarray(l) if needed to use array format
+    #
+    # \returns initializes self.MatInds and self.xBoundaries, self.yBoundaries, self.zBoundaries
+    def initializeGeometryMemory(self):
+        self._sampleNumberOfHyperplanes()
+        self._sampleHyperplaneLocations()
+        
+
+    ## \brief compute number of pseudo-interfaces
+    #
+    # \computes number of pseudo-interfaces in x, y, and z direction for defined geometry 
+    def _sampleNumberOfHyperplanes(self):
+        xLength = abs( self.xbounds[0]-self.xbounds[1] );  xAve = self.rhoB * xLength
+        yLength = abs( self.ybounds[0]-self.ybounds[1] );  yAve = self.rhoB * yLength
+        zLength = abs( self.zbounds[0]-self.zbounds[1] );  zAve = self.rhoB * zLength
+        assert xLength==yLength and xLength==zLength  #logic currently only for a cube
+        self.xPInterfaces = self.Rng.poisson(xAve)
+        self.yPInterfaces = self.Rng.poisson(yAve)
+        self.zPInterfaces = self.Rng.poisson(zAve)
+        self.MatInds = np.full( ( self.xPInterfaces+1 , self.yPInterfaces+1 , self.zPInterfaces+1 ),None )
+        #print(self.xPInterfaces,self.yPInterfaces,self.zPInterfaces, "planes sampled in x, y, and z")
+
+
+    ## \brief sample location of pseudo-interfaces and constructs Box-Poisson geometry
+    #
+    # \returns boundary location vectors for each cartesianal direction 
+    def _sampleHyperplaneLocations(self):
+        self.xBoundaries = [self.xbounds[0],self.xbounds[1]]
+        self.yBoundaries = [self.ybounds[0],self.ybounds[1]]
+        self.zBoundaries = [self.zbounds[0],self.zbounds[1]]
+        
+        for x in range(0,self.xPInterfaces):
+            location = self.Rng.uniform(self.xbounds[0],self.xbounds[1])
+            self.xBoundaries.append(location)
+
+        for y in range(0,self.yPInterfaces):
+            location = self.Rng.uniform(self.ybounds[0],self.ybounds[1])
+            self.yBoundaries.append(location)
+
+        for z in range(0,self.zPInterfaces):
+            location = self.Rng.uniform(self.zbounds[0],self.zbounds[1])
+            self.zBoundaries.append(location)
+            
+        self.xBoundaries.sort()
+        self.yBoundaries.sort()
+        self.zBoundaries.sort()
+        
+    ## \brief Samples a new point according to the selected rules
+    #
+    # \returns scattering ratio of cell material type
+    def samplePoint(self):
+        #find index of nearest point
+        x = self.Part.x
+        y = self.Part.y
+        z = self.Part.z
+
+        xIndex = bisect.bisect(self.xBoundaries,x) - 1
+        yIndex = bisect.bisect(self.yBoundaries,y) - 1
+        zIndex = bisect.bisect(self.zBoundaries,z) - 1
+        #check if material has been sampled. If not sample material type.
+        self.CurrentMatInd = self.MatInds[xIndex,yIndex,zIndex]
+        if self.CurrentMatInd==None:
+            self.CurrentMatInd = int(self.Rng.choice(p=self.prob))
+            self.MatInds[xIndex,yIndex,zIndex] = self.CurrentMatInd
+        
+    def samplePointsFast(self, points):
+        #find index of nearest point 
+        xinds = np.searchsorted(self.xBoundaries, points[:,0], side='right')-1
+        yinds = np.searchsorted(self.yBoundaries, points[:,1], side='right')-1
+        zinds = np.searchsorted(self.zBoundaries, points[:,2], side='right')-1
+                
+        array_of_indices = np.vstack((xinds, yinds, zinds)).T
+        #check if material has been sampled. If not, sample material type.
+        matTypes = self.MatInds[xinds, yinds, zinds]
+        unassigned_inds = array_of_indices[matTypes == None]
+        if len(unassigned_inds) > 0:
+            nbx = len(self.xBoundaries)
+            nby = len(self.yBoundaries)
+            nbz = len(self.zBoundaries)
+            nbyz = nby*nbz            
+            single_ids = np.int64(xinds)*nbyz + np.int64(yinds)*nbz + np.int64(zinds)
+            try:
+                uniques = pd.unique(single_ids)
+            except:
+                uniques = np.unique(single_ids)
+            ynb_plus_z = uniques % nbyz
+            zi = ynb_plus_z % nbz
+            yi = (ynb_plus_z - zi)/nbz
+            xi = (uniques - ynb_plus_z)/nbyz
+            
+            xi = np.int16(xi)
+            yi = np.int16(yi)
+            zi = np.int16(zi)   
+            #uniques = np.array(list(uniques.keys()))
+            #uniques = np.unique(unassigned_inds, axis=0)
+            new_assignments = self.Rng.RngObj.choice(self.nummats, size=len(uniques), p=self.prob)            
+            #self.MatInds[uniques[:,0], uniques[:,1], uniques[:,2]] = new_assignments
+            self.MatInds[xi, yi, zi] = new_assignments
+            matTypes = self.MatInds[xinds, yinds, zinds]
+        return matTypes.astype(np.int16)
+