Latin Hypercube Neutronic Sampling

.circleci/config.yml

@@ -11,6 +11,10 @@ jobs:
                  command: |
                      export APP_ENV
                      APP_ENV=docker
+                     pip3 install pyDOE
+                     cd neutronics/
+                     pytest
+                     cd ../optimization/
                      pytest
 workflows:  
     version: 2

neutronics/base_input.txt

@@ -1,4 +1,4 @@
-MCNP6 homog depltion r= ${r_core} frac= ${cool_frac}
+${model_information}
 c Cell Card
 1 1 -${fuel_rho} -1 4 -5 imp:n=1 
 2 2 -1.7 (-2 1 -7 6):(-1 5 -7):(6 -1 -4) imp:n=1 
@@ -15,7 +15,7 @@ c Surface Card
 c Data Card
 c
 c Burnup
-burn time=0.25 0.5 9R 10 5R 100 3R 365.25 8R power=${thermal_power} pfrac=1 29R 
+burn time=0.25 65 300 365.25 8R power=${core_power} pfrac=1 11R 
         bopt=1 14 -1 mat= 1
       omit= -1 128  66159 67163 67164 67166 68163 68165 68169 69166   
          69167 69171 69172 69173 70168 70169 70170 70171 70172 70173  

neutronics/mcnp_inputs.py

@@ -151,12 +151,15 @@ def write_mat_string(self):
         # write mcnp-form string
         self.fuel_string = self.homog_mat.mcnp().strip('\n')
 
-    def write_input(self):
+    def write_input(self, file_num, header):
         """ Write MCNP6 input files.
         This function writes the MCNP6 input files for the leakage experiment using
         the template input string. It writes a bare and reflected core input file
         for each core radius.
-
+
+        Arguments:
+        ----------
+            file_num (int): file number to define a filename
         Returns:
         --------
             filename (str): name of written MCNP6 input file
@@ -166,7 +169,8 @@ def write_input(self):
         self.write_mat_string()
         input_tmpl = open('base_input.txt')
         templ = Template(input_tmpl.read())
-        file_string = templ.substitute(cool_frac = self.vfrac_cermet,
+        file_string = templ.substitute(model_information = header,
+                                       cool_frac = self.vfrac_cermet,
                                        r_core = self.r,
                                        core_z = self.z,
                                        r_refl = self.r + self.refl_t,
@@ -175,11 +179,10 @@ def write_input(self):
                                        fuel_string = self.fuel_string,
                                        fuel_rho = self.rho,
                                        fuel_vol = self.core_vol,
-                                       refl_vol = self.core_vol,
-                                       thermal_power = self.Q_therm)
+                                       refl_vol = self.refl_vol,
+                                       core_power = self.Q_therm)
         # write the file
-        filename = 'r_{0}_{1}.i'.format(round(self.vfrac_cermet, 3), 
-                                                 round(self.r, 3))
+        filename = '{0}.i'.format(file_num)
         ifile = open(filename, 'w')
         ifile.write(file_string)
         ifile.close()

neutronics/neutronic_sweeps.py

@@ -0,0 +1,143 @@
+"""This module contains functions to sweep parameter-space, creating MCNP
+depletion inputs to calculate keff at EOL.
+* AR
+* core_z
+* cool_r
+* PD
+* power
+* enrich
+"""
+
+from pyDOE import lhs
+import os
+import glob
+import numpy as np
+import tarfile
+
+from mcnp_inputs import HomogeneousInput, build_pyne_matlib
+
+# reactor parameters to be sampled. Add parameters with caution
+_dimensions = ['AR', 'core_r', 'cool_r', 'PD', 'power', 'enrich']
+
+# number of samples for each sampled dimension
+nsamples = 100
+
+# Set dimension values. A tuple with range boundaries will sample the dimension
+# evenly using LHS, a float will set the dimension to the constant value
+dims = {'core_r'  : (20, 50),         
+        'AR'      : (0.7, 1.3),
+        'PD'      : (1.4, 1.6),        
+        'enrich'  : (0.3, 0.9),
+        'cool_r'  : 0.5,
+        'power'   : 150
+       }
+
+def get_sampling_params():
+    """Decide which parameters are constants and which are ranges to be sampled.
+    Convert constants (float) to ranges (tuples) with zero width.    
+    Returns:
+    --------
+        sampled (list): list of keys for sampled parameters
+    """
+    for item in _dimensions:
+        if type(dims[item]) != tuple:
+            dims[item] = (dims[item], dims[item])
+
+    sampled = _dimensions
+
+    return sampled
+
+def gen_hypercube(samples, N, seed=4654562):
+    """Generate N-dimensional latin hypercube to sample dimensional reactor
+    space.
+
+    Arguments:
+    ----------
+        samples (int): number of test cases to try
+        N (int): number of dimensions to test
+
+    Returns:
+    --------
+        cube (ndarray): normalized, N-D latin hypercube
+    """
+
+    np.random.seed(seed)
+    hypercube = lhs(N, samples=samples)
+
+    return hypercube
+
+def fill_data_array(swept, cube):
+    """Fill an ndarray with the sampling set generated by lhs and bounded by the
+    ranges provided for each sampled dimension.
+
+    Arguments:
+    ----------
+        swept (list): keys for sampled dimensions
+        const (list): keys for constant dimensions
+        cube (ndarray): latin hypercube 
+
+    Returns:
+    --------
+        rxt_confs (ndarray): array of data for every reactor configuration. This
+        data is used to build MCNP6 input files.
+    """
+    # initialize array
+    rxt_confs = np.zeros(nsamples, dtype={'names' : list(dims.keys()),
+                                          'formats' : ['f8']*len(dims)})
+
+    # for all samples in latin hypercube
+    for didx, parm in enumerate(_dimensions):
+        l_limit = dims[parm][0]
+        u_limit = dims[parm][1]
+        # uniform distribution
+        width = u_limit - l_limit
+
+        rxt_confs[parm] = np.add(l_limit, cube[:,didx]*width)
+
+    return rxt_confs
+
+def write_inputs(sampling_data):
+    """Write MCNP depletion inputs for each reactor configuration.
+
+    Arguments:
+    ----------
+        sampling_data (ndarray): array of reactor configurations, one for each
+        MCNP6 input file.
+    """
+    # build PyNE material library
+    pyne_mats = build_pyne_matlib()
+    # initialize tarball to save input files
+    tarputs = tarfile.open('smpl_mcnp_depl_inps.tar', 'w')
+    # write HTC input list
+    htc_inputs = open('input_list.txt', 'w')
+
+    # generate inputs
+    for num, sample in enumerate(sampling_data):
+        input = HomogeneousInput(sample['core_r'],
+                                 sample['core_r']*sample['AR'],
+                                 sample['power'], pyne_mats)
+        homog_comp = input.homog_core(sample['enrich'],
+                                      sample['cool_r'],
+                                      sample['PD'])
+        input.write_mat_string()
+
+        # identifying header string for post-processing
+        header_str = ''
+        for param in _dimensions:
+            header_str += str(round(sample[param], 5)) + ','
+        # write the input and tar it
+        filename = input.write_input(num, header_str)
+        tarputs.add(filename)
+        os.remove(filename)
+        htc_inputs.write('{0}.i\n'.format(num))
+
+    htc_inputs.close()
+    tarputs.add('input_list.txt')
+    tarputs.close()
+    os.remove('input_list.txt')
+
+if __name__=='__main__':
+    swept = get_sampling_params()
+    cube = gen_hypercube(nsamples, len(swept))
+    data = fill_data_array(swept, cube)
+    write_inputs(data)