Weak-Form Neural Networks Solve Multigroup Neutron Diffusion Equations Without Spatial Derivatives

Computational modeling of nuclear reactor cores relies fundamentally on multigroup neutron diffusion equations to evaluate core power distributions, assess operational margins, and confirm critical safety parameters. In nuclear reactor analysis, steady-state eigenvalue equations determine the effective multiplication factor and group-wise neutron flux profiles across complex core geometries, while transient formulations track temporal power evolution following control rod movements or localized material property changes. Standard engineering practices have historically depended on mature numerical discretization frameworks, such as the finite-difference method, finite-element formulations, and nodal diffusion schemes. Although these established numerical approaches provide robust engineering solutions, increasing geometrical dimensionality, intricate core partitioning, and multi-physics coupling impose substantial administrative and computational workloads when constructing meshes, transferring cross-program data, and managing repeated solution workflows. Consequently, computational physicists have sought flexible, differentiable continuous-field representations capable of unifying forward simulations, parameter identification, and multi-scale data assimilation within modern scientific computing pipelines.

Physics-informed neural networks have emerged as a promising methodology for continuous partial differential equation solutions, embedding governing equations, initial states, and boundary constraints directly into optimization loss functions. However, when conventional strong-form physics-informed neural networks are applied to heterogeneous multigroup neutron diffusion problems, significant computational bottlenecks emerge. Standard strong-form formulations require the evaluation of second-order spatial derivatives of neural network outputs via automatic differentiation at collocation points distributed across the computational domain. In reactor cores characterized by piecewise-discontinuous diffusion coefficients and sharp material interfaces, classical solution regularity deteriorates significantly. The resulting disparities in numerical loss scales and conflicting gradient vectors among distinct residual terms frequently induce optimization stiffness, slow training convergence, or convergence to suboptimal local solutions. Furthermore, three-dimensional extensions of strong-form networks introduce prohibitive computational complexity, creating an unfavorable trade-off between the density of spatial collocation points and overall training duration.

To overcome the operational limitations and derivative burdens inherent in strong-form neural approximations, a research team from China Nuclear Power Engineering Co., Ltd. developed a novel physics-informed neural network framework based on elementwise weak-form discrete physical residuals. Rather than evaluating high-order spatial derivatives through point-by-point automatic differentiation, the method leverages standard finite-element shape functions to construct continuous flux interpolation fields and assemble discrete physical operator matrices prior to neural training. The neural network predicts nodal flux values directly, while preassembled sparse matrices evaluate diffusion-loss, intergroup scattering, and fission source operators through algebraic matrix-vector multiplications. By shifting derivative calculations entirely into preassembled finite-element operators, the training loop relies on automatic differentiation strictly for computing the gradients of the physical loss with respect to network parameters, thereby eliminating spatial derivative evaluations of network outputs during iterative optimization.

The mathematical formulation of the steady-state two-group neutron diffusion eigenvalue problem models fast and thermal neutron populations while neglecting up-scattering phenomena. In this system, all fission neutrons are assumed to be born within the fast energy group, with fission events induced primarily by thermal-group neutrons. The physical system is governed by group-dependent diffusion coefficients, removal cross sections, intergroup scattering cross sections, and thermal fission-production cross sections. To formulate the weak discrete physical residual, each term in the multigroup diffusion equation is multiplied by an arbitrary test function that vanishes along homogeneous Dirichlet boundaries. Integrating by parts over the spatial domain via Green’s identity converts the second-order diffusion derivatives into inner products of the flux and test-function gradients. Because test functions vanish on zero-flux outer boundaries and normal neutron currents are zero on reflective symmetry planes, the corresponding surface boundary integrals vanish, naturally satisfying normal current balance across material interfaces without explicit penalty functions.

