Voronoi tessellation for spatial domain partitions

For reconstruction using movable sensors, Voronoi tessellation is an essential step that maps the observations to the entire spatial domain. For a given space Ω, which is generally in 2D, a set of points $\left\{r_i\mathrm{,}i=\mathrm{1,...,}{n}_{\text{obs}}\right\}\in \text{Ω}$. The tessellation approach optimally partitions the given space Ω into n_obs regions $G=\left\{g_1\mathrm{,}{g}_2\mathrm{,...,}{g}_{n_{\text{obs}}}\right\}$ using the boundaries determined by the measure d among the given points. Using the measure d, Voronoi tessellation can be expressed as $g_i=\left\{r\in \text{Ω}|d\left(r\mathrm{,}{r}_i\right)<d\left(r\mathrm{,}{r}_j\right)\mathrm{,}j\ne i\right\}\mathrm{.}$ (1) In this study, the Euclidean measure was used, and the Voronoi boundaries between the points were their bisectors. Figure 2 shows an example of 81 points and the related Voronoi partitions in the Euclidean measure. The Voronoi tessellation conveniently projects sparse sensor observations to the global physical field. Thus, learning the map from sparse observations to the global field via CNNs is possible. More importantly, this partitioning process can tolerate local perturbations in the sensor locations. For more details on the mathematical theory of the Voronoi tessellation, we refer readers to [40-43].

Fig. 2Full-screen View

(Color online) Example of 81 points and the related Voronoi tessellation using Euclidean measure

Input and output of the machine learning model

To reconstruct the physical field using machine learning, we prepared the input data for the model using the following process:

(i) Determine the sensor locations $r_{y_i}\mathrm{,}i=\mathrm{1,...,}{n}_{\text{obs}}$ in the reactor core. $r_{y_i}$ may vibrate from its nominal location $r_{y_i}^{\text{nominal}}$. $r_{y_i}=r_{y_i}^{\text{nominal}}+\delta r_{y_i}$, where $\delta r_{y_i}$ presents the amplitude of the vibration. In the following, for reasons of convenience, we denote $\delta r_{y_i}$ as δ.

