ANL11-A2 reactor: multigroup Neutron Diffusion

Aim of this notebook: generate the parametric snapshots for the ANL11-A2 reactor (2D) from the Argonne Code Center - Supplement 2.

This notebook may took a while to run, consider converting it into a .py script.

import numpy as np
import os
from IPython.display import clear_output

from pyforce.tools.functions_list import FunctionsList
from pyforce.tools.backends import LoopProgress
from pyforce.tools.write_read import StoreFunctionsList

import pickle
import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning)

Importing geometry and meshing

The geometry and the mesh can be either created with GMSH or imported from “ANL11A2_octave.msh”, if available.

from neutronics import create_anl11a2_mesh

domain, ct, ft = create_anl11a2_mesh(use_msh=False, save_mesh=True)

fuel1_marker    = 1
fuel2_marker    = 2
fuel_rod_marker = 3
refl_marker     = 4

void_marker     = 10
sym_marker      = 20

clear_output()

Definition of the physical parameters

Since there are different subdomains in \(\Omega\), the values of the parameters changes according to the region. No thermal feedback effects are considered.

The different values describe how the neutrons interact in the different regions of the domain.

regions = [fuel1_marker, fuel2_marker, fuel_rod_marker, refl_marker]

neutronics_param = dict()

neutronics_param['Energy Groups'] = 2

neutronics_param['xs_a'] = [np.array([0.01,  0.01,   0.01,  0.]),
                            np.array([0.085,  0.08,  0.13,  0.01])]
neutronics_param['nu_xs_f'] = [np.array([0., 0., 0., 0.]),
                               np.array([0.135, 0.135,  0.135, 0.])]
neutronics_param['chi'] = [np.array([1,1,1,1]),
                           np.array([0,0,0,0])]
neutronics_param['B2z'] = [np.array([0.8e-4]*4),
                           np.array([0.8e-4]*4)]
neutronics_param['xs_s'] = [[np.array([0]*4), np.array([0.02,  0.02,   0.02,  0.04])],
                            [np.array([0]*4), np.array([0]*4)]]

# The values in the reflector regions are the one of the benchmark
neutronics_param['D'] = [np.array([1.5, 1.5, 1.5, 2.]),
                         np.array([0.4,   0.4,    0.4,   0.3])]

# Normalisation parameters for the flux
nu_value = 2.41 # Number of neutrons per fission
Ef = 200e6 * 1.602e-19 # Energy of fission in Joules
reactor_power = 1e9 # 1 GW

# Void boundary conditions
albedo = [0.4692] * 2

Generation of the snapshots

This section generates the snapshots employing a Finite Element solution. The PDEs are dependent of the parameters \(D_1\) and \(\Sigma_{a,2}\), in the reflector (marker=4) and fuel+rod region (marker=3) respectively:

\begin{equation*} D_1^{\text{refl}} \in [0.45, 3.05] \qquad \qquad \Sigma_{a,2}^{\text{rod}} = [0.0775, 0.1525] \end{equation*}

The variational forms of the neutron diffusion equations are generated to prepare for the numerical solution and the subsequent generation of the snapshots. The class steady_neutron_diff implemented in neutronics.py is used to create the snapshots. The details of the numerical solution are out of the scope of this notebook, interested readers can have a look at the Github repository OFELIA, developed within the ERMETE-Lab.

from neutronics import steady_neutron_diff

# The class is initialized with the mesh, the cell zones, the face tags, the parameters, the region and boundary markers
# An option to used the albedo boundary condition is used
neutr_steady_problem = steady_neutron_diff(domain, ct, ft, neutronics_param, regions, void_marker,
                                           albedo = albedo)

# The variational forms are assembled and the numerical solver is set up to be direct (LU factorization)
neutr_steady_problem.assembleForm(direct=True)

Train Set

The parameters and the FunctionsList used to store the snapshots are defined.

delta_D1 = 0.05
delta_xs_a2 = 0.02

D_1_refl_train  = np.arange(0.5  - 0.1/2,   3.00 + 0.1/2,   delta_D1)
xs_a2_rod_train = np.arange(0.08 - 0.005/2, 0.15 + 0.005/2, delta_xs_a2)

# Number of training snapshots
Ns = len(D_1_refl_train) * len(xs_a2_rod_train)

# Create structures to save the data
train_param = list()
param_idx = 0
phi_1_train = FunctionsList(neutr_steady_problem.V.sub(0).collapse()[0])
phi_2_train = FunctionsList(neutr_steady_problem.V.sub(1).collapse()[0])

For each parameter, the snapshots are generated and stored in the FunctionsList.

progress_bar = LoopProgress(msg='Generating the training snapshots', final = Ns)

for ii in range(len(D_1_refl_train)):
    for jj in range(len(xs_a2_rod_train)):

        # Update the diffusion coefficient in the reflector and fuel1 region (respectively)
        neutr_steady_problem.phys_param['D'][0][3] = D_1_refl_train[ii]
        neutr_steady_problem.phys_param['xs_a'][1][2] = xs_a2_rod_train[jj]

        # Solve the eigenvalue problem
        res = neutr_steady_problem.solve(power = reactor_power, Ef=Ef, nu = nu_value, maxIter = 500, verbose=False)

        # Store the snapshots
        phi_1_train.append(res[0][0].flatten())
        phi_2_train.append(res[0][1].flatten())

        # Store parameters
        train_param.append(dict())
        train_param[param_idx]['D_1_refl']   = D_1_refl_train[ii]
        train_param[param_idx]['xs_a_2_rod'] = xs_a2_rod_train[jj]
        train_param[param_idx]['k_eff'] = res[1]

        # Update bar
        param_idx += 1
        progress_bar.update(1, percentage=False)

