Added new command line options. Now can run particle/batch/inactive optimization. Also added more depletion capability

BEAVRS/make_beavrs.py

@@ -1,14 +1,34 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 
-import openmc, os
+import json, openmc, os, shutil, sys
+import numpy as np
 import openmc.deplete
 from beavrs.builder import BEAVRS
 from optparse import OptionParser
+import matplotlib.pyplot as plt
 
 try:
-    os.mkdir('build')
+    os.mkdir(os.path.dirname(os.path.realpath(__file__)) + '/build')
 except OSError: pass
-os.chdir('build')
+os.chdir(os.path.dirname(os.path.realpath(__file__)) + '/build')
+
+def no_overwrite(filename):
+    if os.path.isfile(filename):
+        j = 0
+        new_filepath = filename*1
+        path, new_name = os.path.split(new_filepath)
+        name_len = len(new_name)
+        while os.path.isfile(new_filepath):
+            if new_name[-2] == 'h':
+                new_name = f'#{new_name}.{j}'
+            else:
+                j += 1
+                new_name = f'{new_name[:name_len+1]}.{j}'
+
+            new_filepath = f"{path}/{new_name}"
+
+        print(f"Backed up file {filename} to {new_filepath}")
+        os.replace(filename, new_filepath)
 
 usage = """usage: %prog [options]"""
 p = OptionParser(usage=usage)
@@ -20,38 +40,116 @@ p.add_option('-s', '--symmetric', action='store_true', dest='is_symmetric',
              + ' not symmetric by default')
 p.add_option('-e', '--deplete', action='store_true', dest='run_deplete', default=False,
              help='Run the depletion code for this setup')
+p.add_option('-i', '--insert', action='store', type='int', dest='S', default=200,
+             help='Set where the control rods should be. Acceptable values are 15-574')
+p.add_option('-o', '--optimal', action='store_true', default=False,
+             help='Find the best combination of particles, inactive and active batches to run')
 (options, args) = p.parse_args()
 
 if not len(args) == 0:
     p.print_help()
 
-chain_path = "/opt/xdata/endfb-vii.1-hdf5/chain_endfb71_pwr.xml"
-
 b = BEAVRS(is_symmetric=options.is_symmetric, is_2d=options.is_2d)
-b.write_openmc_geometry()
-b.write_openmc_materials()
-b.write_openmc_plots()
-b.write_openmc_settings()
+
+model = openmc.Model(geometry=b.openmc_geometry, materials=b.materials,
+                         settings=b.settings, tallies=b.common_tallies, plots=b.plots)
+
+model.materials.cross_sections = "/opt/xdata/endfb-vii.1-hdf5/cross_sections.xml"
+model.geometry.merge_surfaces = True
+model.settings.confidence_intervals = False
 
 # openmc.calculate_volumes()
+# openmc.plot_geometry()
+
+if options.optimal:
+
+    particles = [5000, 10000, 20000, 50000, 100000]
+    inactive = [10, 20, 50, 100, 150]
+    active = [50, 100, 150, 300, 500]
+
+    for p in particles:
+        for i in inactive:
+            for a in active:
+                print(f"Particles: {p} Inactive Cycles: {i} Active Cycles: {a}\n\n")
+                model.settings.batches = i+a
+                model.settings.inactive = i
+                model.settings.particles = p
+                model.export_to_xml()
+                openmc.run()
+                print("\n\n\n")
+                sys.stdout.flush()
 
 ###############################################################################
-#                      Depletion settings
+#                      Depletion Settings
 ###############################################################################
 
 if options.run_deplete:
-    model = openmc.model.Model(geometry=b.openmc_geometry, materials=b.materials, settings=b.settings, tallies=b.tallies, plots=plots)
+    model.settings.batches = 350
+    model.settings.inactive = 50
+    model.settings.particles = 20000
+    model.settings.seed = np.random.randint(1e14, dtype=np.uint64)
+    model.export_to_xml()
 
-    ## setting the transport operator
-    operator = openmc.deplete.CoupledOperator(model,chain_path,diff_burnable_mats=False,normalization_mode='fission-q',fission_q={"U235": 202.27e6})
+    # Get fission Q values from JSON file generated by get_fission_qvals.py
+    with open('../beavrs/serpent_fissq.json', 'r') as f:
+        serpent_fission_q = json.load(f)
+
+    chain_file = "/opt/xdata/endfb-vii.1-hdf5/chain_endfb71_pwr.xml"
+
+    try:
+        os.makedirs('depletion_results')
+    except FileExistsError:
+        pass
 
     ## setting the system linear power [W]
-    power = [45.50,91.01,136.51,182.02,227.52,273.03,324.45,324.45,273.03,273.03,227.52,182.02,136.51,91.01,45.50]
-    time_steps = [8640.00,8640.00,8640.00,8640.00,8640.00,8640.00,4268160.00,45792000.00,1728000.00,46915200.00,8640.00,8640.00,8640.00,8640.00,8640.00]
+    """powers = [45.50,91.01,136.51,182.02,227.52,273.03,324.45,324.45,273.03,273.03,227.52,182.02,136.51,91.01,45.50]
+    time_steps = [8640.00,8640.00,8640.00,8640.00,8640.00,8640.00,4268160.00,45792000.00,1728000.00,46915200.00,8640.00,8640.00,8640.00,8640.00,8640.00]"""
+
+    ## cumulative steps in days
+    time_steps_cum = np.array([3/24, 6/24, 9/24, 12/24, 15/24, 18/24, 21/24, 1, 10, 30, 60, 90, 365, 548])
+    time_steps = np.diff(time_steps_cum, prepend=0.0)
+
+
+    ## Power in Watts
+    full_power = 3411
+    powers = [full_power if x > 1 else full_power*x for x in time_steps_cum]
 
     print(time_steps)
 
-    ## depleting usin a first-order predictor algorithm
-    integrator = openmc.deplete.PredictorIntegrator(operator, time_steps, power, timestep_units =  's')
+    i = 0
+    for power, time_step in zip(powers, time_steps):
+
+        print(f'Current loop info: {power} {time_step}')
+        ## setting the transport operator
+        operator = openmc.deplete.CoupledOperator(model, chain_file,
+                                        #diff_burnable_mats=False, normalization_mode='fission-q',
+                                        fission_q=serpent_fission_q,
+                                       fission_yield_mode="average")
+
+        ## depleting using a first-order predictor algorithm
+        integrator = openmc.deplete.PredictorIntegrator(operator, [time_step], [power], timestep_units='d')
+        integrator.integrate()
+
+        results = openmc.deplete.Results('depletion_results.h5')
+        time, n_Xe135 = results.get_atoms('1', 'Xe135')
+        print(n_Xe135)
+        days = 24*60*60
+        plt.plot(time/days, n_Xe135)
+        plt.xlabel('Time [d]')
+        plt.ylabel('Xe135 [atoms]')
+
+        filename = f'depletion_results/depletion_results_t{i}.h5'
+        no_overwrite(filename)
+        os.replace('depletion_results.h5', filename)
+        shutil.copy2('materials.xml', f'depletion_results/materials_{i}.xml')
+        i += 1
+
+        # Get materials at the end of the last simulation
+        model.materials = results.export_to_materials(len(time_steps))
+
+
+
+
+
+
 
-    integrator.integrate()