Hi Zhiying,
I was able to get it to work with my code using the following method,
time_steps = [5/24] * 8
current_time = 0.0
for i, dt in enumerate(time_steps, start=1):
current_time += dt
step_id = f"{current_time:.10f}_days" #Assigns the step id to the number of days
model = openmc.Model(
geometry=geometry,
materials=materials,
settings=settings,
# tallies=tallies_reactionR_file
)
operator = openmc.deplete.CoupledOperator(model, "chain.xml", normalization_mode="energy-deposition")
integrator = openmc.deplete.PredictorIntegrator(operator, [time_steps[i-1]], power, timestep_units='d', solver='cram16').integrate()
rename_and_move_openmc_outputs(step_id)
filename = f"statepoint.20_{step_id}.h5"
# construct the statepoint path to call in later statepoints
statepoint_path = os.path.join(parent_folder_name, step_id, filename)
results = openmc.deplete.Results(f'{parent_folder_name}/{step_id}/depletion_results_{step_id}.h5')
materials = results.export_to_materials(1)
#Arguments for the boron criticality function
Boron_function_args = {
#Dimensions in cm
'r_fuel': r_fuel,
'r_clad': r_clad,
'pitch': pitch,
'pressure_tube_ir': pressure_tube_ir,
'pressure_tube_or': pressure_tube_or,
'calandria_ir': calandria_ir,
'calandria_or': calandria_or,
'fuel_bundle_length': fuel_bundle_length,
'sheath_length': sheath_length,
'r_coolant': r_coolant,
'Calandria_vessel_ir': Calandria_vessel_ir,
'Calandria_vessel_or': Calandria_vessel_or,
'Calandria_vessel_length': Calandria_vessel_length,
'reflective_distance_x': reflective_distance_x,
'reflective_distance_y': reflective_distance_y,
'void_z': void_z,
#Temperatures in Kelvin
'moderator_temp': moderator_temp,
'fuel_temp': fuel_temp,
'coolant_temp': coolant_temp,
'clad_temp': clad_temp,
'pressure_tube_temp': pressure_tube_temp,
'anulus_temp': anulus_temp,
'calandria_tube_temp': calandria_tube_temp,
'calandria_vessel_temp': calandria_vessel_temp,
'adjuster_temp': adjuster_temp,
'LCZ_temp': LCZ_temp,
#Densities in g/cm3
'moderator_dens': moderator_dens,
'fuel_dens': fuel_dens,
'coolant_dens': coolant_dens,
'clad_dens': clad_dens,
'pressure_tube_dens': pressure_tube_dens,
'anulus_dens': anulus_dens,
'calandria_tube_dens': calandria_tube_dens,
'calandria_vessel_dens': calandria_vessel_dens,
'adjuster_dens': adjuster_dens,
'LCZ_dens': LCZ_dens,
#D2O Atomic Concentraion in homogenized adjuster rods between 0 and 1
'D2OC': D2OC,
#Settings
'batches': batches,
#Need at least 100 entries removed, so 10 batches*10 gens, or 100 batches*1 gen
'inactive_batches': inactive_batches,
'generations': generations,
'particles_boron': particles_boron,
#Misc
'boron_enrichment': boron_enrichment,
'fuel_materials' : [],
'fuel_array' : fuel_array
}
#Parametized model for boron criticality search
def Boron_build_model(ppm_Boron, materials, **kwargs):
openmc.reset_auto_ids() # clears ID counters and registries
r_fuel = kwargs["r_fuel"]
r_clad = kwargs["r_clad"]
pitch = kwargs["pitch"]
pressure_tube_ir = kwargs["pressure_tube_ir"]
pressure_tube_or = kwargs["pressure_tube_or"]
calandria_ir = kwargs["calandria_ir"]
calandria_or = kwargs["calandria_or"]
fuel_bundle_length = kwargs["fuel_bundle_length"]
sheath_length = kwargs["sheath_length"]
r_coolant = kwargs["r_coolant"]
Calandria_vessel_ir = kwargs["Calandria_vessel_ir"]
Calandria_vessel_or = kwargs["Calandria_vessel_or"]
Calandria_vessel_length = kwargs["Calandria_vessel_length"]
reflective_distance_x = kwargs["reflective_distance_x"]
reflective_distance_y = kwargs["reflective_distance_y"]
void_z = kwargs["void_z"]
#Temperatures in Kelvin
moderator_temp = kwargs["moderator_temp"]
fuel_temp = kwargs["fuel_temp"]
coolant_temp = kwargs["coolant_temp"]
clad_temp = kwargs["clad_temp"]
pressure_tube_temp = kwargs["pressure_tube_temp"]
anulus_temp = kwargs["anulus_temp"]
calandria_tube_temp = kwargs["calandria_tube_temp"]
calandria_vessel_temp = kwargs["calandria_vessel_temp"]
adjuster_temp = kwargs["adjuster_temp"]
LCZ_temp = kwargs["LCZ_temp"]
#Densities in g/cm3
moderator_dens = kwargs["moderator_dens"]
fuel_dens = kwargs["fuel_dens"]
coolant_dens = kwargs["coolant_dens"]
clad_dens = kwargs["clad_dens"]
pressure_tube_dens = kwargs["pressure_tube_dens"]
anulus_dens = kwargs["anulus_dens"]
calandria_tube_dens = kwargs["calandria_tube_dens"]
calandria_vessel_dens = kwargs["calandria_vessel_dens"]
adjuster_dens = kwargs["adjuster_dens"]
LCZ_dens = kwargs["LCZ_dens"]
#D2O Atomic Concentraion in homogenized adjuster rods between 0 and 1
D2OC = kwargs["D2OC"]
#Settings
batches = kwargs["batches"]
#Need at least 100 entries removed, so 10 batches*10 gens, or 100 batches*1 gen
inactive_batches = kwargs["inactive_batches"]
generations = kwargs["generations"]
particles_boron = kwargs["particles_boron"]
#Misc
boron_enrichment = kwargs["boron_enrichment"]
#arrays used for fuel
fuel_materials = []
fuel_array = kwargs["fuel_array"]
# Remove old moderator completely
materials = openmc.Materials([m for m in materials if m.name != "Moderator"])
# Reconstruct a new moderator material
moderator = openmc.Material(4, "Moderator")
moderator.add_nuclide('H2', 1.995)
moderator.add_nuclide('H1', 0.005)
moderator.add_nuclide('O16', 1.0)
moderator.add_element('B', ppm_Boron, enrichment=boron_enrichment, enrichment_target='B10', percent_type='ao')
moderator.set_density('g/cm3', moderator_dens)
moderator.add_s_alpha_beta('c_D_in_D2O')
materials.append(moderator)
materials.cross_sections = '/home/smrcenter/openmc/endfb80/cross_sections.xml'
materials.export_to_xml()
#Making planes and cylinders for creating fuel geometry
surf_fuel = openmc.ZCylinder(r=r_fuel)
back_surf_fuel = openmc.ZPlane(-fuel_bundle_length/2)
front_surf_fuel = openmc.ZPlane(fuel_bundle_length/2)
back_surf_clad = openmc.ZPlane(-sheath_length/2)
front_surf_clad = openmc.ZPlane(sheath_length/2)
clad_outside_surface = openmc.ZCylinder(r=r_clad)
coolant_inner_surf = openmc.ZCylinder(r=r_coolant)
pt_inner = openmc.ZCylinder(r=pressure_tube_ir)
pt_outer = openmc.ZCylinder(r=pressure_tube_or)
calandria_inner = openmc.ZCylinder(r=calandria_ir)
calandria_outer = openmc.ZCylinder(r=calandria_or)
#Pitch box
left = openmc.XPlane(-pitch/2)
right = openmc.XPlane(pitch/2)
bottom = openmc.YPlane(-pitch/2)
top = openmc.YPlane(pitch/2)
back = openmc.ZPlane(-sheath_length/2)
front = openmc.ZPlane(sheath_length/2)
pitch_box = +left & -right & +bottom & -top & +back & -front
#Making Fuel cells
for i, code in enumerate(fuel_array):
exec(f"fuel_cell_{code} = openmc.Cell(fill=Fuel_{code}, region=-surf_fuel & +back_surf_fuel & -front_surf_fuel)")
exec(f"fuel_cell_{code}.temperature = fuel_temp")
# Create clad side cell for each material
exec(f"clad_side_cell_{code} = openmc.Cell(fill=Sheath, region= +surf_fuel & -clad_outside_surface & +back_surf_fuel & -front_surf_fuel)")
exec(f"clad_side_cell_{code}.temperature = clad_temp")
# Create clad back cell for each material
exec(f"clad_back_cell_{code} = openmc.Cell(fill=Sheath, region=-clad_outside_surface & -back_surf_fuel & +back_surf_clad)")
exec(f"clad_back_cell_{code}.temperature = clad_temp")
# Create clad front cell for each material
exec(f"clad_front_cell_{code} = openmc.Cell(fill=Sheath, region=-clad_outside_surface & +front_surf_fuel & -front_surf_clad)")
exec(f"clad_front_cell_{code}.temperature = clad_temp")
exec(f"coolant_cell_inner_{code} = openmc.Cell(fill=Coolant, region=-coolant_inner_surf & +back_surf_clad & -front_surf_clad)")
exec(f"coolant_cell_inner_{code}.temperature = coolant_temp")
exec(f"coolant_cell_outer_{code} = openmc.Cell(fill=Coolant, region=+coolant_inner_surf & -pt_inner & +back_surf_clad & -front_surf_clad)")
exec(f"coolant_cell_outer_{code}.temperature = coolant_temp")
exec(f"bundle_universe_{code} = openmc.Universe(cells=(coolant_cell_inner_{code}, coolant_cell_outer_{code}))")
exec(f"pin_universe_{code} = openmc.Universe(cells=(fuel_cell_{code}, clad_side_cell_{code}, clad_back_cell_{code}, clad_front_cell_{code}))")
num_pins = [1, 6, 12, 18]
angles = [0, 0, 15, 0]
ring_radii = np.array([0.0, 1.4885, 2.8755, 4.3305])
for i, (r, n, a) in enumerate(zip(ring_radii, num_pins, angles)):
for j in range(n):
# Determine location of center of pin
theta = (a + j/n*360.) * pi/180.
x = r*cos(theta)
y = r*sin(theta)
pin_boundary = openmc.ZCylinder(x0=x, y0=y, r=r_clad)
eval(f"coolant_cell_inner_{code}").region &= +pin_boundary
eval(f"coolant_cell_outer_{code}").region &= +pin_boundary
# Create each fuel pin -- note that we explicitly assign an ID so
# that we can identify the pin later when looking at tallies
pin = openmc.Cell(fill=eval(f"pin_universe_{code}"), region=-pin_boundary)
pin.translation = (x, y, 0)
eval(f"bundle_universe_{code}").add_cell(pin)
exec(f"bundle_{code} = openmc.Cell(fill=bundle_universe_{code}, region=-pt_inner)")
for i, code in enumerate(fuel_array):
exec(f"pressure_tube_{code} = openmc.Cell(fill=Pressure_tube, region=+pt_inner & -pt_outer & +back_surf_clad & -front_surf_clad)")
exec(f"pressure_tube_{code}.temperature = pressure_tube_temp")
exec(f"CO2_layer_{code} = openmc.Cell(fill=CO2, region =+pt_outer & -calandria_inner & +back_surf_clad & -front_surf_clad)")
exec(f"CO2_layer_{code}.temperature = anulus_temp")
exec(f"calandria_{code} = openmc.Cell(fill=Calandria_tube, region=+calandria_inner & -calandria_outer & +back_surf_clad & -front_surf_clad)")
exec(f"calandria_{code}.temperature = calandria_tube_temp")
for i, code in enumerate(fuel_array):
exec(f"pitch_cell_{code} = openmc.Cell(fill=Moderator, region= +calandria_outer & +left & -right & +bottom & -top & +back & -front)")
for i, code in enumerate(fuel_array):
bundle = eval(f"bundle_{code}")
pressure_tube = eval(f"pressure_tube_{code}")
CO2_layer = eval(f"CO2_layer_{code}")
calandria = eval(f"calandria_{code}")
pitch_cell = eval(f"pitch_cell_{code}")
exec(f"pitch_universe_{code} = openmc.Universe(cells=[bundle, pressure_tube, CO2_layer, calandria, pitch_cell])")
#Pitch
for code in (fuel_array):
Pitch_Universe = eval(f"pitch_universe_{code}")
globals()[code] = Pitch_Universe
lat3d = openmc.RectLattice()
lat3d.lower_left = (-57.15/2, -57.15/2, -24.765)
lat3d.pitch = (pitch, pitch, 49.53)
lat3d.universes = [
[[AAA, BBB],
[CCC, DDD]]]
#defining the reflective boundaries
reflective_distance_x = 28.575
reflective_distance_y = 28.575
void_z = 24.765
x_reflective_plane_left = openmc.XPlane(-reflective_distance_x, boundary_type = 'reflective')
x_reflective_plane_right = openmc.XPlane(+reflective_distance_x, boundary_type = 'reflective')
y_reflective_plane_bottom = openmc.YPlane(-reflective_distance_y, boundary_type = 'reflective')
y_reflective_plane_top = openmc.YPlane(+reflective_distance_y, boundary_type = 'reflective')
z_void_back = openmc.ZPlane(-void_z, boundary_type = 'reflective')
z_void_front = openmc.ZPlane(+void_z, boundary_type = 'reflective')
array_cell = openmc.Cell(fill = lat3d, region = +x_reflective_plane_left & -x_reflective_plane_right & +y_reflective_plane_bottom & -y_reflective_plane_top & +z_void_back & -z_void_front)
root_universe = openmc.Universe(cells=[array_cell])
geometry = openmc.Geometry()
geometry.root_universe = root_universe
geometry.export_to_xml()
bounds = [-reflective_distance_x, -reflective_distance_y, -void_z, reflective_distance_x, reflective_distance_y, void_z]
settings = openmc.Settings()
uniform_dist= openmc.stats.Box(bounds[:3], bounds[3:])
settings.source = openmc.IndependentSource(
space=uniform_dist,
)
settings.batches = batches
settings.inactive = inactive_batches
settings.particles = particles
settings.generations_per_batch = generations
settings.temperature['method'] = 'interpolation'
settings.temperature['tolerance'] = 400
settings.temperature['range'] = [0.0, 2400.0]
settings.output = {'tallies': True}
settings.export_to_xml()
chain = openmc.deplete.Chain.from_xml("/home/smrcenter/openmc/endfb80/chain_casl_pwr.xml")
chain.export_to_xml("chain.xml")
model = openmc.Model(geometry, materials, settings)
return model
crit_ppm, guesses, keffs = openmc.search_for_keff(lambda ppm,**kwargs: Boron_build_model(ppm, materials, **kwargs), bracket=[0, 5.0e-5], model_args=Boron_function_args, tol=5e-2, print_iterations=True, run_args={'output': False})
print('Critical Boron Concentration: {:.6e}'.format(crit_ppm))
ppm_Boron = crit_ppm
openmc.reset_auto_ids() # clears ID counters and registries
if os.path.exists("model.xml"):
os.remove("model.xml")
materials = openmc.Materials([m for m in materials if m.name != "Moderator"])
Moderator = openmc.Material(4, "Moderator")
# add nuclides and boron using crit_ppm here
Moderator.add_nuclide('H2', 1.995)
Moderator.add_nuclide('H1', 0.005)
Moderator.add_nuclide('O16', 1.0)
Moderator.add_element('B', ppm_Boron, enrichment=boron_enrichment,
enrichment_target='B10', percent_type='ao')
Moderator.set_density('g/cm3', moderator_dens)
Moderator.add_s_alpha_beta('c_D_in_D2O')
materials.append(Moderator)
materials.export_to_xml()
You can ignore the kwargs I use and some of the other stuff in my function like the renaming etc., they are just so I didn’t have to change parameters inside and outside the function whenever I changed something with my initial model, the key part that made it work was
# Remove old moderator completely
materials = openmc.Materials([m for m in materials if m.name != "Moderator"])
# Reconstruct a new moderator material
moderator = openmc.Material(4, "Moderator")
moderator.add_nuclide('H2', 1.995)
moderator.add_nuclide('H1', 0.005)
moderator.add_nuclide('O16', 1.0)
moderator.add_element('B', ppm_Boron, enrichment=boron_enrichment, enrichment_target='B10', percent_type='ao')
moderator.set_density('g/cm3', moderator_dens)
moderator.add_s_alpha_beta('c_D_in_D2O')
materials.append(moderator)
materials.cross_sections = '/home/smrcenter/openmc/endfb80/cross_sections.xml'
materials.export_to_xml()
Which is directly calling the materials that you exported from the results, but you need to rewrite the portion that will be adjusted by the criticality function, so in my case that was the boron concentration within the moderator, therefore I reconstructed my moderator. This should work for any material besides possibly your fuel, which you might need to do a little more work since you can’t just reconstruct it using your initial material composition. For your case adjusting the height of the control rods, I don’t know if you would need to construct a material again. The next important thing is including the lambda in the line,
crit_ppm, guesses, keffs = openmc.search_for_keff(lambda ppm,**kwargs: Boron_build_model(ppm, materials, **kwargs), bracket=[0, 5.0e-5], model_args=Boron_function_args,
tol=5e-2, print_iterations=True, run_args={'output': False})
In order to be able to call the materials into the function when it’s expecting only the two positional arguments ppm and **kwargs. Lastly, in order for it to properly update your materials outside of the function you need to reconstruct your material that you adjusted after the criticality search is done using the same code, otherwise it won’t update for the next depletion. I prefer to use the initial guess method as I find it faster, but I found it giving me errors where it has issues with it going into the negatives and cancelling, although I will switch to brentq instead of bisect for the bracket method. Let me know if this works for you, though you might have other issues if you’re having it adjust your geometry and not your material variables.
