Geometry module

This module is primarily used for constructing the geometric structure and material composition of models. It utilizes the Constructive Solid Geometry (CSG) method for geometric descriptions and employs half-space surface methods to build complex computational models, filling the corresponding geometric regions with specified material nuclide compositions and proportions. It also supports ultrafine millimeter-level voxel modeling[13]. Due to the tedious and error-prone nature of manual geometry model creation in MC simulations, this module enables the use of CAD software to construct the geometric input for MagicMC, ensuring consistency between CAD models and MagicMC models. The STEP standard neutral file is used to convert CAD models into MagicMC models[14]. The STEP standard provides a feasible solution for data exchange and sharing between different CAD systems through system-independent neutral files. The specific implementation method is shown in Fig. 2. Firstly, we establish a set of standard modeling processes based on advanced 3D visual modeling platforms. Second, we conduct research on the conversion from Boundary Representation (BREP) to Constructive Solid Geometry (CSG), constructing three-dimensional solid natural half-spaces and generating characteristic points. Finally, based on the CSG expressions and combined with the surface information of three-dimensional solids, we convert this data into input files for the MagicMC platform.

Fig. 2Full-screen View

CAD-based conversion method for the MagicMC model

Source module

This module is designed to establish the initial conditions and properties of particles in general or critical sampling sources, involving time distribution, spatial distribution, angular distribution, and energy distribution. The spatial distribution types are divided into point and volume sources (e.g., rectangular, cylindrical, spherical, and spherical shell sources). Point sources specify particle positions in space using three coordinates (X, Y, and Z). To describe and sample volume sources, spatial location information is required. For example, to uniformly sample source particles in a parallelepiped, the user must specify the diagonal coordinates of the body in the input card to define the spatial distribution of the source particles. The particle angular distribution is primarily isotropic or anisotropic and is defined by three parameters: the polar angle $(-1\le \cos \theta \le 1$), azimuthal angle $(0\le \varphi \le 2\pi )$, and reference vector $(u\mathrm{,}v\mathrm{,}w)$. Each specific distribution, defined by a probability density function (PDF), has its own sampling method. The energy distribution primarily includes monoenergetic, Watt, Maxwell, discrete, and histogram distributions, where the energy range of incident neutrons is from 0 to 20 MeV, and the energy range of incident photons is from 1 keV to 100 GeV. For a monoenergetic distribution, users set a single fixed value for particle energy (MeV), such as setting the photon energy of ¹³⁷Cs to 0.662 MeV. The Watt distribution indicates that particle energy follows the Watt fission spectrum, described by the distribution function $p(E)\text{d}E=c{e}^{-E/a}\sin {\left(bE\right)}^{1/2}\text{d}E$ with parameters $a$, $b$. The Maxwell distribution is represented by $p(E)\text{d}E=c{E}^{1/2}{e}^{-E/a}\text{d}E$, which requires setting the parameter $a$. The discrete distribution specifies source energies at discrete points, each with an associated sampling probability. The histogram distribution involves setting each energy bin and its corresponding probability value, with the energy boundary points monotonically increasing. During sampling, the energy bins are first sampled, followed by uniform sampling within the selected energy bin. This method, based on a discrete description of the probability distribution in each energy bin, allows for a more flexible simulation of actual distribution scenarios by dividing the energy range into different bins. Additionally, particle types, such as photons or neutrons, can be set, as shown in Fig. 3.

Fig. 3Full-screen View

Parameters definition for source particles

Particle transport module

