简介：资源为基于深度学习的核反应堆中子学模拟项目，面向核工程、计算物理与人工智能交叉方向的毕业设计、课程设计及期末大作业场景。内容聚焦中子扩散方程与中子输运理论，借助神经网络求解有效增殖因子、中子通量分布及多维扩散方程，包含Reactor-Effective-Multiplication-Factor、Multi-Dimensional-Neutron-Diffusion、Differential-Order-Theory-in-Neutron-Transport等模块，共39个文件，以28个Python源码为主，配套5个XML工程配置、3个dat训练/验证数据及README说明，压缩包仅268KB，结构清晰便于直接运行与二次开发。项目覆盖硬边界条件、逆问题、多搜索等典型实验设置，适合用来理解物理信息神经网络在核工程中的应用。目前已有102人学习下载，对需要快速上手深度学习核模拟课题的研究者具有较高参考价值。
1. 核反应堆中子学模拟为什么开始用深度学习
做反应堆物理的人都知道，有效增殖因子 k_eff 和中子通量分布不是“测出来”的，而是“算出来”的。传统蒙特卡罗方法要跟踪几十万条中子历史，算一次临界问题动辄几十分钟甚至数小时；确定论方法虽然快一些，但网格离散和能群划分稍不留神就会引入不可忽略的截断误差。这个压缩包给出的思路是把中子扩散方程当成 PDE 约束，用神经网络去逼近通量场，同时把 k_eff 也作为可学习参数训练出来。对做毕业设计或课程设计的人来说，它最大的价值不是提供了多少个现成结果，而是把“深度学习如何嵌入核工程物理方程”这件事拆成了可运行的脚本：从单能群扩散到多维硬边界条件，再到逆问题和参数扫描，文件命名已经透出清晰的实验路线。
2. 从扩散方程到 PINN：先看懂 MultiDimDiffusionEquation 系列
2.1 方程形式与边界条件的物理含义
反应堆中子学模拟里最常见的是多群扩散方程。以稳态单群为例，方程可以写成：
[ -\nabla \cdot (D \nabla \phi) + \Sigma_a \phi = \frac{1}{k} \nu \Sigma_f \phi ]
其中 (\phi) 是中子通量，(D) 是扩散系数，(\Sigma_a) 是吸收截面，(\nu\Sigma_f) 是裂变产生截面。传统有限差分法把区域离散成网格，逐点求解线性方程组；而 PINN 的做法是让神经网络直接输出 (\phi(x,y,z))，然后把上述方程代入损失函数，用自动微分计算 (\nabla \phi) 和 (\nabla^2 \phi)。这种方式对网格不敏感，也不用构造刚度矩阵，特别适合处理复杂几何边界。
压缩包里的MultiDimDiffusionEquation3_3_1.py到MultiDimDiffusionEquation3_3_6.py这一串文件，数字后缀对应了不同的问题维度、能群数和边界条件组合。从命名规律看，3_3_1、3_3_2 这类文件应该是最基础的验证用例，而后面的 3_3_5、3_3_6 逐渐增加了 hardBC（硬边界）和 MSearch（参数搜索）。我一般会先跑MultiDimDiffusionEquation3_3_1.py，确认环境通，再去看带hardBC的版本，因为硬边界实现方式直接影响收敛速度和最终精度。
2.2 硬边界条件（hardBC）是如何写进网络结构的
所谓硬边界，是指边界条件不是通过损失项“软性”约束，而是直接改网络输出结构，让解在边界上自动满足 Dirichlet 条件。比如对于矩形区域 ([0,1]\times[0,1])，四条边上的通量为零，可以构造一个距离函数 (d(x,y)=x(1-x)y(1-y))，它在边界上恒为 0，内部大于 0。那么网络最终输出就是 (u(x,y)=d(x,y)\cdot N(x,y))，无论 (N) 输出什么，边界值都固定为零。
def forward(self, x, y): # 原始网络输出，形状为 [batch, 1] raw = self.net(torch.cat([x, y], dim=1)) # 构造距离函数，在边界处为 0 d = x * (1 - x) * y * (1 - y) # 乘以距离函数，强制实现零边界条件 return d * raw这段代码的关键在于d的计算方式。x和y是归一化到 ([0,1]) 的坐标，所以x*(1-x)在 x=0 和 x=1 时都是 0，y*(1-y)同理。两者相乘后，四条边全部满足 Dirichlet 条件。如果边界值不是零，可以在后面加上边界常数项，比如return d * raw + bc_value，就变成非零边界。这种做法比在损失函数里加边界项更稳定，因为网络永远不需要去“学习”边界上的值，优化器只需要关注内部残差。对于核反应堆这种对边界通量精度要求高的场景，建议优先采用 hardBC 版本。
2.3 损失函数构造：方程残差、边界残差与 k_eff 的耦合
PINN 的中子学模拟不是简单的回归，它要把物理方程“硬塞”进损失。以稳态扩散方程为例，损失项通常包括三部分：内部方程残差、边界残差、以及针对 k_eff 的额外约束。如果直接把 k_eff 当作可训练参数，那么网络每更新一步，k_eff 也随之变化，训练过程容易震荡。常见做法是先把 k_eff 固定为一个初始猜测值，训练若干轮后再更新。
def compute_loss(model, x, y, D, Sigma_a, nuSigma_f, k_eff): x.requires_grad_(True) y.requires_grad_(True) phi = model(x, y) # 自动微分求一阶和二阶导数 phi_x = torch.autograd.grad(phi, x, grad_outputs=torch.ones_like(phi), create_graph=True)[0] phi_xx = torch.autograd.grad(phi_x, x, grad_outputs=torch.ones_like(phi_x), create_graph=True)[0] phi_y = torch.autograd.grad(phi, y, grad_outputs=torch.ones_like(phi), create_graph=True)[0] phi_yy = torch.autograd.grad(phi_y, y, grad_outputs=torch.ones_like(phi_y), create_graph=True)[0] # 方程残差：-D*(phi_xx+phi_yy) + Sigma_a*phi - (1/k)*nuSigma_f*phi residual = -D * (phi_xx + phi_yy) + Sigma_a * phi - (1.0 / k_eff) * nuSigma_f * phi loss = torch.mean(residual ** 2) return loss这里的create_graph=True是为了保留二阶导的计算图，让 loss 能继续反向传播。k_eff作为标量传入，在训练循环里可以每 N 轮用当前网络解重新估计一次。实际上，当方程本身不含外源项时，裂变源项与吸收项之比就是 k_eff 的迭代值。有些脚本会单独写一个ReactorEffectiveMultiplicationFactor.py，它的思路是在给定通量分布后，用体积分计算产生率与吸收率的比值，把神经网络解当作权函数，这样比单纯训练更稳健。
3. 运行 ReactorEffectiveMultiplicationFactor.py：从文件清单到可复现实验
3.1 先按文件名建立实验地图
拿到压缩包后别急着跑，先花十分钟把脚本归类。文件名里的关键词基本反映了实验变量：hardBC代表硬边界条件，InverseProblem代表反问题，Parallel_MSearch代表并行参数搜索，Single_MSearch代表单参数扫描。把它们之间的关系整理清楚，后续调参就能直接瞄准目标。
| 文件名前缀 | 功能定位 | 典型用途 |
|---|---|---|
MultiDimDiffusionEquation3_3_1/3_3_2 | 基础多维扩散方程求解 | 验证方程离散与网络结构 |
MultiDimDiffusionEquation3_3_3/3_3_4 | 加入不同边界组合 | 对比边界条件对精度的影响 |
MultiDimDiffusionEquation3_3_5/3_3_6 | Single/Parallel 版本 | 单工况或多工况参数扫描 |
*_hardBC.py | 硬边界条件实现 | 提升边界精度，适合毕设展示 |
*_InverseProblem.py | 反问题求解 | 从通量分布反推材料参数 |
ReactorEffectiveMultiplicationFactor.py | k_eff 计算 | 深度学习预测有效增殖因子 |
这个表是我按文件名语义推断的，实际运行时以源码注释为准。但无论如何，_hardBC版本和_InverseProblem版本是值得重点研究的，因为它们不再是一般的 PDE 拟合，而是把核工程里最关心的“算得准”和“反演参数”都覆盖到了。loss.dat、train.dat、test.dat是训练过程的输出文件，分别记录损失下降、训练集指标和测试集误差，是判断收敛的重要依据。
3.2 环境准备与最小运行示例
这类脚本通常基于 PyTorch 或 TensorFlow。如果是 PyTorch，推荐用 Anaconda 创建独立环境，避免和系统 Python 打架。
conda create -n pinn python=3.9 -y conda activate pinn pip install torch numpy matplotlib scipytorch用于神经网络搭建和自动微分，numpy用来做数据读取和数值计算，matplotlib用于绘制通量分布和损失曲线，scipy可以在需要对比传统数值解时用。装完后先跑一个最基础的脚本验证安装：
python MultiDimDiffusionEquation3_3_1.py如果终端里能看到损失从初始值逐步下降，并且在几千步后稳定在一个较低的数值，说明环境没问题。如果报错，优先检查 PyTorch 版本和 CUDA 是否可用。我遇到过最常见的问题是torch.autograd.grad在二阶导时因create_graph=False报错，或者运行目录不对导致找不到loss.dat，所以建议所有脚本都在压缩包解压后的根目录下执行，保持相对路径一致。
3.3 从 train.dat 和 loss.dat 里读训练状态
训练结束后，loss.dat会逐行保存每个 epoch 的损失值。你可以用一段简洁的 Python 脚本读取并绘制曲线：
import numpy as np import matplotlib.pyplot as plt # 读取两列：epoch 和 loss data = np.loadtxt('loss.dat') epoch = data[:, 0] loss = data[:, 1] # 使用对数坐标，因为物理损失常常横跨多个数量级 plt.semilogy(epoch, loss) plt.xlabel('epoch') plt.ylabel('loss') plt.title('PINN training loss') plt.savefig('loss_curve.png', dpi=150)np.loadtxt要求文件列数一致，如果loss.dat里还有额外信息，比如 k_eff 的迭代值，可以调整索引列。对数坐标很关键，因为 PINN 的损失经常从 (10^{-2}) 掉到 (10^{-5})，线性坐标会让人误以为后期没在下降。如果曲线出现平台期，你可以去看train.dat里的预测通量误差，判断是网络容量不够还是方程残差权重太低。
4. 参数调优与 MSearch：并行扫描背后的设计思路
4.1 为什么需要 MSearch 而不是手工调参
深度学习解 PDE 最麻烦的是超参数敏感。学习率、隐藏层数、激活函数、边界权重都会影响最终精度，手工调参一轮可能要跑十几分钟，效率极低。压缩包里带MSearch的文件，实际上是把常用的候选参数组合写成了循环，每个组合独立训练，最后比较 loss 筛选最优解。Parallel_MSearch版本应该是用多进程或 GPU 并行加速这个扫描过程，适合在服务器上批量跑。
import itertools # 候选参数组合 grid = { 'lr': [1e-3, 5e-4, 1e-4], 'n_hidden': [4, 8], 'n_layers': [3, 5] } best_loss = float('inf') best_param = None for lr, hidden, layers in itertools.product(grid['lr'], grid['hidden'], grid['layers']): model = build_pinn(n_hidden=hidden, n_layers=layers) loss = train(model, lr=lr, epochs=5000) if loss < best_loss: best_loss = loss best_param = (lr, hidden, layers) print('best:', best_param, 'loss:', best_loss)itertools.product会穷举所有组合，所以组合数等于各参数数量之积。这里 3×2×2=12 组实验，如果每组 5000 轮，单卡跑完可能要数小时。Parallel_MSearch的核心改进就是把这个循环拆成多个并发任务，比如用 Python 的concurrent.futures.ProcessPoolExecutor分配不同参数到不同 CPU 核或 GPU。对课程设计而言，跑完基本就能在报告里画一张“超参数搜索对比表”，比单纯贴最终结果更有说服力。
4.2 从 loss.dat 判断参数是否合适
训练完成后不要只盯着最终 loss，要看曲线的形状。学习率过大时 loss 会先快速下降，然后震荡甚至发散；学习率过小则下降得异常缓慢，后期可能停留在一个较高的平台。针对核反应堆中子学问题，我一般建议观察前 500 轮的 loss 曲线：如果前 100 轮能下降两个数量级以上，说明网络结构没问题，后面只需要调整边界权重或迭代策略。
如果某个组合的 loss 从一开始就卡在 (10^{-1}) 左右不动，先检查方程残差的计算。常见问题有：坐标归一化没做好，导致距离函数在边界外的值失真；或者自动微分时create_graph没有保留，二阶导计算错误。另一个容易被忽略的点是 k_eff 的初始值。如果初始 k_eff 与实际值相差过大，方程残差的量级会失衡，模型会花很大力气去补偿这个错误。一个实用技巧是先用粗略的有限差分或经验值估算 k_eff 的初值，再用 PINN 精修。
4.3 Single 与 Parallel 版本分别解决什么问题
Single_MSearch应该是在单个 GPU 上串行执行参数扫描，适合代码调试和小规模实验。Parallel_MSearch版本则适合需要同时提交多个工况的场景，比如要对比不同几何尺寸或不同材料截面下的通量分布。如果你在自己的机器上跑Parallel版本遇到 CPU 占用过高或显存溢出，可以把并行数量调小，或者干脆先跑Single版本，把逻辑跑通再上多进程。
在工程实践中，我更喜欢用Parallel_MSearch_hardBC做批量实验，因为硬边界能保证每个候选模型的边界精度一致，最后比较内部残差才有意义。如果拿软边界模型做扫描，边界误差会和内部误差混在一起，选出“最优”参数也不代表方程求解最准。这一点在写毕业设计时可以重点说明，评审老师会觉得你理解了边界处理与超参数扫描之间的关系。
5. 把 PINN 中子学模型改成自己的实验：逆问题与验证技巧
压缩包里最容易被低估的是*_InverseProblem.py。这类脚本解决的是反问题：已知中子通量分布，反推材料参数（如扩散系数或吸收截面）。对于核工程应用，这比正问题更有吸引力，因为反应堆运行期间很难直接测量内部截面，但可以通过探测器获得通量分布。利用 PINN 做反问题的做法并不复杂，预测值和观测值的误差作为损失项，同时保留方程残差作为物理约束。
def inverse_loss(model, x, y, flux_observed, D_param): phi_pred = model(x, y) # 数据拟合损失：与观测通量对比 data_loss = torch.mean((phi_pred - flux_observed) ** 2) # 方程残差损失：把 D_param 当作可训练参数 residual = -D_param * (phi_pred_xx + phi_pred_yy) + Sigma_a * phi_pred - source physics_loss = torch.mean(residual ** 2) return data_loss + physics_loss这里的D_param是设置为requires_grad=True的变量，训练时同时更新网络权重和D_param。如果你想改造成自己的问题，只需要替换flux_observed的来源，比如从实验数据文件读取，或者用传统方法先算一组高精度解作为“伪观测值”。验证时可以用蒙特卡罗或有限差分解对比，计算相对误差。我自己常做的一个验证是：先用MultiDimDiffusionEquation3_3_6_Single_hardBC.py算出通量，把它当作标准答案，再用InverseProblem脚本反推回材料参数，看能不能收敛到预设值。这样既不需要真实实验数据，又能量化反演精度，适合毕设复现。最后一类脚本是Differential-Order-Theory-in-Neutron-Transport/DiffTransport.py，它探讨的是中子输运中的差分阶理论，属于相对进阶的方向，建议在跑通扩散方程后再研究。试验时可以只修改DiffTransport.py顶部的问题尺寸参数，其余逻辑保持默认，这能最快验证理论推导与数值实现是否一致。
本文还有配套的精品资源，点击获取
