PWR Pin 2D Axial Symmetry - Solution
[1]:
import os
import time
import sys
mesh_path = '../../../mesh/'
benchmark_path = '../../../BenchmarkData/'
sys.path.append('../../../models/')
from ofelia.aux_obj import *
import matplotlib.pyplot as plt
from matplotlib import cm
Setting OpenMC materials
[2]:
# from ofelia.constants import *
from ofelia.materials import * # importing parameters and OpenMC materials
path_res = './results_'+str(n_div)+'div_'+str(Power/1000)+'kW/'
if not os.path.exists(path_res):
os.system("mkdir "+path_res)
print('The number of axial divisions is '+ str(n_div))
print('The pin power is set to ' +str(Power/1000)+' kW')
xml_materials = updateXML(mat_dict, n_div)
The number of axial divisions is 8
The pin power is set to 65.0 kW
Setting simulation parameters for FEniCSx
[3]:
from ofelia.thermal import *
from dolfinx.io import gmshio
from mpi4py import MPI
import gmsh
from IPython.display import clear_output
gdim = 2
model_rank = 0
mesh_comm = MPI.COMM_WORLD
mesh_factor = 0.05
# Initialize the gmsh module
gmsh.initialize()
# Load the .geo file
gmsh.merge(mesh_path+'pin2D_pwr.geo')
gmsh.model.geo.synchronize()
gmsh.option.setNumber("Mesh.MeshSizeFactor", mesh_factor)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.optimize("Netgen")
# Domain
domain, ct, ft = gmshio.model_to_mesh(gmsh.model, comm = mesh_comm, rank = model_rank, gdim = gdim )
gmsh.finalize()
clear_output()
fdim = gdim - 1
domain.topology.create_connectivity(fdim, gdim)
region_markers = [1, 2, 3] # 1: fuel, 2: helium, 3: cladding
# Boundary Markers
water_out_marker = 10
inl_marker = 20
out_marker = 30
# The thermal conducticity of the helium gap is increased to consider the radiation HT
th_input = thermal_inputs(T_w_average, pressure)
physical_param = {'th_cond': np.array([3, 0.25 * 3, 16.23]) / 100., # W/cm-K
'htc': th_input.compute_htc(pitch, clad_or, u_in)}
# Defining and assembling the linear structures for the Thermal Calculations
TH = thermal_solver(domain, ct, ft,
physical_param, region_markers, water_out_marker)
TH.assemble()
Vnorms = norms(TH.funSpace, domain)
RUN
[4]:
os.system("cat ../../../models/ofelia/header")
print("\n -- Simulation Start --\n")
# Picard settings initialization
i = 0
error = 100
alpha = 1
s1 = s1_val
sum_s = s1
# Guess temperature for the Picard Loop
Tguess_w = Tin
Tguess_f = 273.15 + 600
Tguess = []
Tguess.append(np.ones((n_div,)) * Tguess_f)
Tguess.append(np.ones((n_div,)) * Tguess_w)
# Defining class to convert power from OpenMC to FEniCSx
ExtractPower = extract_power(n_div, Power, meshSize, pin_length, pin_r, J_to_eV, tallyDict)
# The division in slices acts only on the active region
slices = np.linspace(-l_active/2, l_active/2, n_div+1)
# regLogic = parameters.assign(is_fuel + 2 * is_helium + 3 * is_clad + 4 * is_water + 5 * is_graph + 273.15)
# TH.plot2D(regLogic, L, pitch, 'pictures/regions')
#Define input folder path
path = os.getcwd()
phiE_list = [[], []] # mean, unc
q3_relaxed_list = []
q3_std_list = []
z_list = []
alpha_list = [1]
pop_list = [s1]
errors_q3_list = [1]
T_old = Function(TH.funSpace)
T_old.interpolate(lambda x: x[0] * 0. + Tin)
## Set computational time
comp_time = dict()
start_time = time.time()
start_cpu_time = time.process_time()
# os.system('python ../../../models/ofelia/build_xml.py')
from ofelia.build_xml import *
RemoveFolders(path=os.getcwd())
print('Tguess: fuel = {:.2f}'.format(np.mean(Tguess[0])-273.15)+' °C, water = {:.2f}'.format(np.mean(Tguess[1])-273.15)+' °C')
print('Tolerance for convergence: {:.2e}'.format(tol))
#print('Under relaxation constant: {:.2f}'.format(alpha))
print(' ')
while error > tol:
print(f'Iter #{i+0:02}')
input_folder = path + '/build_xml/it' + str(i)
#Run the simulation
print(' Running OpenMC')
openmc.run(cwd=input_folder, output = False)
#Rename the output with the info
old_sp_name = 'statepoint.'+str(batches)+'.h5'
old_sp_path = input_folder + '/' +old_sp_name
new_sp_name = 'statepoint.'+str(batches)+'.it' + str(i) + '.h5'
new_sp_path = input_folder + '/' +new_sp_name
os.system('mv ' + old_sp_path + ' ' + new_sp_path)
#Collecting the output (for simplicity, the k)
sp=openmc.StatePoint(new_sp_path)
# Get the energy flux from the iteration
phiE_list, EnergyStairs = ExtractPower.getSpectrum(sp, phiE_list)
# Get the axial power density from the iteration
q3_unrelaxed, z, q3std, A = ExtractPower.eval(sp, i)
q3_integral = np.trapz(q3_unrelaxed, z)
uq3_integral = np.trapz(q3std, z)
if i == 0:
print('\n Input Pin Power: {:.2f}'.format(Power/1000) + ' kW')
print(' Pin Power calculated (MC): {:.4f}'.format(A*q3_integral/1000) + ' +- {:.2f}'.format(A*uq3_integral/1000) +' kW\n')
q3_relaxed_list.append(q3_unrelaxed)
q3_relaxed = q3_unrelaxed
q3_std_list.append(q3std)
z_list.append(z)
# Under Relaxation
sn = (s1+ np.sqrt(s1*s1 + 4 * s1 * sum_s))/2
sum_s=sn
error = 1
if i>0:
sn = (s1+ np.sqrt(s1*s1 + 4 * s1 * sum_s))/2
sum_s = sum_s + sn
alpha = sn/sum_s
alpha_list.append(alpha)
pop_list.append(sn)
q3_relaxed = q3_unrelaxed * alpha + q3_relaxed_list[i-1] * (1-alpha)
errors_q3 = np.sqrt(np.trapz( (q3_relaxed - q3_relaxed_list[i-1])**2, z) / np.trapz( (q3_relaxed)**2, z) )
errors_q3_list.append(errors_q3)
q3_relaxed_list.append(q3_relaxed)
q3_std_list.append(q3std)
#Print the k value at the iteration
k_eff = sp.keff
keff = k_eff.n
ukeff = k_eff.std_dev
print(' keff: {:.6f}'.format(keff) + ' +- {:.6f}'.format(ukeff))
print(' alpha = {:.4f}'.format(alpha))
print(' Next pop size = ' + str(round(sn))+' (neutrons/cycle)\n')
# np.savetxt(path + '/build_xml/it' + str(i) + '/q3_it' + str(i)+'.txt', np.vstack([q3, z]))
# power_fun = importq3(path + '/build_xml/it' + str(i) + '/q3_it' + str(i)) # W/cm3
print(' Running FEniCSx')
power_fun, Tb_fun, res_bulk = th_input.mapping_q3_Tb(z, q3_relaxed, Tin, total_length, l_active, fuel_or)
T_sol = TH.solve(power_fun, Tb_fun)
# Computing regional average temperature
average_T = []
average_T.append(TH.computeSolidAverageT(region_markers[0], T_sol, slices))
dimz_full = z.shape[0]
if z.shape[0] < 100:
dimz_full = 100
else:
dimz_full = z.shape[0]+100
average_T.append(th_input.computeWaterAverageT(slices, dimz_full))
# Computing regional errors
errors_mat = np.zeros((len(average_T,)))
for mat_ID in range(len(average_T)):
errors_mat[mat_ID] = np.linalg.norm(average_T[mat_ID] - Tguess[mat_ID])/np.linalg.norm(Tguess[mat_ID])
# Check convergence: maxIter vs error < tol?
i += 1
error = np.mean(errors_mat)
error_L2 = Vnorms.L2norm(T_sol - T_old) / Vnorms.L2norm(T_old)
Tguess = average_T.copy()
T_old.x.array[:] = T_sol.x.array
print(f' Temperature regional average error: {float(error):.2e} | relative L2 error: {float(error_L2):.2e}')
if i > 1:
print(f' Power density: relative L2 error {float(errors_q3):.3e}')
error += errors_q3
print(f' Average T_fuel: {float(np.mean(Tguess[0])-273.15):.2f} °C | Average T_water: {float(np.mean(Tguess[1])-273.15):.2f} °C\n')
error = errors_q3_list[i-1]
if errors_q3_list[i-1] > tol:
xml_materials.update(Tguess[0], Tguess[1], sn, i)
if i >= maxIter:
error = 0
print(' ')
print('###############################################################')
print(' ')
print('Warning: Maximum Number of iteration reached!')
##################################################################
## ##
## OOO FFFFF EEEEEE L III A ##
## O O F E L I A A ##
## O O F E L I AAAAA ##
## O O FFFF EEEEEE L I A A ##
## O O F E L I AAAAAAA ##
## O O F E L I A A ##
## OOO F EEEEEEE LLLLL III A A ##
## ##
##################################################################
## ##
## Openmc-FEnicsx for muLtiphysics tutorIAl ##
## ------- ##
## ERMETE Lab - Politecnico di Milano ##
## ##
## GitHub: https://github.com/ERMETE-Lab/MP-OFELIA ##
## ##
##################################################################
-- Simulation Start --
Tguess: fuel = 600.00 °C, water = 290.00 °C
Tolerance for convergence: 5.00e-02
Iter #00
Running OpenMC
Input Pin Power: 65.00 kW
Pin Power calculated (MC): 67.0096 +- 1.75 kW
keff: 1.380605 +- 0.001630
alpha = 1.0000
Next pop size = 1214 (neutrons/cycle)
Running FEniCSx
Temperature regional average error: 1.60e-01 | relative L2 error: 8.60e-01
Average T_fuel: 741.78 °C | Average T_water: 307.00 °C
Iter #01
Running OpenMC
keff: 1.360200 +- 0.001300
alpha = 0.5357
Next pop size = 1400 (neutrons/cycle)
Running FEniCSx
Temperature regional average error: 9.17e-02 | relative L2 error: 1.83e-01
Power density: relative L2 error 3.985e-01
Average T_fuel: 744.00 °C | Average T_water: 310.47 °C
Iter #02
Running OpenMC
keff: 1.355326 +- 0.001248
alpha = 0.4111
Next pop size = 1824 (neutrons/cycle)
Running FEniCSx
Temperature regional average error: 1.27e-02 | relative L2 error: 2.56e-02
Power density: relative L2 error 5.613e-02
Average T_fuel: 743.73 °C | Average T_water: 309.98 °C
Iter #03
Running OpenMC
keff: 1.354819 +- 0.001075
alpha = 0.3352
Next pop size = 2238 (neutrons/cycle)
Running FEniCSx
Temperature regional average error: 8.20e-03 | relative L2 error: 1.67e-02
Power density: relative L2 error 3.630e-02
Average T_fuel: 744.02 °C | Average T_water: 310.21 °C
Post-Processing and data storing
[5]:
import pickle
from dolfinx.io import XDMFFile
plt.rcParams.update({
"text.usetex": True,
"font.family": "serif"
})
plt.rcParams['text.latex.preamble'] = r'\usepackage{amssymb} \usepackage{amsmath} \usepackage{amsthm} \usepackage{mathtools}'
comp_time['Wall time'] = (time.time() - start_time) / 60
comp_time['CPU time'] = (time.process_time() - start_cpu_time) / 60
print('Computational wall time = {:.4f}'.format(comp_time['Wall time'])+' minutes')
print('Computational CPU time = {:.4f}'.format(comp_time['CPU time'])+' minutes')
# Temperature storage in .h5 format
with XDMFFile(domain.comm, path_res+'T_sol.xdmf', "w") as loc:
T_sol.name = "T"
loc.write_mesh(domain)
loc.write_function(T_sol)
# Data store using numpy arrays
res2d = TH.extract_2D_data(T_sol, total_length, clad_or)
with open(path_res+'results.pin', 'wb') as f:
pickle.dump([comp_time, res_bulk, Tguess, res2d], f)
# Write axial power to txt
zWrite = np.linspace(-total_length/2, total_length/2, int(1e3))
q_to_write = np.zeros((len(zWrite), 2))
q_to_write[:, 0] = zWrite[:].flatten()
q_to_write[:, 1] = power_fun(zWrite, 0.).flatten()
np.savetxt(path_res+'axial_powerDensity.txt', q_to_write, delimiter=',')
# Write L2 error in taking the axial steps
# step_errors = np.zeros((4,))
# step_errors[0], step_errors[2] = TH.FuelAxError(T_sol, average_T, L_active, fuel_or, slices)
# step_errors[1], step_errors[3] = TH.WaterAxError(T_sol, average_T, L_active, clad_or, pitch, slices)
# np.savetxt('T_L2error_axsteps.txt', step_errors, delimiter=',')
# Under relaxation plots
fig = plt.figure (figsize = (6,4))
plt.plot(range(len(pop_list)), pop_list, '--o', label=r'$S_n$')
plt.ylabel('MC population (n/cycle)', fontsize = 20)
plt.xlabel('Iteration $i$', fontsize = 20)
plt.xticks(range(len(errors_q3_list)))
plt.grid(which='major',linestyle='-')
plt.grid(which='minor',linestyle='--')
plt.tight_layout()
#plt.legend(fontsize = 15)
fig.savefig(path_res+'PopSizeIt.pdf', format='pdf', dpi=600, bbox_inches='tight')
plt.close()
fig = plt.figure (figsize = (6,4))
plt.plot(range(len(alpha_list)), alpha_list, '--*', label=r'$\alpha$')
plt.ylabel('Under relaxation factor', fontsize = 20)
plt.xticks(range(len(errors_q3_list)))
plt.xlabel('Iteration $i$', fontsize = 20)
plt.grid(which='major',linestyle='-')
plt.grid(which='minor',linestyle='--')
plt.tight_layout()
#plt.legend(fontsize = 15)
fig.savefig(path_res+'AlphaIt.pdf', format='pdf', dpi=600, bbox_inches='tight')
plt.close()
fig = plt.figure (figsize = (6,4))
plt.semilogy(range(len(errors_q3_list)), errors_q3_list, '--x', label=r"$ \frac{\lVert (q'''_i - q'''_{i-1}) \rVert_{L^2}}{ \lVert q'''_i \rVert_{L^2}}$")
plt.axhline(y = tol, color = 'r', linestyle = '--', label="Threshold")
plt.xticks(range(len(errors_q3_list)))
plt.xlabel('Iteration $i$', fontsize = 20)
plt.ylabel(r'Relative $L^2$ Error', fontsize = 20)
plt.grid(which='major',linestyle='-')
plt.grid(which='minor',linestyle='--')
plt.tight_layout()
plt.legend(fontsize = 15)
fig.savefig(path_res+'ErrorsIt.pdf', format='pdf', dpi=600, bbox_inches='tight')
plt.close()
# Plot flux spectrum
fig,ax = plt.subplots(figsize = (10,6))
iCol = np.linspace(0,1,len(phiE_list[0]))
colors = [cm.jet(r) for r in iCol]
for idx in range(len(phiE_list[0])):
plt.stairs(phiE_list[0][idx], EnergyStairs, color = colors[idx], linewidth = 1.5, label = "it "+str(idx))
plt.stairs(phiE_list[0][idx]+1*phiE_list[1][idx], EnergyStairs, baseline=phiE_list[0][idx]-1*phiE_list[1][idx], fill=True, alpha=0.2, facecolor = colors[idx])
plt.xscale('log')
plt.xlabel('Energy (eV)', fontsize = 20)
plt.ylabel(r'Neutron fluence per unit lethargy $\left(\frac{\text{n}}{\text{cm}^2 \, \cdot \, \text{part}} \right)$', fontsize = 16)
plt.grid(linestyle = ':')
plt.tight_layout()
plt.legend(fontsize = 18, loc = 'upper left', framealpha = 1, ncol = 3)
plt.tick_params(axis='both', which='both', labelsize=17)
fig.savefig(path_res+'phiE.pdf', format='pdf', dpi=600, bbox_inches='tight')
plt.close(fig)
Computational wall time = 2.9910 minutes
Computational CPU time = 0.1665 minutes