ANL11-A2 reactor: multigroup Neutron Diffusion

This notebooks implements the generation of the snapshots arising from 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.

[1]:

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

import gmsh
from dolfinx.io.gmshio import model_to_mesh, read_from_msh

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

import pickle

from neutronics import steady_neutron_diff

Mesh and parameters

The geometry and the main physical parameters will be assigned.

Importing geometry and meshing

The geometry and the mesh are imported from “ANL11A2_octave.msh”, generated with GMSH.

[2]:

gdim = 2

model_rank = 0
mesh_comm = MPI.COMM_WORLD

use_msh = False

# If the snapshots are generated from a different computer than the one used for the offline/online simulation, there may be issues in the number of dofs
if (use_msh):
    domain, ct, ft = read_from_msh("ANL11A2_octave.msh", comm = mesh_comm, rank = model_rank, gdim = gdim)
else:
    # Initialize the gmsh module
    gmsh.initialize()

    # Load the .geo file
    gmsh.merge('ANL11A2_octave.geo')
    gmsh.model.geo.synchronize()

    # Set algorithm (adaptive = 1, Frontal-Delaunay = 6)
    gmsh.option.setNumber("Mesh.Algorithm", 6)

    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.optimize("Netgen")

    # Domain
    domain, ct, ft = model_to_mesh(gmsh.model, comm = mesh_comm, rank = model_rank, gdim = gdim )
    gmsh.finalize()

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

void_marker     = 10
sym_marker      = 20

tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)

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.

[ ]:

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 = 1
Ef = 1
reactor_power = 1

# Void boundary conditions
albedo = [0.4692] * 2

Generation of the snapshots

This section generates the snapshots from the Finite Element solver, dependent on 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.5, 3] \qquad \qquad \Sigma_{a,2}^{\text{rod}} = [0.08, 0.15] \end{equation*}

The variational forms are generated to prepare for the numerical solution and the subsequent generation of the snapshots.

[4]:

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

neutr_steady_problem = steady_neutron_diff(domain, ct, ft, neutronics_param, regions, void_marker,
                                           albedo = albedo)
neutr_steady_problem.assembleForm()

Train Set

The train set of snapshots is computed and stored.

[8]:

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

# 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])

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: 2.000 / 2.00 - 13.351 s/it

The train snapshots and the correspondent parameters are stored using pyforce and pickle.

[ ]:

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

StoreFunctionsList(neutr_steady_problem.domain, phi_1_train, 'phi_1', path+'train_snap_phi_1')
StoreFunctionsList(neutr_steady_problem.domain, phi_2_train, 'phi_2', path+'train_snap_phi_2')

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

Test Set

The test set of snapshots is computed and stored.

[ ]:

D_1_refl_test  = np.arange(0.50, 3.001, 0.1)
xs_a2_rod_test = np.arange(0.08, 0.151, 0.05)

# 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])

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)

The test snapshots and the correspondent parameters are stored using pyforce and pickle.

[ ]:

StoreFunctionsList(neutr_steady_problem.domain, phi_1_test, 'phi_1', path+'test_snap_phi_1')
StoreFunctionsList(neutr_steady_problem.domain, phi_2_test, 'phi_2', path+'test_snap_phi_2')

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