(ii) Calculate the Voronoi tessellation $s_i\in {ℝ}^{n_x\times n_y}$ using $r_{y_i}\mathrm{,}i=\mathrm{1,...,}{n}_{\text{obs}}$. The Voronoi tessellation first partitions the reactor domain Ω into n_obs regions $G=\left\{g_1\mathrm{,...,}{g}_{n_{\text{obs}}}\right\}$ with $\text{Ω}={\cup }_{i=1}^{n_{\text{obs}}}{g}_i$, such that each region gi contains one sensor located at $r_{y_i}$. si is defined as follows: $s_i\left(r\right)=\{\begin{array}{ll}1\hfill & \text{if}r\in g_i\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\mathrm{.}$ (2)

(iii) Prepare the Voronoi mask field ${\varphi }_{\text{m}}\in {ℝ}^{n_x\times n_y}$ using Voronoi tessellation $s_i\mathrm{,}i=\mathrm{1,...,}{n}_{\text{obs}}$. The element ø_m(r) satisfies ${\varphi }_{\text{m}}\left(r\right)=y_{r_i}\mathrm{,}\text{if}r\in g_i\mathrm{.}$ (3) The Voronoi mask field ϕ_m and related target neutronic field ϕ_n are provided to a machine learning model ${\mathcal{F}}_{\delta }$ such that ${\mathcal{F}}_{\delta }:{\varphi }_{\text{m}}\mapsto {\varphi }_{\text{n}}$, where δ indicates that the model is trained for a given δ. The final output of the model is then denoted by ${\varphi }_{\text{r}}={\mathcal{F}}_{\delta }\left({\varphi }_{\text{m}}\right)$, where the subscript ‘r’ denotes the reconstructed neutronic field. As the specified input vector contains the observed values and position information of the sensors, the proposed model can handle arbitrary sensor locations and an arbitrary number of sensors.

The perturbations of sensor locations and its impact on the effectiveness of the network were also investigated in this study. The typical amplitude δ of the vibration of sensors in a reactor core is less than 1 cm [6], and we investigated the cases of δ= 1 cm, 3 cm, and 5 cm. The effects of the number and locations of sensors in these scenarios are described in the numerical results section.

To construct the training $\left\{{\varphi }_{\text{m}}^{\left(\text{train}\right)}\mathrm{,}{\varphi }_{\text{n}}^{\left(\text{train}\right)}\right\}$ and test $\left\{{\varphi }_{\text{m}}^{\left(\text{test}\right)}\mathrm{,}{\varphi }_{\text{n}}^{\left(\text{test}\right)}\right\}$ sets for the learning process, a physical model of the underlying problem $\mathcal{M}\left(\varphi \left(\mu \right)\mathrm{,}\mu \right)=\mathrm{0,}$ (4) was solved numerically. Here, $\mu \in \mathcal{D}\subset {ℝ}^p$ denotes the p-dimensional parameter of the model, while $\mathcal{D}$ denotes the feasible parameter domain. The training and test sets were accumulated by solving eq:manifold over a discrete set ${\mathcal{D}}^{\left(\text{discrete}\right)}$ which is representative of $\mathcal{D}$.

Learning the map using convolutional neural networks

After preparing the input data for the model, a CNN model can be used to learn the map from observations to the field, in the same manner as handling images [44-46]. In this study, the channel of the CNN was set to one, and for each layer, the extraction of t key features of the input data through filtering operations was expressed as $q_{ijm}^{\left(l\right)}=\sigma \left({\displaystyle \sum _{p=0}^{H-1}{\displaystyle \sum _{c=0}^{H-1}{q}_{i+p-C\mathrm{,}j+c-C\mathrm{,}m}^{\left(l-1\right)}}}{h}_{pcm}+b_m\right)\mathrm{.}$ (5) Here, $C=\text{floor}(H/2)$, $q^{\left(l-1\right)}$ and $q^{\left(l\right)}$ are the input and output data at layer l, respectively; hpcm represents a filter of size H × H and bm is the bias. l_max, H, and m denote the number of layers, filter size, and number of filters, respectively. The output of each filter operation is fed to the activation function $\sigma (\cdot )$ as the output of the neurons. In this study, we selected a rectified linear unit (ReLU) σ(z) = max(0,z) as the activation function [47]. We used the ADAM optimizer with an early stopping criterion for training and threefold cross-validation [48-50]. Furthermore, the L₂ norm was used to measure the errors in the learning process. The detailed parameters of the CNN used in this study are summarized in Table 1.

The Voronoi mask field ${\varphi }_{\text{m}}\in {ℝ}^{n_x\times n_y}$ was used to model $\mathcal{F}$; that is, $q^{\left(0\right)}={\varphi }_{\text{m}}$, and the output of the model was a high-resolution neutronic field $q^{\left(l_{\text{max}}\right)}$. The learning process is formulated as follows. $w=\arg \underset{w}{{\displaystyle \min }}{\Vert q-\mathcal{F}\left({\varphi }_{\text{m}}\mathrm{,}w\right)\Vert }_2\mathrm{,}$ (6) where w denotes the model parameters, specifically the filters of the CNN in this study. After completing the training process, the field was reconstructed by feeding the observations y to the model $\mathcal{F}$, i.e., $\varphi \left(y\right)=\mathcal{F}\left({\varphi }_{\text{m}}\left(y\right)\mathrm{,}w\right)$.

Physical model

In this section, the reconstruction method for nuclear reactor cores is tested. We considered a typical benchmark in nuclear reactor physics, namely the 3D IAEA benchmark problem [51] prepared by the Computational Benchmark Problems Committee of the Mathematics and Computation Division of the American Nuclear Society. This benchmark was selected because it is adapted from realistic reactors, and its geometry and composition are much more complex than those of single-region or two-region problems. After completing this test, the method was directly tested based on real reactor calculations.

To test the algorithm, we used the 2D IAEA case, which represents the midplane z=190 cm of the 3D IAEA benchmark; see [51, page 437] for a detailed description. The 2D geometry of the reactor is shown in Fig. 3, where only one-quarter is shown because the rest can be inferred from the symmetry along the x and y axes. This quarter is denoted by Ω and comprises four subregions with different physical properties: the first three subregions form the core domain Ω_1,2,3, while the fourth subregion is the reflector domain Ω₄. Certain Newman boundary conditions are satisfied on the x and y axes considering symmetry, and the zero-boundary condition is satisfied on the external border, see Fig. 3.

Fig. 3Full-screen View

(Color online) Geometry of the 2D IAEA nuclear core, upper octant: region assignments, lower octant: fuel assembly identification (from [51])

The neutron fields consist of fast and thermal fluxes, that is, $\varphi =\left({\varphi }_1\mathrm{,}{\varphi }_2\right)$ are modeled using a two-group neutron diffusion equation with suitable boundary conditions. The fluxes are the solutions to the following eigenvalue problem (see [52, 53]): Specifically, the flux ϕ satisfies the following eigenvalue problem: Find $\left(\lambda \mathrm{,}\varphi \right)\in ℂ\times L^{\infty }\left(\text{Ω}\right)\times L^{\infty }\left(\text{Ω}\right)$, such that $\{\begin{array}{l}-\nabla \left(D_1\nabla {\varphi }_1\right)+\left({\displaystyle \sum {}_{\text{a}\mathrm{,1}}}+{\displaystyle \sum {}_{1\to 2}}+D_1{B}_{z1}^2\right){\varphi }_1=\frac{1}{\lambda }\nu {\displaystyle \sum {}_{\text{f}\mathrm{,2}}}{\varphi }_2\hfill \\ -\nabla \left(D_2\nabla {\varphi }_2\right)+\left({\displaystyle \sum {}_{\text{a}\mathrm{,2}}}+D_2{B}_{z2}^2\right){\varphi }_2={\displaystyle \sum {}_{1\to 2}}{\varphi }_1\hfill \end{array},$ (7) with the zero boundary condition ϕi = 0, i=1,2 on the external border $\partial \text{Ω}$ and the Newman boundary conditions $\partial \left({\varphi }_i\right)/\partial \left(n\right)=0$, i=1, 2 on axes. The generated nuclear reactor rate $P=\nu {\displaystyle \sum {}_{\text{f},1}}{\varphi }_1$ reflects the power distribution over the core. The following parameters were used in the above equation:

• Di: diffusion coefficient of group i with $i\in \left\{\mathrm{1,2}\right\}$;

• ${\displaystyle \sum {}_{\text{a}\mathrm{,}i}}$: macroscopic absorption cross section of group i;

• ${\displaystyle \sum {}_{\text{1}\to \text{2}}}$: macroscopic scattering cross section from group 1 to 2;

• ${\displaystyle \sum {}_{\text{f,}i}}$: macroscopic fission cross section of group i;

• ν: average number of neutrons emitted per fission.

The axial buckling $B_{zi}^2=0.8\times {10}^{-4}$ for all the regions and energy groups (Table 2). The nominal values of the coefficients in diffusion model (7) are listed in Sect. 2.

Under certain mild conditions of the parameters, the maximum eigenvalue λ_max is real and strictly positive (see [54, Chapter XXI]). The associated eigenfunction ϕ, which is also real and positive at each point x, and is the flux of interest. In neutronics, the inverse of λ_max is conventionally used and is called the multiplication factor $k_{\text{eff}}:=\frac{1}{{\lambda }_{\text{max}}}\mathrm{.}$ (8) For each parameter setting, k_eff is determined by the solution to the eigenvalue problem (7). The maximum eigenvalue λ_max is computed using the well-known power method (see [52]). In this study, Freefem++ [55] was used to solve the 2D IAEA benchmark problem. The spatial approximation used ${\mathbb{P}}_1$ finite elements with a grid size of h = 1 cm; thus, the resolution of the field was nx × ny = 171 × 171.

Field reconstruction

To simulate the variation in the neutronic fields with respect to the parameter variations, we considered the parameters in Sect. 2 as uncertain parameters. Specifically, we assumed the following: $\mu =({\mu }_1\mathrm{,}{\mu }_2\mathrm{,...,}{\mu }_n)\in \mathcal{D}=[{\mu }_{i\mathrm{,}\text{nominal}}\cdot \mathrm{0.8,}{\mu }_{i\mathrm{,}\text{nominal}}\cdot {1.2]}^n$ (9) where $\mathcal{D}$ is the parameter domain. We randomly generated 10000 samplings of μ in $\mathcal{D}$ and solved eq:diffusion to obtain a collection of 10000 samples of neutronic fields $\mathcal{M}=\left\{{\varphi }_1\left(\mu \right)\mathrm{,}{\varphi }_2\left(\mu \right)\mathrm{,}P\left(\mu \right)|\mu \in \mathcal{D}\right\}$. Among then, 8000, 1000, and 1000 samples were used for training, validation, and testing, respectively. We trained three CNN models with the same input data, that is, the thermal flux ø₂ which is measurable with in-core detectors, and the output fields, which are the fast flux ø<₁, thermal flux ø₂, and reaction rate P (see Fig. 4 for a schematic).

Fig. 4Full-screen View

(Color online) Schematic for the reconstruction of neutronic fields

To synthesize the observations, we assumed that an in-core sensor was located at the center of each assembly to acquire the thermal flux. We further assumed that these sensors could move in a local square with a width δ cm, centered at the center of each assembly. We performed numerical tests for the cases δ = 1, 3, and 5 cm to investigate the effect of different levels of vibration of the sensors. Thus, observations were generated randomly from windows of width δ centered at their nominal locations (see Fig. 5). Then, the model $\mathcal{F}$ was trained based on set $\mathcal{M}$, and schematics of the training process are shown in Figs. 1 and 4.

Fig. 5Full-screen View

Different amplitudes of vibrations of a sensor near the nominal location in an assembly

Error metrics

Before presenting the numerical results, we define several metrics to evaluate the quality of the field reconstructions. The normalized root-mean-square residual of the difference between the reconstruction ø_r and test field ø_t is $e_2\left(\varphi \right):=\frac{\parallel {\varphi }_{\text{r}}-{\varphi }_{\text{t}}{\parallel }_2}{\parallel {\varphi }_{\text{t}}{\parallel }_2}\mathrm{.}$ (10) In nuclear engineering, the error in the reconstructed field in L_∞ is another import metric that reflects the worst case scenario for each reconstruction. $e_{\infty }\left(\varphi \right):=\frac{\parallel {\varphi }_{\text{r}}-{\varphi }_{\text{t}}{\parallel }_{L_{\infty }}}{\parallel {\varphi }_{\text{t}}{\parallel }_{L_{\infty }}}\mathrm{.}$ (11) The total average root-mean-square residual and standard deviation over the given test set $\mathcal{M}$ are defined as $\begin{array}{ll}E\left(e_{\chi }\left(\varphi \right)\right)\hfill & :={\text{average}}_{\varphi \in \mathcal{M}}\left(e_{\chi }\left(\varphi \right)\right)\hfill \\ \text{STD}\left(e_{\chi }\left(\varphi \right)\right)\hfill & :={\text{standard deviation}}_{\varphi \in \mathcal{M}}\left(e_{\chi }\left(\varphi \right)\right)\hfill \end{array}\mathrm{,}$ (12) where χ denotes the L₂ or L_∞ norm.

Furthermore, the average assembly field (fluxes and power rate) and related errors were investigated. The average assembly power is defined as ${\varphi }_{\text{ass}\mathrm{,}k}=\frac{1}{v_k}{\displaystyle {\int }_{v_k}}\varphi \text{d}v\mathrm{,}$ (13) where ϕ denotes ϕ₁, ϕ2 or P, vk denotes the volume of the k-th subassembly, and k denotes the fuel assemblies, as shown in the lower octant of Fig. 3. The average error e₂(ϕ_ass), maximum relative error $e_{\infty }\left({\varphi }_{\text{ass}}\right)$, related average $E\left(e_{\chi }\left({\varphi }_{\text{ass}}\right)\right)$, and standard deviation $\text{STD}\left(e_{\chi }\left({\varphi }_{\text{ass}}\right)\right)$ of the errors over the given test set $\mathcal{M}$ can also be similarly defined.

Therefore, we introduce the structural similarity (SSIM) [56] index to measure the field reconstruction. In contrast to the general L₂ error, the SSIM index measures similarity by comparing two images based on luminance, contrast, and structural similarity information.

Performance for the benchmark problem

In Fig. 6, we present the error distributions of the reconstructed fields for different vibration amplitudes, namely, δ = 1,3, and 5 cm, for the 2D IAEA benchmark problem. The reconstruction of ø₁ using observations from ø₂ performs better than those of ø₂ and P. Furthermore, the reconstruction error increases with the amplitude of the vibration, that is, when δ varies from 1 cm to 5 cm. The largest error appears around the interface of the fuel and reflector because of discontinuities in the materials, particularly for the fields of the thermal flux ø₂ and the power rate P.

Fig. 6Full-screen View

(Color online) Reconstructed fields for different vibration amplitudes, δ = 1,3,5 cm for the 2D IAEA benchmark problem

The average assembly values of the reconstructed fluxes and power rate are shown in Fig. 7, and the same conclusions can be drawn. Because of the averaging process, the assembly exhibits considerably smaller relative errors than the pin-wise case.

Fig. 7Full-screen View

(Color online) Reconstructed fields in assembly wise for different vibration amplitudes, and δ = 1,3, and 5 cm for the 2D IAEA benchmark problem

Three main conclusions can be drawn from analysis of the numerical results.

(i) The proposed V-CNN can reconstruct the multi-field using observations only from thermal flux;

(ii) The reconstruction errors for the assembly are significantly lower than 5%, which is acceptable for engineering applications, i.e., less than 10%, which is an acceptable criterion in reactor physics (more information is available in [57]);

(iii) Even with a movement of amplitude δ = 5 cm for the sensors, the proposed V-CNN reconstructs the field with an error less than 10%.

Average performance over a test set

To investigate the generalizability of the field reconstruction method, we analyzed the error performance for 1000 test samples. The error metrics were the average relative L₂ error E(e₂(ϕ)), average relative L_∞ error $E\left(e_{\infty }\left(\varphi \right)\right)$, average relative assembly L₂ error $E\left(e_2\left({\varphi }_{\text{ass}}\right)\right)$, average relative assembly L_∞ error $E\left(e_{\infty }\left({\varphi }_{\text{ass}}\right)\right)$, average SSIM, $E(\text{SSIM}(\varphi ))$, and the standard deviation of the aforementioned errors.

Table 3 illustrates the numerical results of the errors for the reconstruction of ø₁ over the 1000 test samples. All error metrics exhibit good agreement between the reconstructed and original fields. The maximum errors, that is, the L_∞ errors in both the pin-wise and assembly wise cases, are below 2%. The good performance is attributed to the smoothness of the fast flux owing to its relatively longer diffusion length than the thermal flux. Thus, the fast flux is less affected by the heterogeneity of the materials in this benchmark (Fig. 6a for example).

Tables 4 and 5 illustrate the numerical results of the errors for the reconstruction of ø₂ and P over the 1000 test samples. The average L_∞ errors in the pin-wise of the thermal flux and power rate were below 10% when the vibration amplitude was less than 3 cm. When the amplitude of the vibration increases, the average L_∞ error exceeds 10%, which unsuitable for practical engineering applications. However,, the L_∞ errors are much smaller for the assembly. The worst case occurs for δ = 5 cm when reconstructing ø₂, which leads to an error of $E\left(e_{\infty }\left({\varphi }_{\text{ass}}\right)\right)=0.0288$ and standard deviation $\text{STD}\left(e_{\infty }\left({\varphi }_{\text{ass}}\right)\right)=0.0109$. These results further confirm its acceptability for engineering applications. Although the relative L_∞ errors are 10%, most of the points appear around the interface between the fuel and reflector. A relatively large error in this domain is not crucial for safety analysis. Thus, the SSIM indices in all the three tables exceeded 0.99, which further demonstrates the excellent performance for all the field reconstructions.

Robustness analysis

The robustness of the reconstruction with respect to the number of observations n_obs and amount of training data n_snapshot was examined. In Fig. 8, we show the dependence of the relative reconstruction errors in L₂ norm and L_∞ norm on n_snapshot = 128, 1280, 4096, 8192, 15743 and on n_obs=25, 45, 81 to recover the thermal flux ϕ₂ over the test set. The number of observations n_obs=25, 45, and 81 correspond to sparsities of 0.0855%, 0.154%, and 0.277%, respectively, against the number of grid points in the field. The two figures demonstrate the robustness of the proposed method with respect to the sparsity of the sensors and training data. Few observations and training data lead to low-level reconstruction errors. Furthermore, the addition of training data improves the reconstruction accuracy more significantly than the addition of sensors. This result demonstrates that the proposed field reconstruction framework is robust against sensor failures and confirms its potential for practical applications.

Fig. 8Full-screen View

Dependence of the relative reconstruction errors in the L₂ and L_∞ norms on the number of training snapshots n_snapshot = 128, 1280, 4096, 8192, 15743 and the number of observations n_obs=25, 45, 81 for recovering the thermal flux ϕ₂

To investigate the robustness of the recovery with respect to the observation noise, a noise ϵσ that was randomly sampled in the range (-σ,σ) was added to each clean observation y, thereby generating a noisy observation yo=y(1+ϵ_σ) for each sensor. The dependence of the relative reconstruction errors in the L₂ and L_∞ norms of different noise levels, that is, σ=0.01,0.02,0.03,0.04,0.05 for recovering the thermal flux ϕ₂ are shown in Table 6. The test was performed using 81 sensors with a vibration amplitude δ = 5 cm. The errors were first averaged over 100 repeated random observation samplings for each field reconstruction, and then averaged over the test set. On average, the reconstruction error changes significantly for noisy observations. The reconstruction error exhibits a slow linear growth trend with respect to the noise level. Although the L₂ error is below 10%, which is satisfactory for nuclear engineering applications, the L_∞ error remains at approximately 30%, which is unsatisfactory. This provides a direction for further research on reducing the L_∞ error.