Flattop k-eff sensitivity example
This example demonstrates some basic sensitivity calculation capabilities for the Popsy (Flattop) critical experiment (PU-MET-FAST-006 in the ICSBEP handbook).
Description
The experimental configuration consists of a plutonium alloy sphere (radius 4.5332 cm) reflected by a natural uranium blanket (outer radius 24.1420 cm). The multiplication factor of the system is close to unity (being a critical experiment). This input calculates the sensitivity of the multiplication factor to perturbations applied to the cross sections of 235U, 238U, 239Pu and 240Pu. A processing script for MATLAB/OCTAVE is provided for printing and plotting out the sensitivities.
Input
Input for Popsy (Flattop) sensitivity calculation example:
% --- Core outer surface and blanket outer surface surf s1 sph 0.0 0.0 0.0 4.5332 surf s2 sph 0.0 0.0 0.0 24.1420 % --- Divide space to core, blanket and outside cell c1 0 core -s1 cell c3 0 blanket s1 -s2 cell c4 0 outside s2 % --- Core is made from plutonium alloy mat core sum 94239.03c 3.6697e-2 94240.03c 1.8700e-3 94241.03c 1.1639e-4 31000.03c 1.4755e-3 % --- Blanket is made from natural uranium mat blanket sum 92234.03c 2.6438e-6 92235.03c 3.4610e-4 92238.03c 4.7721e-2 % --- Cross section libraries %set acelib "<path-to-acelib>" % --- Neutron population set pop 10000 100000 50 1 100 100 % --- Use a central point source as the initial source src point sp 0.0 0.0 0.0 % --- Write output each 100 cycles set outp 100 % --- Modify event buffer size to 3 events per neutron set nbuf 3 3 % --- Use unresolved resonance sampling for the fast system set ures 1 2 % --- Define Vitamin-J energy grid (one of the defaults) ene myvit-j 4 nj23 % %%% % %%% --- Sensitivity options % %%% % --- Run sensitivity calculation using Vitamin-J energy grid sens opt egrid myvit-j % --- Use 15 latent generations for the sensitivity calculations sens opt latgen 15 % --- Use 10 generations for Iterated Fission Probability set ifp 10 % %%% % %%% --- Sensitivity responses % %%% % --- Calculate sensitivity of k-effective to perturbations sens resp keff % %%% % %%% --- Sensitivity perturbations % %%% % --- Separate perturbations for the different nuclides and total sens pert zailist 922350 922380 942390 942400 total % --- Do not calculate material-wise perturbations (only total) sens pert matlist total % --- Perturb cross sections separately for each sum reaction mode (fission, capture etc.) sens pert xs all % --- Score sensitivity results directly to a matrix for each particle % Assume 0.2 matrices needed per generation per particle sens opt direct 0.2
Output
The contents of some long arrays have been replaced with "...".
Output for Popsy (Flattop) sensitivity calculation example:
% Number of different bins in sensitivity calculation: SENS_N_MAT = 1; SENS_N_ZAI = 5; SENS_N_PERT = 7; SENS_N_ENE = 175; % Materials included in sensitivity calculation: SENS_MAT_LIST = [ 'total ' ]; % Indices for different materials: iSENS_MAT_TOT = 1; % Nuclides included in sensitivity calculation: SENS_ZAI_LIST = [ 0 % total 922350 922380 942390 942400 ]; % Indices for different ZAIs: iSENS_ZAI_TOT = 1; iSENS_ZAI_922350 = 2; iSENS_ZAI_922380 = 3; iSENS_ZAI_942390 = 4; iSENS_ZAI_942400 = 5; % Reactions included in sensitivity calculation: SENS_PERT_LIST = [ 'total xs ' 'ela scatt xs ' 'sab scatt xs ' 'inl scatt xs ' 'capture xs ' 'fission xs ' 'nxn xs ' ]; % Indices for different perturbations: iSENS_PERT_TOT_XS = 1; iSENS_PERT_ELA_XS = 2; iSENS_PERT_SAB_XS = 3; iSENS_PERT_INL_XS = 4; iSENS_PERT_CAPT_XS = 5; iSENS_PERT_FISS_XS = 6; iSENS_PERT_NXN_XS = 7; % Sensitivity calculation energy group boundaries: SENS_E = [ ... ]; % Sensitivity calculation energy group lethargy widths: SENS_LETHARGY_WIDTHS = [ ... ]; % Sensitivities with 15 latent generations: % Effective multiplication factor: ADJ_PERT_KEFF_SENS = [ ... ]; ADJ_PERT_KEFF_SENS = reshape(ADJ_PERT_KEFF_SENS, [2, SENS_N_ENE, SENS_N_PERT, SENS_N_ZAI, SENS_N_MAT]); ADJ_PERT_KEFF_SENS = permute(ADJ_PERT_KEFF_SENS, [5, 4, 3, 2, 1]); ADJ_PERT_KEFF_SENS_E_INT = [ ... ]; ADJ_PERT_KEFF_SENS_E_INT = reshape(ADJ_PERT_KEFF_SENS_E_INT, [2, SENS_N_PERT, SENS_N_ZAI, SENS_N_MAT]); ADJ_PERT_KEFF_SENS_E_INT = permute(ADJ_PERT_KEFF_SENS_E_INT, [4, 3, 2, 1]);
Data-processing and plotting
Script for printing and plotting the results tested with MATLAB R2014a and OCTAVE 4.0.0:
%% --- Flag to print out zero values
printzeros = false;
%% --- Input file name
fname = 'input';
%% --- Run the sensitivity output file
run([fname '_sens.m']);
%% --- Print some basic output data
disp(sprintf('Results from Serpent sensitivity calculation using:\n'));
%% --- Materials
disp([num2str(SENS_N_MAT) ' perturbed materials:']);
for i=1:1:SENS_N_MAT
disp(SENS_MAT_LIST(i,:));
end
%% --- Nuclides
disp(sprintf(['\n' num2str(SENS_N_ZAI) ' perturbed nuclides (ZAIs):']));
for i=1:1:SENS_N_ZAI
disp(sprintf('%d',SENS_ZAI_LIST(i,:)));
end
%% --- Perturbations
disp(sprintf(['\n' num2str(SENS_N_PERT) ' different perturbations:']));
for i=1:1:SENS_N_PERT
disp(SENS_PERT_LIST(i,:));
end
%% --- Energy grid
disp(sprintf(['\nEnergy grid with ' num2str(SENS_N_ENE) ' bins.']));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% -- Output sensitivities one by one %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(sprintf('\n********************\n\nSensitivity results:\n'))
txt = sprintf('%-20s %-6s %-20s : %12s | 2 sigma confidence interval', 'Material', 'ZAI', ' Perturbation', 'sensitivity');
disp(txt)
%% --- This script only prints/plots the k-eff sensitivities
%% You can change the following three lines to plot something else
res = ADJ_PERT_KEFF_SENS;
resE = ADJ_PERT_KEFF_SENS_E_INT;
resStr = 'K-eff';
for iMat=1:1:SENS_N_MAT
matStr = SENS_MAT_LIST(iMat,:);
for iZai=1:1:SENS_N_ZAI
myZai = SENS_ZAI_LIST(iZai);
%% --- Turn ZAI into a string for printing
if (myZai == 0)
zaiStr = sprintf('%-6s ', 'total');
else
zaiStr = sprintf('%6d ', SENS_ZAI_LIST(iZai));
end
%% --- Last loop is over perturbations
for iPert=1:1:SENS_N_PERT
%% --- Turn perturbation into a string for printing
pertStr = sprintf('%-20s ', SENS_PERT_LIST(iPert,:));
val = resE(iMat,iZai,iPert,1);
err = resE(iMat,iZai,iPert,2);
%% --- Calculate 2 sigma confidence interval assuming
%% normally distributed random variable
minV = min([val*(1-2*err), val*(1+2*err)]);
maxV = max([val*(1-2*err), val*(1+2*err)]);
%% --- Create a formatted string for printing out the sensitivity
txt = sprintf('%-20s %-6s %-20s: %+12.8f | [%+12.8f, %+12.8f]', matStr, zaiStr, pertStr, val, minV, maxV);
%% --- Print formatted string
if (printzeros || (val ~= 0) || (err ~= 0))
disp(txt);
%% --- Move on to plotting a stairstep plot
%% --- Convert energy grid to eV
xvals = SENS_E*1e6;
%% --- Get energy dependent sensitivity and relative error
yvals = squeeze(res(iMat, iZai, iPert, :, 1))./SENS_LETHARGY_WIDTHS';
yerr = squeeze(res(iMat, iZai, iPert, :, 2));
%% --- We'll need to re-append the last value if we want to
%% get the stairstep plot correct
yvals = [yvals; yvals(end)];
yerr = [yerr; yerr(end)];
%% --- Create the figure for plotting
fig1 = figure('visible', 'off');
%% --- Plot the sample mean
h1 = stairs(xvals, yvals, 'k-');
set(h1, 'LineWidth', 2);
%% --- Plot the sample mean +- 2*sqrt(sample variance)
hold on;
h2 = stairs(xvals, yvals.*(1+2*yerr), 'r-');
h3 = stairs(xvals, yvals.*(1-2*yerr), 'r-');
%% --- Fill the area of non-zero-sensitivity
N = size(xvals,2)*2-2;
xx = zeros(N,1);
yy = zeros(N,1);
xx(1:2:N) = xvals(1:1:end-1);
xx(2:2:N) = xvals(2:1:end);
yy(1:2:N) = yvals(1:1:end-1);
yy(2:2:N) = yvals(1:1:end-1);
h4 = fill(xx, yy, [1.0 0.8 0.8]);
hold off;
%% --- Make the picture a bit nicer
set(gca,'xscale','log');
set(gca,'XTick', [1e4 1e5 1e6 1e7]);
grid on;
box on;
xlim([1e4, SENS_E(end)*1e6]);
ylim([-0.1, 0.5]);
%% --- Label the axes and add a title
xlabel('Energy (eV)');
ylabel('Sensitivity per unit lethargy');
txt = sprintf('Response: %s\nMaterial: %s\nZAI: %s\n perturbation %s', resStr, matStr, zaiStr, pertStr);
title(txt);
%% --- Save the plot to file
filename = sprintf('%s_sens_%ld_%ld_%ld.png', resStr, iMat, iZai, iPert);
print(filename,'-dpng');
%% --- Close figure
close(fig1);
end
end
%% --- Add empty line after each ZAI
disp(sprintf(' '));
end
end
Results
Runtime of the input on a computer cluster node with a Intel Xeon E5-2690 v2 processor (10 physical cores, 20 logical cores) using 20 OpenMP threads took approximately 4 hours 50 minutes and used (allocated) 6117 MB of RAM. The run produced an 133 KB [input]_sens.m file.
Executing the processing and plotting script produced the following output using OCTAVE 4.0.0
Output from the processing script:
Results from Serpent sensitivity calculation using: 1 perturbed materials: total 5 perturbed nuclides (ZAIs): 0 922350 922380 942390 942400 7 different perturbations: total xs ela scatt xs sab scatt xs inl scatt xs capture xs fission xs nxn xs Energy grid with 175 bins. ******************** Sensitivity results: Material ZAI Perturbation : sensitivity | 2 sigma confidence interval total total total xs : +0.90527000 | [ +0.90411125, +0.90642875] total total ela scatt xs : +0.16188300 | [ +0.16084695, +0.16291905] total total inl scatt xs : +0.07863120 | [ +0.07815941, +0.07910299] total total capture xs : -0.05447680 | [ -0.05459665, -0.05435695] total total fission xs : +0.71923200 | [ +0.71897308, +0.71949092] total total nxn xs : +0.00125871 | [ +0.00122598, +0.00129144] total 922350 total xs : +0.00818980 | [ +0.00809971, +0.00827989] total 922350 ela scatt xs : +0.00081898 | [ +0.00074855, +0.00088941] total 922350 inl scatt xs : +0.00038832 | [ +0.00035803, +0.00041861] total 922350 capture xs : -0.00054016 | [ -0.00054913, -0.00053120] total 922350 fission xs : +0.00752267 | [ +0.00747753, +0.00756781] total 922350 nxn xs : +0.00000661 | [ +0.00000463, +0.00000859] total 922380 total xs : +0.21901900 | [ +0.21805532, +0.21998268] total 922380 ela scatt xs : +0.13649800 | [ +0.13562441, +0.13737159] total 922380 inl scatt xs : +0.06499990 | [ +0.06459690, +0.06540290] total 922380 capture xs : -0.03974930 | [ -0.03986060, -0.03963800] total 922380 fission xs : +0.05727050 | [ +0.05714450, +0.05739650] total 922380 nxn xs : +0.00096484 | [ +0.00093783, +0.00099186] total 942390 total xs : +0.65488800 | [ +0.65424621, +0.65552979] total 942390 ela scatt xs : +0.02212830 | [ +0.02159722, +0.02265938] total 942390 inl scatt xs : +0.01241460 | [ +0.01216631, +0.01266289] total 942390 capture xs : -0.01294550 | [ -0.01299469, -0.01289631] total 942390 fission xs : +0.63329000 | [ +0.63301135, +0.63356865] total 942390 nxn xs : +0.00027191 | [ +0.00025614, +0.00028768] total 942400 total xs : +0.02006720 | [ +0.01991469, +0.02021971] total 942400 ela scatt xs : +0.00129135 | [ +0.00117513, +0.00140757] total 942400 inl scatt xs : +0.00063176 | [ +0.00056606, +0.00069746] total 942400 capture xs : -0.00102246 | [ -0.00103534, -0.00100958] total 942400 fission xs : +0.01916650 | [ +0.01909367, +0.01923933] total 942400 nxn xs : +0.00000945 | [ +0.00000567, +0.00001323]
These results were obtained using ENDF/B-VII.1 based cross section libraries.
A number of produced figures are also shown here for ease of comparison:
