Kickstart

At this moment you have already downloaded the files onto your computer and added the necessary paths (check install) to Matlab path.

In this tutorial you will use the example case within stop/examples/sphere to calculate the criticality radius of a plutonium sphere, and plot the radial flux distribution. So let's get started.

Change directory to sphere folder from within Matlab command window:

>> pwd
/home/user/stop
>> cd examples/sphere
>> ls
data_for_sphere.m  sphere.spi

sphere.spi is a parametrized serpent input file. It contains 2 parameter: TITLE and RAD.

>> dbtype sphere.spi 4:16
4     
5     set title <TITLE>   
6     
7     % --- sph cell:
8     
9     surf 1 sph 0 0 0 <RAD>
10    
11    cell 1 0 plutonium -1
12    cell 2 0 outside    1
13    
14    mat plutonium   -19.84      %g/cm3
15    94239.03c           -1      %Pu-239 wt%
16    

The only thing you need to change in sphere.spi is the path to XS library:

>> dbtype sphere.spi 17:19
17    % --- Cross section library file path:
18    
19    set acelib "/opt/xs/jeff31/sss_jeff31.xsdata"

Now sphere.spi is ready and we start to use Matlab and stop wrapper. Let us test the code. Type in the command window:

>> sim = stop('sphere.spi')

This command uses stop class and creates sim object with the following initial parameters:

sim = 
  stop handle
  Properties:
       serpExe: 'sss'
    saveResPar: {2x1 cell}
          name: []
        values: [1x1 struct]
        isTest: 0
    isContinue: 0
        isEcho: 1
    serpParInp: [1x1 struct]
       saveDir: '/home/user/stop/examples/sphere'

Detailed description of properties of stop wrapper is given in stop-API.

---

WARNING: One needs to specify the name or absolute path of serpent executable if it's not the default sss.

sim.serpExe = '/path/to/serpent/sss' or sim.serpExe = 'yoursssname'

---

To check the sim object just type:

>> sim.serpParInp
ans = 
        file: [1x2719 char]
         dir: '/home/green/stop/examples/sphere'
        name: 'sphere.spi'
    fullname: '/home/green/stop/examples/sphere/sphere.spi'

where serpParInp stands for serpent parametrized input.

Check sphere.spi file content by typing:

>> sim.serpParInp.file

To see the list of parameters in sphere.spi type:

>> sim.displayParameters()
List of unique Parameters in file "/home/green/stop/examples/sphere/sphere.spi"
    'RAD'
    'TITLE'
List of unique equations
    'RAD'
    'TITLE'

To perform calculation with a specific parameters just type:

>> sim.name= 'case'; % simulation name is 'case'
>> sim.run('TITLE', sim.name,'RAD', 3)

Later part of this tutorial is a step by step description of data_for_sphere.m.

To perform multiple simulations with different RAD parameter we use for loop as follows:

>> R = linspace(3,6,20);% create linearly spaced vector between 3 and 6 cm.
>> for i = 1:length(R)
sim.name = ['case_' num2str(i)];
sim.run('TITLE',sim.name,'RAD',R(i));
end

After 2-3 minutes the calculation is over. The code will output 'allresults.mat' file which contains all results.

To load the data into maltab type:

>> clear all; clc
>> load allresults.mat
>> who sphere_case_*

Each of sphere_case_* variables is a data structure and has the following fields:

>> sphere_case_1
sphere_case_1 = 
          serpInp: [1x1 struct]
           stdout: [1x128 char]
    isSimComplete: 1
              res: [1x1 struct]
              det: [1x1 struct]

res and det are data structures which contain data from sphere_case_1_res.m and sphere_case_1_det0.m files. To access their content just type:

>> sphere_case_1.res
ans = 
            ANA_KEFF: [0.6382 0.0029]
                ...
            IMP_KEFF: [0.6391 0.0018]
                ...

To analyse the data and extract k_eff and r(radius of the sphere) from all 20 simulations we use the following code:

name = who('sphere_case_*') % get names of all variables starting with 'sphere_case_'
>> for i = 1:length(name)
data = eval(name{i});
r(i) = data.serpInp.values.RAD;
k_eff(i) = data.res.IMP_KEFF(1,1);  
end

To find the critical radius type:

>> [k_eff' r']
ans =
    0.6391    3.0000
    0.9087    4.4211
    0.9354    4.5789
    0.9633    4.7368
  > 0.9930    4.8947 <
    1.0189    5.0526

To extract radial distribution of neutron flux (actually, neutron flux times volume of spherical layer)

>> [~, idx] = min(abs(k_eff-1.0)); % index of k_eff closest to 1.0
>> data = eval(name{idx});
>> rad_flux = data.det.DET1(:,11);
>> rad_pos = data.det.DET1Z(:,3);
>> plot(rad_pos, rad_flux);