diff --git a/BEAVRS/make_beavrs.py b/BEAVRS/make_beavrs.py index 210c2b8..9fe0eaf 100755 --- a/BEAVRS/make_beavrs.py +++ b/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()