The particle transport module is the core module of the basic unit, linking upwards to the geometry module, source module, and database module, and supporting the tally and output modules below. The transport module is independently encapsulated and connected through interfaces. This module primarily handles neutron and photon transport calculations. The main particle transport tracking in this section is similar to that of mainstream Monte Carlo programs (Fig. 4)[15-17], and the specific flow is as follows: (i) Perform source particle sampling to determine the initial state of the particle $S_m=\left(r_m\mathrm{,}{E}_m\mathrm{,}{\text{Ω}}_m\right)$. The particle departs with this initial state, and the transport length is sampled according to the macroscopic cross-section ${\displaystyle \sum {}_{\text{t}}}\left(E_m\right)$ of the material at the current particle collision point to obtain $L=-\ln \zeta /{\displaystyle \sum {}_{\text{t}}}\left(E_m\right)$. (ii) Calculate the distance d from the particle to the nearest geometric boundary along its direction of motion. (iii) Determine the relationship between the transport length L and the nearest boundary distance d. The position of the (m+1)-th collision is given by $r_{m+1}=r_m+L_m{\text{Ω}}_m$. (iv) To determine the energy and direction of the particle after a collision, it is necessary to sample the particle interactions to identify the type of reaction occurring with a specific nuclide in the current material. First, using the macroscopic cross-sections and total cross-section ${\displaystyle \sum {}_{\text{1}}}\mathrm{,}{\displaystyle \sum {}_{\text{2}}}\mathrm{,}\dots \mathrm{,}{\displaystyle \sum {}_n}\mathrm{,}{\displaystyle \sum {}_{\text{t}}}$ of n nuclides in the current material, a random number ζ uniformly distributed in [0,1] is sampled to identify the interacting nuclide. For example, $0\le \zeta \le {\displaystyle \sum {}_{\text{1}}}/{\displaystyle \sum {}_{\text{t}}}$ indicates that the particle interacts with the first nuclide. After sampling the colliding nuclide, the microscopic cross-section σ is used to determine the specific type of interaction between the particle and the current nuclide. The main neutron interactions include elastic scattering, inelastic scattering, and absorption, where absorption is divided into (n, α), (n, p), (n, γ), and (n, f), etc. The total microscopic cross-section of a neutron is ${\sigma }_{\text{t}}={\sigma }_{\text{e}}+{\sigma }_{\text{in}}+{\sigma }_{\text{a}}$, and the type of neutron-matter interaction is determined based on the random number ζ in the range of $[\mathrm{0,}{\sigma }_{\text{e}}/{\sigma }_{\text{t}}]$, $[{\sigma }_{\text{e}}/{\sigma }_{\text{t}}\mathrm{,}{\sigma }_{\text{e}}+{\sigma }_{\text{in}}/{\sigma }_{\text{t}}]$, and $[{\sigma }_{\text{e}}+{\sigma }_{\text{in}}/{\sigma }_{\text{t}}\mathrm{,1}]$ for neutron elastic scattering, inelastic scattering, and absorption, respectively, with the current nuclide. Photon-matter interaction types include the photoelectric effect, electron-pair production, and Compton scattering, with similar determination methods. (v) After determining the colliding nuclide and type of reaction, the energy and direction of the particle $f(E_m\mathrm{,}{\text{Ω}}_m\to E_{m+1}\mathrm{,}{\text{Ω}}_{m+1})$ can be established.

Fig. 4Full-screen View

Particle transport tracking process

(1) Dynamic Monte Carlo Particle Transport Method

During nuclear reactor operations, transient processes such as power regulation and rod ejection accidents cause changes in geometry models, which usually cannot be simulated with conventional Monte Carlo codes. Based on the MagicMC Monte Carlo simulation framework, the dynamic behaviors of moving geometric components in the reactor can be modeled in real-time by updating the spatial positions of dynamic geometric components according to the particle time-state parameters. The specific implementation involves updating the dynamic geometric lattice elements, corresponding to vector VVV and rotation matrix RRR, based on the neutron time variables and dynamic geometric motion parameters (translation direction, translation speed, rotation direction, rotation speed, etc.). In MagicMC, particle tracking is divided into two categories: tracking in dynamic geometries and tracking in fixed geometries. The tracking process in a fixed geometry is similar to that in standard Monte Carlo codes. For tracking in a dynamic geometry, the positions of geometric cells are first updated based on the motion parameters of the dynamic components. The particle tracking process then proceeds similarly to that of standard Monte Carlo codes. This process is illustrated in Fig. 5.

Fig. 5Full-screen View

Particle tracking method in MagicMC

(2) CPU-GPU Heterogeneous Parallel Transport Method

