Updates to allow error handling for failed mcnp runs without killing program.

Author: jamesbevins

Code/Coeus_local.py

@@ -39,10 +39,10 @@
 import sys
 sys.path.insert(0,os.path.abspath(os.getcwd())+'/Sampling')
 
+#-------------------------------------------------------------------------------------------------------------#  
 
 def main(**kwargs):
 
-#-------------------------------------------------------------------------------------------------------------#  
 ### Local Function definitions
 
     # Print MCNP input Files

Code/Gnowee_Utilities.py

@@ -601,27 +601,33 @@ def Calc_Fitness(ids, pop, obj, min_fiss=0, max_w=1000):
     rundir=os.path.abspath(os.path.join(os.path.abspath(os.getcwd()),os.pardir))+"/Results/Population/"
 
     for i in ids:
+        tmp_fit=1E15
         index = next((c for c, parent in enumerate(pop) if parent.ident == i), -1)
         (tally,fissions,weight)=Read_MCNP_Output(rundir+str(i)+'/tmp/ETA.out', '24', '14')
-        tally=to_NormDiff(tally)
-        tmp_fit=RelativeLeastSquares(tally,obj)
-        module_logger.debug("Parent ID # {} has fitness = {} from RLS".format(i,tmp_fit))
-        
-        # Check constraints
-        weight=weight/1000         # conversion to kg
-        if fissions[0] == 0.0 and fissions[0] < min_fiss:
-            module_logger.warning("WARNING: No fissions occured for the ETA design in parent #{}".format(i))
-            tmp_fit+=1E15
-        elif fissions[0] > 0 and fissions[0] < min_fiss:
-            tmp_fit+= 0.1*(min_fiss/fissions[0]-1)
-        elif fissions[0] > min_fiss:
-            tmp_fit-= 0.01*(fissions[0]/min_fiss-1)
-        module_logger.debug("fissions[0] = {} and min_fiss = {} ".format(fissions[0],min_fiss))
-        module_logger.debug("Parent ID # {} has fitness = {} from RLS+fissions".format(i,tmp_fit))
+        try:
+            tally=to_NormDiff(tally)
+            tmp_fit=RelativeLeastSquares(tally,obj)
+            module_logger.debug("Parent ID # {} has fitness = {} from RLS".format(i,tmp_fit))
+
+            # Check constraints
+            weight=weight/1000         # conversion to kg
+            if fissions[0] == 0.0 and fissions[0] < min_fiss:
+                module_logger.warning("WARNING: No fissions occured for the ETA design in parent #{}".format(i))
+                tmp_fit+=1E15
+            elif fissions[0] > 0 and fissions[0] < min_fiss:
+                tmp_fit+= 0.1*(min_fiss/fissions[0]-1)
+            elif fissions[0] > min_fiss:
+                tmp_fit-= 0.01*(fissions[0]/min_fiss-1)
+            module_logger.debug("fissions[0] = {} and min_fiss = {} ".format(fissions[0],min_fiss))
+            module_logger.debug("Parent ID # {} has fitness = {} from RLS+fissions".format(i,tmp_fit))
+
+            if weight > max_w:
+                tmp_fit+=1E15
+            module_logger.info("Parent ID # {} has fitness = {} from RLS+fissions+weight".format(i,tmp_fit))
 
-        if weight > max_w:
-            tmp_fit+=1E15
-        module_logger.info("Parent ID # {} has fitness = {} from RLS+fissions+weight".format(i,tmp_fit))
+        except:
+            module_logger.info("WARNING: Parent ID # {} MCNP run failed.".format(i,tmp_fit))
+            pass
 
         # Save fitness
         pop[index].fit=tmp_fit  

Code/MCNP_Utilities.py

@@ -980,7 +980,7 @@ def Print_MCNP_Input(eta,geom,settings,mats,num,adv_print=True):
             inp_file.write("c  Physics  \n")
             inp_file.write("{}".format(settings.phys))
             inp_file.write("NPS {}\n".format(settings.nps))
-            inp_file.write("RAND GEN=2 STRIDE=1529171\n".format(settings.nps))
+            inp_file.write("RAND GEN=2 STRIDE=1529\n".format(settings.nps))
 
             # Print Material Cards
             inp_file.write("c ****************************************************************************\n")

Code/Utilities.py