To manage the mutual coupling between the effective multiplication factor and the spatial neutron flux distribution, the framework incorporates an outer-inner neural source iteration algorithm. In reactor eigenvalue problems, the effective multiplication factor governs the magnitude of the fission source in the fast-group balance equation, while the fast- and thermal-group fluxes jointly determine the spatial distribution of fission and slowing-down sources. In each outer iteration, normalized two-group fluxes from the preceding step are held fixed to establish the source terms for the fast and thermal equations. An inner optimization stage subsequently executes fifty network parameter updates using weak-form discrete loss evaluations to approximate the corresponding fixed-source problem. Upon obtaining updated nodal flux vectors via forward evaluation, the effective multiplication factor is recalculated through a weighted least-squares minimization of the discrete fast-group balance residual, followed by under-relaxation using a factor of 0.2 to prevent inter-iteration numerical oscillations without jointly optimizing the eigenvalue alongside network weights.

Geometrical representation and boundary enforcement within the network are achieved through specialized trial-function architectures and problem-specific coordinate embeddings. Both two-dimensional and three-dimensional configurations employ fully connected networks featuring three hidden layers with 64 neurons per layer, utilizing hyperbolic tangent activation functions and a final Softplus mapping to guarantee flux non-negativity. For the two-dimensional TWIGL benchmark, network inputs consist of normalized Cartesian coordinates augmented with continuous sine and cosine harmonic features, multiplied by an analytical cosine envelope that enforces zero flux along vacuum outer boundaries and zero normal derivatives along reflective symmetry axes. For the three-dimensional IAEA pressurized water reactor benchmark, the model uses a 32-dimensional feature vector containing normalized coordinates, harmonic sine-cosine functions, four Gaussian features centered at control rod positions, and three hyperbolic tangent layer functions aligned with axial reflector and control rod interfaces, complemented by an axial sine envelope to enforce boundary conditions.

The performance of the weak-form discrete neural framework was progressively verified on standard reactor benchmarks, beginning with the two-dimensional TWIGL quarter-core problem. Discretized using linear P1 triangular elements on a 40 by 40 physical training mesh, the network was trained for 16,000 optimization steps using the AdamW optimizer with cosine annealing, reaching a best physical residual loss of 1.73 × 10−5. In the steady-state evaluation, the neural network yielded an effective multiplication factor of 0.9132141, differing by only 0.04 pcm from the reference finite-element method result of 0.9132137. The root-mean-square errors for both the fast- and thermal-group fluxes remained below 2.4 × 10−4 relative to the reference finite-element distribution, while the maximum absolute nodal flux errors across the entire heterogeneous quarter-core domain remained strictly below 2.0 × 10−3, demonstrating high local and global solution fidelity across disparate material compositions.

Beyond steady-state eigenvalue solutions, the researchers validated the method’s practical utility within a hybrid steady-state-transient workflow by transferring the trained neural field to time-dependent reactor kinetics. The steady-state nodal fluxes obtained on the TWIGL mesh were projected onto an 81 by 81 node grid to initialize a conventional Crank-Nicolson finite-element transient calculation with a single delayed-neutron precursor group. The hybrid simulation investigated both step and ramp perturbations involving reductions in the thermal absorption cross section of the core’s primary material over a 0.5-second transient window with a time step of 0.001 seconds. In both perturbation scenarios, the transient power curves computed from the neural network initialization closely tracked pure finite-element reference calculations, exhibiting root-mean-square errors in relative core power of 5.27 × 10−4 for the step transient and 4.57 × 10−4 for the ramp transient, with maximum absolute deviations from published POLCA-T benchmark data remaining below 0.1 percent.

