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 ²³⁵U, ²³⁸U, ²³⁹Pu and ²⁴⁰Pu. 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

% --- 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 MT number

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:

%% --- 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