@@ -509,7 +509,7 @@ def Run_Transport(lst,nps=[],code='mcnp6'):
             path=os.path.abspath(os.path.join(os.path.abspath(os.getcwd()), os.pardir))+"/Results/Population/"+str(i)+"/"
             sub.Popen("cp {} {}".format(path+"tmp/output/wwinp",path+"wwinp"),cwd=path,stdout=sub.PIPE,shell=True)
             sub.Popen("cp {} {}".format(path+"tmp/output/inp_edits.txt",path+"inp_edits.txt"), cwd=path,stdout=sub.PIPE,shell=True)
-            sub.Popen("rm -rf {}tmp/*".format(path,i),cwd=path,stdout=sub.PIPE,shell=True)
+            sub.Popen("rm -rf {}tmp/*".format(path,i),cwd=path,stderr=sub.STDOUT,stdout=sub.PIPE,shell=True)
 
     module_logger.info('Total transport time was {} sec'.format(time.time() - start_time))
 
@@ -611,6 +611,30 @@ def Build_Batch(lst,tasks,code,suf=""):
 
     return fname
 
+#-------------------------------------------------------------------------------------------------------------#  
+def to_Norm(spectrum):
+    """
+    Normalizes a MCNP tallied flux 
+   
+    Parameters
+    ==========
+    spectrum : array
+        The input flux spectrum
+
+    Optional
+    ========
+        
+    Returns
+    =======
+    result : array
+        The output normalized differential flux spectrum
+    """ 
+    
+    flux=np.zeros(len(spectrum[:,0]))
+    result=flux/np.sum(flux)
+        
+    return result
+
 #-------------------------------------------------------------------------------------------------------------#  
 def to_NormDiff(spectrum):
     """
@@ -639,25 +663,12 @@ def to_NormDiff(spectrum):
     # Calculate the differential flux
     diff[0]=(spectrum[0,1])/(spectrum[0,0])
     for i in range(1,len(spectrum[:,0])):
-        diff[i]=(spectrum[i,1]+spectrum[i-1,1])/(spectrum[i,0]-spectrum[i-1,0])
-    
-    # Calculate the integral of the differential flux using simpson's method
-    intdiff[0]=0.5*spectrum[0,0]*diff[0]
-    for i in range(1,len(spectrum[:,0])):
-        intdiff[i]=0.5*(spectrum[i,0]-spectrum[i-1,0])*(diff[i]+diff[i-1])
+        diff[i]=(spectrum[i,1])/(spectrum[i,0]-spectrum[i-1,0])
 
     # Calculate the normalized differential flux 
-    integral=np.sum(intdiff)
-    for i in range(0,len(spectrum[:,0])):
-        normdiff[i]=diff[i]/integral
-    
-    # Calculate the integral of the normalized differential flux using simpson's method
-    result[0]=0.5*spectrum[0,0]*normdiff[0]
-    for i in range(1,len(spectrum[:,0])):
-        result[i]=0.5*(spectrum[i,0]-spectrum[i-1,0])*(normdiff[i]+normdiff[i-1])
+    result=diff/np.sum(diff)
 
     return result
-
 
 #-------------------------------------------------------------------------------------------------------------#  
 def Uopt(c,d):
@@ -688,7 +699,7 @@ def Uopt(c,d):
 #-------------------------------------------------------------------------------------------------------------#  
 def LeastSquares(c,d):
     """
-    Calculate the U-optimality 
+    Calculate the LeastSquares fit between two arrays 
    
     Parameters
     ==========
@@ -706,7 +717,7 @@ def LeastSquares(c,d):
         The least-squares design based fitness
     """    
 
-    assert len(c)==len(d), "The length of the candidate and objective design must be equal in Uopt."  
+    assert len(c)==len(d), "The length of the candidate and objective design must be equal in LeastSquares."  
 
     return np.sum((d-c)**2)
 

Code/WorkingTemp.ipynb

@@ -280,6 +280,50 @@
     " "
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "2.2960312514\n"
+     ]
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "tally=np.array([7.25982E-07,5.20184E-08,5.79922E-08,8.65185E-08,7.74178E-08,1.00518E-07,9.05579E-08, 6.84944E-08,7.68031E-08,1.19829E-07,1.37603E-07,1.34688E-07,7.17418E-08,1.55737E-07,1.02621E-07, 5.98951E-07,4.16204E-07,4.80461E-07,7.74921E-07,6.96426E-07,3.03626E-07,2.48423E-07,1.16472E-07,2.03025E-07,1.80390E-07,2.21293E-07,1.77128E-07,1.26854E-07,1.68121E-08,8.70282E-08,7.81311E-08,3.20860E-08 ,9.61892E-09,4.57282E-08,2.56482E-08,1.65637E-08,1.85974E-08,1.84712E-08,2.07756E-08,2.32107E-08,7.91602E-09,7.37025E-08,8.14818E-08,2.89143E-09,6.85995E-11, 9.35566E-12])\n",
+    "obj=np.array([4.066918E-13,2.825704E-11,7.962365E-10,7.373641E-09,2.404237E-08,1.185623E-07,4.099819E-07,1.204311E-06,4.002197E-06,2.135014E-05,1.563630E-04,6.279676E-04,5.722449E-04,2.080694E-03,2.725402E-03,1.179853E-02,1.589037E-02,3.674752E-02,4.712247E-02,7.091333E-02,4.365288E-02,4.831482E-02,2.674840E-02,4.052126E-02,3.298227E-02,6.017881E-02,5.755246E-02,6.032581E-02,1.847843E-02,1.261844E-01,4.444340E-02,2.998172E-02,1.113570E-02,3.622714E-02,9.127742E-03,5.905254E-03,3.046808E-03,1.812024E-03,1.438887E-03,1.387365E-03,6.852790E-04,6.874538E-03,1.843700E-02,5.523352E-02,6.085823E-02,9.803895E-03])\n",
+    "\n",
+    "import Utilities as util\n",
+    "tmp_fit=util.RelativeLeastSquares(tally,obj)\n",
+    "\n",
+    "weight=1.687E6\n",
+    "fissions=[0,0]\n",
+    "\n",
+    "min_fiss=0\n",
+    "max_w=5000\n",
+    "\n",
+    "# Check constraints\n",
+    "weight=weight/1000         # conversion to kg\n",
+    "if fissions[0] == 0.0 and fissions[0] < min_fiss:\n",
+    "    tmp_fit+=1E15\n",
+    "elif fissions[0] > 0 and fissions[0] < min_fiss:\n",
+    "    tmp_fit+= 0.1*(min_fiss/fissions[0]-1)\n",
+    "elif fissions[0] > min_fiss:\n",
+    "    tmp_fit-= 0.01*(fissions[0]/min_fiss-1)\n",
+    "\n",
+    "if weight > max_w:\n",
+    "    tmp_fit+=1E15\n",
+    "    \n",
+    "print tmp_fit"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,

Inputs/eta_constraints_BNCT.csv

@@ -0,0 +1,39 @@
+Minimum Fissions,0.0
+ETA Max Weight, 5000
+Source Strength,1.5E15
+/
+TCC to ETA Distance, 6.0
+Debris Shield Thickness,1e-5
+ETA Wall Thickness,10
+Snout Distance, 71
+ETA Back Cover Thickness, 10
+ETA to Snout Mount Thickness,1e-5
+ETA Face Radius, 49
+ETA Cone Outer Radius, 50
+ETA Cone Opening Angle, 5.71
+/
+Debris Shield Material,Stainless Steel 409
+ETA Structural Material,Aluminum alloy 6061-O
+ETA Void Fill Material,Vacuum
+Fissile Foil,HEU
+/
+NAS Thickness, 1e-5
+NAS Radius, 1e-5
+NAS Material,Aluminum alloy 6061-O
+/
+NAS Activation Foils, Zr, Ni, In, Al, Ta
+NAS Activation Foil Thickness,1e-5, 1e-5, 1e-5, 1e-5, 1e-5
+NAS Activation Foil Radius, 1e-5
+TOAD Follows Material, In
+/
+TOAD Material,Aluminum alloy 6061-O
+TOAD Activation Foils,Al,Au,Al,HEU,Al
+TOAD Activation Foil Thickness,1e-5,1e-5,1e-5,1e-5,1e-5
+TOAD Activation Foil Radius,1e-5
+/
+Holder Material,Aluminum alloy 6061-O
+Holder Fill Material, Pb
+Holder Wall Thickness,1e-5
+/
+Max Vertical Components,3
+Max Horizontal Components,7

Inputs/mcnp_settings.csv

@@ -1,6 +1,6 @@
-Physics:
-MODE n
-/
-NPS:
-1E6
-/
+Physics:
+MODE n
+/
+NPS:
+1E6
+/