With the development of new nuclear reactors, high-fidelity and high-resolution modeling methods have become increasingly important. For reactor transients and multiphysics problems, the computational requirements for achieving satisfactory solutions using the Monte Carlo method are quite expensive, and the central processing unit (CPU) computing resources on personal computers have gradually become insufficient. Unlike standard Monte Carlo codes with history-based transport algorithms, MagicMC employs a stack-driven event-based transport algorithm to enable GPU parallelization. The event-based particle transport algorithm divides the entire transport computation into various types of events, including source sampling, cross-sectional calculations, particle forward distance, collision, surface penetration, particle termination, and particle revival events. Based on the type of the next event each particle must process, multiple active particles are distributed into corresponding stacks. Different stack events are executed in order; the stack length of the event is determined before each event starts, and the longest stack event is selected for calculation until the stack is empty. In subsequent iterations, the longest stack event is selected for computation until the stack is empty. At this point, new particle events are emitted into the event stack, and the entire process is repeated until all event stacks are empty, thereby ending the program. By parallel processing particles undergoing the same events, thread divergence is minimized, and GPU memory utilization is maximized. This allows more particles to be tracked simultaneously, significantly reducing computational time, as shown in Fig. 6 [18, 19].

Fig. 6Full-screen View

CPU-GPU heterogeneous parallel Monte Carlo particle transport method

Database module

This module stores nuclear data cross-sections essential for reactor physics and radiation shielding design calculations, primarily based on the CENDL database, and is also compatible with other databases like JENDL and ENDF. It includes microscopic cross-sections for nuclear reactions induced by neutrons across various energy ranges with various nuclides, angular distributions of neutrons for elastic and inelastic scattering, yields, and energy spectra for prompt and delayed fission neutrons, and thermal neutron scattering $S(\alpha \mathrm{,}\beta )$ data for moderators. It also contains cross-sectional data for photon interactions, including the photoelectric effect, Compton effect, and pair production. For activation analysis, the activation database includes cross-sectional data, uncertainty data for neutron-induced reactions, decay data, fission product yield data from neutron-induced reactions, and $\gamma$ absorption data. For calculations such as reactor depletion analysis and spent fuel composition analysis, the module also includes the corresponding depletion chain library, detailed isotopic transformations, and decay channels over time. For medical dosimetry calculations in BNCT, a library of 106 human tissue compositions from ICRU reports and the corresponding Kerma factor library were created and stored, with functionalities to modify and edit the related data.

Tally module

The tallies module primarily accounts for each particle's contribution to the expected solution during its motion, calculating and solving flux density using collision estimation, track length estimation, and absorption estimation methods. Specifically, MagicMC solves for the dynamic parameters and adjoint flux in critical systems by approximating the iterated fission probability using the number of fission neutrons from the last generation of neutrons. By considering first-generation neutrons as source neutrons and tallying all fission neutrons produced by the last generation in the tallies module, the neutron worth of the source neutrons can be obtained. Eventually, by dividing the neutron worth of the source neutrons by the number of source neutrons in the corresponding region and energy group, the iterated fission probability—approximately equal to the adjoint flux—can be obtained. After obtaining the adjoint flux for each region and energy group, the effective delayed neutron fraction ${\beta }_{\text{eff}}$ and the neutron generation time $\text{Λ}$ are calculated. ${\beta }_{\text{eff}}=\frac{{\displaystyle {\sum }_i{\displaystyle {\sum }_m{\displaystyle \int }{\varphi }^{*}}}\left(r\mathrm{,}E\mathrm{,}\text{Ω}\right){\chi }_{\text{d}i}^m\left(r\mathrm{,}E\right)u_{\text{d}i}^m\left(r\mathrm{,}{E}^{\prime }\right)\cdot {\displaystyle {\sum }_{\text{f}}^m\left(r\mathrm{,}{E}^{\prime }\right)}\varphi \left(r\mathrm{,}{E}^{\prime }\mathrm{,}{\text{Ω}}^{\prime }\right)\text{d}{\text{Ω}}^{\prime }\text{d}{E}^{\prime }\text{d}\text{Ω}\text{d}E\text{d}r}{{\displaystyle {\sum }_m{\displaystyle \int }{\varphi }^{*}\left(r\mathrm{,}E\mathrm{,}\text{Ω}\right)}{\chi }_{\text{t}}^m\left(r\mathrm{,}E\right)u_{\text{t}}^m\left(r\mathrm{,}{E}^{\prime }\right)\cdot {\displaystyle {\sum }_{\text{f}}^m\left(r\mathrm{,}{E}^{\prime }\right)}\varphi \left(r\mathrm{,}{E}^{\prime }\mathrm{,}{\text{Ω}}^{\prime }\right)\text{d}{\text{Ω}}^{\prime }\text{d}{E}^{\prime }\text{d}\text{Ω}\text{d}E\text{d}r}$ (2) $\text{Λ}=\frac{{\displaystyle \int }{\varphi }^{*}\left(r\mathrm{,}E\mathrm{,}\text{Ω}\right)\varphi \left(r\mathrm{,}E\mathrm{,}\text{Ω}\right)/v\left(r\mathrm{,}E\right)\text{d}\text{Ω}\text{d}E\text{d}r}{{\displaystyle {\sum }_m{\displaystyle \int }{\varphi }^{*}}\left(r\mathrm{,}E\mathrm{,}\text{Ω}\right){\chi }_{\text{t}}^m\left(r\mathrm{,}E\right)u_{\text{t}}^m\left(r\mathrm{,}{E}^{\prime }\right){\displaystyle {\sum }_{\text{f}}^m\left(r\mathrm{,}{E}^{\prime }\right)}\cdot \varphi \left(r\mathrm{,}{E}^{\prime }\mathrm{,}{\text{Ω}}^{\prime }\right)\text{d}{\text{Ω}}^{\prime }\text{d}{E}^{\prime }\text{d}\text{Ω}\text{d}E\text{d}r}$ (3) Ω, $E$, and $r$ represent the angle, energy, and position, respectively; the superscript $m$ denotes fissile nuclides; $i$ represents the group of delayed neutrons; subscripts d and t represent delayed neutrons and all neutrons, respectively; ${\text{Φ}}^{*}$ is the adjoint flux; $\chi$ is the fission spectrum; $u$ is the number of neutrons produced per fission; ${\displaystyle \sum {}_{\text{f}}^m}$ represents the macroscopic fission cross-section; $\text{Φ}$ is the neutron flux; and $\nu$ is the neutron speed. The process flow for calculating the kinetic parameters using MagicMC is shown in Fig. 7.

Fig. 7Full-screen View

Process flow of the kinetic parameter tally method

Output module

This module primarily displays the results of data tallied by the tallies module, such as particle flux, reaction rate, dose rate, power distribution, effective multiplication factor, adjoint flux, and kinetic parameters ($\alpha$, ${\beta }_{\text{eff}}$, $\text{Λ}$, $\rho$). To visually present results in areas of interest to users in a more intuitive way and facilitate the analysis and optimization of geometric structure, material layout, and radiation shielding effects, the MagicMC employs visualization tools from the auxiliary module. This approach addresses the challenge of visualizing multiple radiation fields simultaneously, as the different radiation field data resolutions (mesh sizes) vary by location. Multiple VTK data fields have been developed for these visualization methods [20], as shown in Fig. 8. The visualization processing flow is as follows:

Fig. 8Full-screen View

Data visualization method with multiple VTK format fields

(1) Based on MagicMC, physical field calculations are performed for different tally regions, followed by mesh tallying to obtain multiple fields with mesh data (at various resolutions).

(2) The physical field mesh tally results for each region are processed by extracting spatial topology and field distribution data, forming a standardized VTK data format.

Auxiliary module

The auxiliary module is designed to enhance the main functionalities and improve the maintainability, user-friendliness, and overall performance of the code bench. It includes libraries related to the functions of the aforementioned modules, such as mathematical libraries for various probability samplings in source sampling (math), file operation functions, and other functional utilities related to code execution. To enhance stability, this module contains an exception-catching mechanism (Fatal Error, Warning) capable of detecting and handling exceptional situations during the operation of MagicMC, thereby increasing the fault diagnosis capability of the code-bench. Visualization tools are mainly used to display the data results passed from the output module in a 2D/3D format, assisting users in more intuitively understanding the output of the code bench. The interface information output conveys information to the user, providing the running status, progress reports, date, and the number of running threads, facilitating user interaction with the code bench.

Variance reduction functionality

In nuclear reactor physics calculations, especially in large-scale radiation shielding transport calculations, the large-scale and extremely complex structure of reactor geometries can cause the neutron flux rate from the core to the exterior to vary significantly–sometimes by more than ten orders of magnitude. When using the Monte Carlo method, this often leads to issues such as low-probability deep penetration problems, high computational resource demands, and convergence difficulties. MagicMC has developed a mesh-based weight-window module to address large-scale deep-penetration challenges [21]. The global variance reduction (GVR) technique, based on the MC anti-forward transport method, was derived from the weight window lower limit bias parameter generation method of MC anti-forward transport. Unlike other variance reduction methods (e.g., CADIS), this technique does not require consideration of code coupling and mesh mapping. Figure 9 shows the computational flow of the global variance reduction technique based on MC anti-forward transport in MagicMC [22].

Fig. 9Full-screen View

Global variance reduction method based on anti-forward flux

The adaptive global variance reduction method is illustrated in Fig. 10. Based on the model data provided by the user, a particle transport simulation is conducted to obtain the forward flux. Using the simulation results, data processing is performed with Formula (4), introducing a smoothing factor SI ($SI\in [\mathrm{0,}1]$) to generate importance weight windows. The smoothing factor SI is used to prevent excessive particle splitting due to large gradients in the lower parameter of the weight window. Typically, the SI value defaults to 1.0. However, when the majority of the simulated geometry is shielded, SI can be set to 0.5–1.0, and when only the far-source region of the simulated geometry is shielded, the SI value can be set between 0.0–0.5. These values for SI were obtained through empirical test calculations. The generated weight window is then used to update the input file for the next generation, and the number of particles is doubled from the previous generation. Following this concept, each generation is iterated in the same manner. Users can select the number of iterations and the initial particle count for the first generation. This functionality also allows for continued computation if system errors occur or if initial particle counts yield insufficient results, ensuring the accuracy of the calculations. $w_{\text{th}}\left(r\right)={\left(\frac{\varphi \left(r\right)}{\varphi \left(\max \right)}\right)}^{SI}$ (4) Here, $w_{\text{th}}\left(r\right)$ represents the lower limit of the weight window in space $r$, $\text{Φ}(r)$ is the flux rate in space $r$, ${\text{Φ}}_{\text{max}}$ represents the maximum flux for the entire model, and $SI$ is the smoothing factor with variable values.

Fig. 10Full-screen View

Adaptive global variance reduction method based on forward flux

High‑resolution activation functionality

During nuclear reactor operation, neutron activation reactions lead to the production of a large number of radioactive nuclides within the reactor's structural materials. This directly impacts radiation safety considerations, including reactor shielding design, maintenance planning, refueling strategies, and decommissioning plans—crucial aspects of reactor design [23-26]. Conventional activation methods, typically cell-based, do not consider the spatial distribution of neutron flux within a cell and thus cannot accurately describe the nonuniform spatial effects of activation in structural materials. This limitation can lead to coarse distributions of source terms and decay $\gamma$-source, thereby limiting the precise evaluation of the radiation dose induced by decay $\gamma$-source. To address these issues, this study developed a mesh-based high-resolution activation module for reactor structural materials using the Rigorous-2-Step (R2S) method, as shown in Fig. 11. This process includes three main steps: neutron transport, mesh-based activation, and decay photon dose calculations.

Fig. 11Full-screen View

Mesh-based activation analysis method based on R2S

(1) Neutron Transport Calculations: A Monte Carlo model of reactor structural materials is established for neutron transport calculations. Mesh tally models are created for structural material components. Based on MagicMC, Monte Carlo neutron transport calculations are performed to obtain mesh-based neutron flux and reaction rates for structural material regions.