Generating the training snapshots: 1.000 / 2.00 - 8.951 s/it

ld: warning: duplicate -rpath '/Users/stefanoriva/miniconda3/envs/ml/lib' ignored
ld: warning: duplicate -rpath '/Users/stefanoriva/miniconda3/envs/ml/lib' ignored

Generating the training snapshots: 2.000 / 2.00 - 10.451 s/it

The neutron flux can reach very high values, which can lead to numerical instabilities. To avoid this, the snapshots are normalized using a min-max normalization, i.e.

\begin{equation*} \tilde{\phi}(\cdot; \boldsymbol{\mu}) = \frac{\phi(\cdot; \boldsymbol{\mu}) - \min\limits_{\boldsymbol{\mu}^*}(\phi(\cdot; \boldsymbol{\mu}^*))}{\max\limits_{\boldsymbol{\mu}^*}(\phi(\cdot; \boldsymbol{\mu}^*)) - \min\limits_{\boldsymbol{\mu}^*}(\phi(\cdot; \boldsymbol{\mu}^*))} \end{equation*}

FunctionsList comes with a method called return_matrix that returns the snapshots in a matrix format, that is as many rows as the spatial mesh and as many columns as the number of snapshots/parameters.

phi_1_train_rescaled = FunctionsList(phi_1_train.fun_space)
phi_2_train_rescaled = FunctionsList(phi_2_train.fun_space)

min_phi_1 = phi_1_train.return_matrix().min()
max_phi_1 = phi_1_train.return_matrix().max()

min_phi_2 = phi_2_train.return_matrix().min()
max_phi_2 = phi_2_train.return_matrix().max()

for mu in range(len(phi_1_train)):
    phi_1_train_rescaled.append((phi_1_train(mu) - min_phi_1) / (max_phi_1 - min_phi_1))
    phi_2_train_rescaled.append((phi_2_train(mu) - min_phi_2) / (max_phi_2 - min_phi_2))

The train snapshots and the correspondent parameters are stored using StoreFunctionsList and pickle, respectively.

path='./Snapshots/'
if not os.path.exists(path):
    os.makedirs(path)

StoreFunctionsList(neutr_steady_problem.domain, phi_1_train_rescaled, 'phi_1', path+'train_snap_phi_1')
StoreFunctionsList(neutr_steady_problem.domain, phi_2_train_rescaled, 'phi_2', path+'train_snap_phi_2')

with open(path+'train.param', 'wb') as f:
    pickle.dump([train_param], f)

Test Set

The parameters and the FunctionsList used to store the snapshots are defined.

D_1_refl_test  = np.arange(0.50, 3.001, delta_D1)
xs_a2_rod_test = np.arange(0.08, 0.151, delta_xs_a2)

# Number of training snapshots
Ns = len(D_1_refl_test) * len(xs_a2_rod_test)

# Create structures to save the data
test_param = list()
param_idx = 0
phi_1_test = FunctionsList(neutr_steady_problem.V.sub(0).collapse()[0])
phi_2_test = FunctionsList(neutr_steady_problem.V.sub(1).collapse()[0])

For each parameter, the snapshots are generated and stored in the FunctionsList.

progress_bar = LoopProgress(msg='Generating the test snapshots', final = Ns)

for ii in range(len(D_1_refl_test)):
    for jj in range(len(xs_a2_rod_test)):

        # Update the diffusion coefficient in the reflector and fuel2 region (respectively)
        neutr_steady_problem.phys_param['D'][0][3] = D_1_refl_test[ii]
        neutr_steady_problem.phys_param['xs_a'][1][2] = xs_a2_rod_test[jj]

        # Solve the eigenvalue problem
        res = neutr_steady_problem.solve(power = reactor_power, Ef=Ef, nu = nu_value, maxIter = 500, verbose=False)

        # Store the snapshots
        phi_1_test.append(res[0][0].flatten())
        phi_2_test.append(res[0][1].flatten())

        # Store parameters
        test_param.append(dict())
        test_param[param_idx]['D_1_refl']   = D_1_refl_test[ii]
        test_param[param_idx]['xs_a_2_rod'] = xs_a2_rod_test[jj]
        test_param[param_idx]['k_eff'] = res[1]

        # Update bar
        param_idx += 1
        progress_bar.update(1, percentage=False)

ld: warning: duplicate -rpath '/Users/stefanoriva/miniconda3/envs/ml/lib' ignored
ld: warning: duplicate -rpath '/Users/stefanoriva/miniconda3/envs/ml/lib' ignored

Generating the test snapshots: 52.000 / 52.00 - 20.585 s/it

Let us rescale also the test snapshots.

phi_1_test_rescaled = FunctionsList(phi_1_test.fun_space)
phi_2_test_rescaled = FunctionsList(phi_2_test.fun_space)

for mu in range(len(phi_1_test)):
    phi_1_test_rescaled.append((phi_1_test(mu) - min_phi_1) / (max_phi_1 - min_phi_1))
    phi_2_test_rescaled.append((phi_2_test(mu) - min_phi_2) / (max_phi_2 - min_phi_2))

The test snapshots and the correspondent parameters are stored using StoreFunctionsList and pickle, respectively.

StoreFunctionsList(neutr_steady_problem.domain, phi_1_test_rescaled, 'phi_1', path+'test_snap_phi_1')
StoreFunctionsList(neutr_steady_problem.domain, phi_2_test_rescaled, 'phi_2', path+'test_snap_phi_2')

with open(path+'test.param', 'wb') as f:
    pickle.dump([test_param], f)