View source for 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 <sup>235</sup>U, <sup>238</sup>U, <sup>239</sup>Pu and <sup>240</sup>Pu. A processing script for MATLAB/OCTAVE is provided for printing and plotting out the sensitivities.

== Input ==

<div class="toccolours mw-collapsible mw-collapsed" style="width:60em;">
'''Input for Popsy (Flattop) sensitivity calculation example:'''
<div class="mw-collapsible-content">
 <nowiki>% --- 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</nowiki>
</div>
</div>

== Output ==

The contents of some long arrays have been replaced with "...".

<div class="toccolours mw-collapsible mw-collapsed" style="width:70em;">
'''Output for Popsy (Flattop) sensitivity calculation example:'''
<div class="mw-collapsible-content">
 <nowiki>
% 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]);</nowiki>
</div>
</div>

== Data-processing and plotting ==

<div class="toccolours mw-collapsible mw-collapsed" style="width:80em;">
'''Script for printing and plotting the results tested with MATLAB R2014a and OCTAVE 4.0.0:'''
<div class="mw-collapsible-content">
 <nowiki>%% --- 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</nowiki>
</div>
</div>

== Results ==