(2) Mesh-based Activation Calculation: For structural material components, a mesh model is created for activation calculations. Mesh activation calculations are performed using mesh-based neutron flux and reaction rates from neutron transport calculations. This provides a high-resolution distribution of mesh-based activation source terms, including mesh decay photon source strength, decay photon energy spectrum, and the main activated nuclides.

(3) Decay Photon Dose Calculation: Using the mesh-based decay photon source derived from mesh-based activation calculations, a mesh-based source sampling method is integrated into MagicMC. A Monte Carlo photon transport calculation is then performed with the mesh-based decay photon source, and the radiation dose rates induced by the decay photon source are obtained through mesh tallying.

Shielding intelligent optimization functionality

Reactor radiation shielding design is a typical multi-objective optimization problem, requiring comprehensive consideration of factors such as weight, volume, and radiation dose rates to achieve an optimal design solution [27-36]. Traditional design methods heavily rely on expert experience and are limited by low efficiency, high uncertainty, and suboptimal solutions, making them unsuitable for rapid design processes. The shielding optimization module in MagicMC includes built-in algorithm interfaces, enabling the automated invocation of multi-objective optimization algorithms, such as NSGA-II, NSGA-III, and MOPSO. This facilitates the iterative optimization of shielding schemes, thereby optimizing the design of shielding structures and materials. The main steps of the shielding intelligent optimization module are as follows:

(1) Pre-processing: Acquiring user-specified structures or regions for optimization, defining optimization objectives (e.g., weight, volume, and dose), constraint parameters, etc.

(2) Optimization Simulation: Utilizing evolutionary algorithms for automated radiation shielding encoding and selecting appropriate algorithms for optimization design. Through the algorithm interface, data transmission between the MagicMC transport calculation and optimization algorithms is completed, iteratively obtaining the optimal solution set.

(3) Post-processing: Recommending design solutions based on evaluation indicators and visually presenting the optimization results.

The shielding optimization module leverages the advantages of various optimization algorithms, effectively balancing multiple design objectives and rapidly determining the optimal parameters for reactor shielding design, including shielding thickness, structural layout, and material composition. Within the defined constraints, it achieves optimization objectives related to minimizing shielding structure volume, space occupation, and protection costs. This module provides valuable guidance for novel reactor-shielding designs. The shielding optimization module is illustrated in Fig. 12.

Fig. 12Full-screen View

Shielding optimization functionality with evolutionary algorithms

Multiphysics coupling functionality

The normal operation of a nuclear reactor involves the coupling of multiple physical fields, spanning disciplines such as neutronics, heat transfer, fluid dynamics, and structural mechanics [37-39]. In transient-oriented neutronics and thermal-hydraulics coupling problems, MagicMC can provide time-dependent calculations of reactor power and its distribution changes. Using neutronics data, users can perform transient thermal calculations to determine the temperature field of a reactor at different times. MagicMC generates temperature-dependent cross-sections online through a mapping relationship. This process calculates the impact of thermal feedback on reactor operation, achieving transient neutronics and thermal-hydraulics coupling. The transient neutronics and thermal-hydraulics coupling based on MagicMC is illustrated in Fig. 13.

Fig. 13Full-screen View

Transient neutronics and thermal-hydraulics coupling method

Monte Carlo codes require substantial computational resources for time-dependent calculations, particularly for long-duration transient processes. This demand, coupled with the need to perform calculations across multiple physical fields, makes the simulation computation time a critical limiting factor. MagicMC includes a module for calculating the average reactor temperature coefficient as part of neutronics-thermal-mechanical coupling calculations for reactors, building on online temperature-dependent cross-sectional processing. This module uses the point kinetics method to evaluate changes in reactor power and power distribution. Integration with ANSYS for finite-element thermal stress calculations yields results such as reactor temperature fields, structural displacements, and stresses. These results are then fed back to both MagicMC and the point kinetics code, completing the neutronics-thermal-mechanical coupling calculations for the reactor [40]. The process is illustrated in Fig. 14.

Fig. 14Full-screen View

Multiphysics coupling of the neutronics-thermal-mechanical coupling method