The computational scalability of the unified weak-form discrete framework was further challenged using the three-dimensional IAEA pressurized water reactor benchmark, which includes an irregular radial core boundary, axial upper and lower reflectors, and partially inserted control rod banks. The three-dimensional domain was discretized into trilinear Q1 hexahedral elements on a uniform 2-centimeter mesh, incorporating one-eighth radial symmetry folding and diagonal symmetry expansion weights. The network was optimized over 24 staged training intervals comprising 12,000 total steps using the L-BFGS optimization algorithm, reducing the normalized physical loss from an initial random state of 5.8477 down to 5.3353 × 10−5. The final converged effective multiplication factor reached 1.0289669, exhibiting an absolute deviation of only 3.67 pcm when compared against a same-mesh Q1 finite-element power iteration calculation that converged to 1.0290036.

Spatial accuracy assessments across the three-dimensional core revealed close agreement between the neural prediction and conventional numerical solutions across both neutron energy groups. The root-mean-square errors of the neural fast- and thermal-group flux distributions relative to the same-mesh finite-element baseline were 3.008 × 10−3 and 1.487 × 10−3, respectively, with corresponding maximum absolute errors of 1.121 × 10−2 and 2.486 × 10−2. Downstream assembly-integrated fission source calculations demonstrated that the physics-informed spatial feature embeddings successfully captured local flux depressions induced by the partially inserted control rods and flux reflections near outer reflector zones. Evaluated against published Argonne National Laboratory benchmark reference data, the neural network achieved a core-wide assembly relative power root-mean-square error of 0.001939, with a mean absolute relative error of 0.1779 percent and a maximum absolute relative error of 0.4647 percent occurring in peripheral fuel assemblies.

Computational timing comparisons performed on a standard personal laptop using a single central processing unit core highlighted the trade-offs between initial neural network training duration and standard numerical execution. For the steady-state two-dimensional TWIGL problem, training the neural network required 110.96 seconds, whereas the direct steady-state finite-element solution required only 0.665 seconds. In the three-dimensional IAEA core simulation, the 12,000-step staged neural training procedure required 9,833.16 seconds, or approximately 2.73 hours of single-core compute time. During transient operations, executing the hybrid Crank-Nicolson time advancement initialized by the neural state required 64.56 seconds for the step perturbation and 64.76 seconds for the ramp perturbation, compared to 51.52 and 53.46 seconds for the baseline finite-element simulations, confirming that transient time-stepping performance remains comparable once the neural initialization is established.

The study authors noted that while current training times for isolated single-condition forward simulations remain longer than those of optimized linear solvers, the distinctive engineering value of the weak-form discrete neural framework lies in its continuous, differentiable field representation and seamless modularity with conventional numerical solvers. By establishing physical operators through finite elements, enforcing natural interface continuity without penalty parameters, and employing neural source iteration for eigenvalue decoupling, the method provides a unified foundation across diverse spatial dimensions and element geometries. The authors anticipate that this differentiable framework will prove advantageous in multi-query applications, such as core parameter optimization sweeps, inverse cross-section estimation, data assimilation, and transfer learning across shifting operating conditions, where the initial offline training cost of the neural field can be amortized over numerous rapid evaluations.

Subject of Research: Weak-form discrete physical-residual neural network methods for multigroup neutron diffusion eigenvalue and transient problems

Article Title: A weak-form discrete physical-residual PINN method for solving multigroup neutron diffusion equations

Article References: Hu, K., Wang, X., Zhang, P., Zhang, Y., Feng, Z., Xiao, X., Xu, Y., & Yuan, Y. (2026). A weak-form discrete physical-residual PINN method for solving multigroup neutron diffusion equations. Frontiers in Energy Research, 14, 2005399. https://doi.org/10.3389/fenrg.2026.2005399

Image Credits: AI Generated

DOI: 10.3389/fenrg.2026.2005399

Keywords: neutron diffusion, physics-informed neural networks, weak-form discretization, finite element method, nuclear reactor physics, TWIGL benchmark, IAEA-3D PWR benchmark, effective multiplication factor, multigroup diffusion, neural source iteration, Crank-Nicolson scheme, eigenvalue problems, computational physics, reactor core simulation, spatial discretization, delayed neutrons