器學(xué)習(xí)的中子學(xué)PINN研究:物理信息神經(jīng)網(wǎng)絡(luò)求解中子輸運(yùn)方程實(shí)戰(zhàn))
簡介這份資源面向核工程、人工智能方向的畢業(yè)設(shè)計(jì)、課程設(shè)計(jì)與期末大作業(yè)學(xué)習(xí)者聚焦物理信息神經(jīng)網(wǎng)絡(luò)PINN在中子學(xué)中的應(yīng)用。內(nèi)容圍繞有效增殖因子計(jì)算、多維中子擴(kuò)散方程的無網(wǎng)格求解以及中子輸運(yùn)微分階理論展開將物理定律約束嵌入神經(jīng)網(wǎng)絡(luò)訓(xùn)練幫助讀者理解如何用機(jī)器學(xué)習(xí)方法替代傳統(tǒng)復(fù)雜數(shù)值求解流程。壓縮包共39個(gè)文件約269KB以28個(gè)Python腳本為核心覆蓋單硬邊界、并行多搜索、反問題等多種實(shí)驗(yàn)配置另含5個(gè)xml與1個(gè)iml等開發(fā)環(huán)境配置、3個(gè)dat訓(xùn)練與損失數(shù)據(jù)文件及README說明目錄按反應(yīng)堆參數(shù)、擴(kuò)散方程、輸運(yùn)理論等模塊組織便于按主題檢索與復(fù)現(xiàn)。目前已有46人學(xué)習(xí)下載適合希望快速搭建PINN中子學(xué)實(shí)驗(yàn)框架、對(duì)照代碼理解物理約束建模與排錯(cuò)思路的讀者參考。1. 從“算不動(dòng)”到“算得快”這份中子學(xué) PINN 資源到底解決了什么問題如果你做過中子學(xué)相關(guān)課題大概率經(jīng)歷過這種絕望輸運(yùn)方程維度一高傳統(tǒng)數(shù)值方法網(wǎng)格量直接爆炸算一個(gè)工況等半天調(diào)個(gè)參數(shù)又要重來。這份《基于機(jī)器學(xué)習(xí)的中子學(xué)PINN研究.zip》就是沖著這個(gè)痛點(diǎn)來的——它把物理信息神經(jīng)網(wǎng)絡(luò)PINN搬進(jìn)中子學(xué)場景用神經(jīng)網(wǎng)絡(luò)逼近中子通量密度分布把物理方程作為約束塞進(jìn)損失函數(shù)讓模型在少量數(shù)據(jù)甚至無標(biāo)簽條件下也能收斂到物理上合理的解。資源里包含完整的代碼實(shí)現(xiàn)、訓(xùn)練腳本和實(shí)驗(yàn)配置適合做畢業(yè)設(shè)計(jì)、課程設(shè)計(jì)或期末大作業(yè)的同學(xué)直接上手復(fù)現(xiàn)也適合想從傳統(tǒng)數(shù)值方法轉(zhuǎn)向機(jī)器學(xué)習(xí)求解偏微分方程的從業(yè)者做技術(shù)驗(yàn)證。它不是純理論推導(dǎo)而是能跑起來、能改參數(shù)、能看到損失曲線下降的工程包。2. 拆開壓縮包之前PINN 求解中子輸運(yùn)的底層邏輯與選型理由2.1 為什么是中子學(xué) PINN而不是純數(shù)據(jù)驅(qū)動(dòng)中子輸運(yùn)方程本質(zhì)上是玻爾茲曼方程在中子場下的具體形式變量包括空間坐標(biāo)、能量、方向和時(shí)間。傳統(tǒng)離散縱標(biāo)法或蒙特卡洛方法精度高但計(jì)算成本隨維度指數(shù)上升。純數(shù)據(jù)驅(qū)動(dòng)模型需要大量標(biāo)注數(shù)據(jù)而中子學(xué)實(shí)驗(yàn)數(shù)據(jù)獲取成本極高很多工況根本拿不到實(shí)測(cè)值。PINN 的核心思路是把方程殘差作為損失的一部分讓網(wǎng)絡(luò)在訓(xùn)練中同時(shí)擬合邊界條件和物理規(guī)律這樣即使沒有標(biāo)簽數(shù)據(jù)也能通過最小化方程殘差得到近似解。常見做法是把中子通量密度作為網(wǎng)絡(luò)輸出輸入是空間點(diǎn)和能量點(diǎn)損失函數(shù)由三部分組成方程殘差項(xiàng)、邊界條件項(xiàng)、初始條件項(xiàng)。這種混合驅(qū)動(dòng)方式在數(shù)據(jù)稀缺場景下優(yōu)勢(shì)明顯也是這份資源選擇 PINN 而不是普通神經(jīng)網(wǎng)絡(luò)的原因。2.2 資源里的網(wǎng)絡(luò)結(jié)構(gòu)與訓(xùn)練流程拆解從代碼結(jié)構(gòu)看資源采用全連接網(wǎng)絡(luò)作為主干激活函數(shù)用 tanh因?yàn)樗亩A導(dǎo)數(shù)連續(xù)適合計(jì)算二階微分方程殘差。輸入維度根據(jù)具體問題設(shè)定一般包括空間坐標(biāo) x、y、z 和能量群 e輸出是中子通量。訓(xùn)練流程分三步先采樣配置點(diǎn)包括內(nèi)部殘差點(diǎn)和邊界點(diǎn)然后前向傳播計(jì)算通量預(yù)測(cè)值最后用自動(dòng)微分求方程殘差反向傳播更新權(quán)重。下面是一個(gè)典型的訓(xùn)練循環(huán)骨架我按資源里的寫法整理出來你可以直接對(duì)照自己的環(huán)境改import torch import torch.nn as nn class PINN(nn.Module): def __init__(self, layers): super().__init__() self.net nn.Sequential() for i in range(len(layers)-1): self.net.add_module(flinear_{i}, nn.Linear(layers[i], layers[i1])) if i len(layers)-2: self.net.add_module(ftanh_{i}, nn.Tanh()) def forward(self, x): return self.net(x) # 損失函數(shù)方程殘差 邊界條件 def loss_fn(model, x_int, x_bc, D, Sigma_a, S): x_int.requires_grad_(True) phi model(x_int) # 自動(dòng)微分求二階導(dǎo)數(shù) phi_x torch.autograd.grad(phi, x_int, torch.ones_like(phi), create_graphTrue)[0] phi_xx torch.autograd.grad(phi_x, x_int, torch.ones_like(phi_x), create_graphTrue)[0] # 擴(kuò)散方程殘差: -D*phi_xx Sigma_a*phi - S 0 residual -D * phi_xx Sigma_a * phi - S loss_res torch.mean(residual ** 2) # 邊界條件通量在邊界處為零或給定值 phi_bc model(x_bc) loss_bc torch.mean(phi_bc ** 2) return loss_res loss_bc這段代碼里D是擴(kuò)散系數(shù)Sigma_a是宏觀吸收截面S是源項(xiàng)這三個(gè)參數(shù)需要根據(jù)你實(shí)際的中子學(xué)問題從截面庫或文獻(xiàn)里取值。x_int是內(nèi)部采樣點(diǎn)x_bc是邊界點(diǎn)采樣數(shù)量直接影響訓(xùn)練穩(wěn)定性和精度。資源里默認(rèn)內(nèi)部點(diǎn) 2000 個(gè)、邊界點(diǎn) 200 個(gè)你可以根據(jù)問題維度調(diào)整。注意create_graphTrue必須開否則二階導(dǎo)數(shù)求不出來這是新手最容易翻車的地方。2.3 環(huán)境配置與依賴版本資源基于 PyTorch 框架建議用 Python 3.8 以上PyTorch 1.10 以上。CUDA 版本根據(jù)你顯卡驅(qū)動(dòng)選沒有 GPU 也能跑但訓(xùn)練時(shí)間會(huì)從幾分鐘拉到半小時(shí)以上。依賴包包括 numpy、matplotlib、scipy用來做數(shù)據(jù)預(yù)處理和結(jié)果可視化。安裝命令如下pip install torch torchvision numpy matplotlib scipy如果你用 conda可以建一個(gè)獨(dú)立環(huán)境避免版本沖突conda create -n pinn_neutron python3.9 conda activate pinn_neutron pip install torch numpy matplotlib scipy裝完之后先跑一個(gè)官方 MNIST 例子驗(yàn)證 PyTorch 能正常調(diào)用 GPU再跑資源里的訓(xùn)練腳本這樣能排除環(huán)境問題。我一般會(huì)先torch.cuda.is_available()確認(rèn)一下返回 False 就老老實(shí)實(shí)用 CPU 跑小規(guī)模測(cè)試。3. 把代碼跑起來從數(shù)據(jù)生成到損失收斂的完整操作鏈3.1 采樣策略內(nèi)部點(diǎn)、邊界點(diǎn)、初始點(diǎn)怎么分配PINN 訓(xùn)練效果很大程度上取決于采樣點(diǎn)分布。資源里采用均勻采樣加拉丁超立方采樣兩種方式均勻采樣適合規(guī)則區(qū)域拉丁超立方適合高維或復(fù)雜幾何。內(nèi)部點(diǎn)數(shù)量建議是網(wǎng)絡(luò)參數(shù)量的 5 到 10 倍太少會(huì)導(dǎo)致方程殘差約束不足太多會(huì)拖慢訓(xùn)練速度。邊界點(diǎn)要覆蓋所有邊界每個(gè)邊界至少 50 個(gè)點(diǎn)否則邊界條件容易失效。初始點(diǎn)只在瞬態(tài)問題里需要穩(wěn)態(tài)問題可以忽略。下面是一個(gè)采樣函數(shù)示例import numpy as np from scipy.stats import qmc def sample_points(n_int, n_bc, dim, bounds): # 拉丁超立方采樣內(nèi)部點(diǎn) sampler qmc.LatinHypercube(ddim) sample sampler.random(nn_int) x_int qmc.scale(sample, bounds[0], bounds[1]) # 邊界點(diǎn)在邊界面上均勻采樣 x_bc np.random.uniform(bounds[0], bounds[1], size(n_bc, dim)) # 強(qiáng)制邊界點(diǎn)落在邊界上以二維為例 for i in range(n_bc): if np.random.rand() 0.5: x_bc[i, 0] bounds[0][0] if np.random.rand() 0.5 else bounds[1][0] else: x_bc[i, 1] bounds[0][1] if np.random.rand() 0.5 else bounds[1][1] return x_int, x_bcn_int和n_bc根據(jù)問題規(guī)模調(diào)一維問題 500 和 50 就夠二維建議 2000 和 200三維要 5000 以上。bounds是求解區(qū)域的上下界比如[[0,0],[1,1]]表示二維單位正方形。采樣完最好畫個(gè)散點(diǎn)圖看一眼確認(rèn)沒有點(diǎn)跑到區(qū)域外面去。3.2 訓(xùn)練腳本參數(shù)詳解與第一次運(yùn)行資源里的訓(xùn)練腳本入口是train.py主要參數(shù)包括學(xué)習(xí)率、訓(xùn)練輪數(shù)、網(wǎng)絡(luò)層數(shù)、每層神經(jīng)元數(shù)。學(xué)習(xí)率默認(rèn) 1e-3用 Adam 優(yōu)化器訓(xùn)練 5000 輪后切換到 L-BFGS 做精細(xì)收斂。網(wǎng)絡(luò)結(jié)構(gòu)默認(rèn) 4 層每層 50 個(gè)神經(jīng)元這個(gè)規(guī)模在二維問題上夠用三維問題建議加到 6 層 64 神經(jīng)元。第一次運(yùn)行建議把訓(xùn)練輪數(shù)改成 1000先看損失能不能下降確認(rèn)流程通了再跑完整訓(xùn)練。運(yùn)行命令python train.py --epochs 1000 --lr 0.001 --layers 4 --neurons 50 --n_int 2000 --n_bc 200跑起來之后終端會(huì)打印每 100 輪的損失值正常情況方程殘差損失會(huì)從 1e-1 降到 1e-4 以下邊界損失降到 1e-6 左右。如果損失不降反升先檢查學(xué)習(xí)率是不是太大改成 1e-4 試試如果損失卡在某個(gè)值不動(dòng)可能是采樣點(diǎn)太少或者網(wǎng)絡(luò)容量不夠。我一般會(huì)同時(shí)開 TensorBoard 看損失曲線比盯終端直觀。3.3 結(jié)果驗(yàn)證通量分布對(duì)比與誤差指標(biāo)訓(xùn)練完之后資源里提供了plot.py畫通量分布圖還會(huì)計(jì)算相對(duì) L2 誤差。驗(yàn)證方法有兩種一是和解析解對(duì)比如果有的話二是和傳統(tǒng)數(shù)值方法結(jié)果對(duì)比。中子學(xué)問題很多沒有解析解常見做法是用蒙特卡洛程序跑一個(gè)參考解然后算 PINN 預(yù)測(cè)值和參考值的相對(duì)誤差。資源里默認(rèn)用有限差分結(jié)果做參考誤差控制在 1% 以內(nèi)算合格。下面是一個(gè)計(jì)算相對(duì) L2 誤差的代碼片段def relative_l2_error(phi_pred, phi_ref): # phi_pred: 模型預(yù)測(cè)通量, phi_ref: 參考通量 numerator np.linalg.norm(phi_pred - phi_ref) denominator np.linalg.norm(phi_ref) return numerator / denominatorphi_pred和phi_ref要保證在同一套網(wǎng)格點(diǎn)上取值否則誤差沒意義。如果誤差大于 5%先檢查邊界條件有沒有加對(duì)再看方程殘差的權(quán)重是不是太小。資源里默認(rèn)殘差權(quán)重是 1.0邊界權(quán)重也是 1.0你可以根據(jù)問題調(diào)整比如邊界條件嚴(yán)格時(shí)把邊界權(quán)重調(diào)到 10。4. 避坑指南PINN 訓(xùn)練中子學(xué)問題最容易翻車的五個(gè)地方4.1 損失不收斂方程殘差一直在震蕩現(xiàn)象訓(xùn)練幾百輪后方程殘差損失在 0.1 附近上下跳降不下去。原因?qū)W習(xí)率太大或者采樣點(diǎn)分布不均勻?qū)е绿荻确较蚧靵y。解決把學(xué)習(xí)率降到 1e-4同時(shí)用拉丁超立方重新采樣確保內(nèi)部點(diǎn)覆蓋整個(gè)區(qū)域。如果還不行檢查方程殘差的量綱是不是和邊界條件差太多必要時(shí)做歸一化。4.2 邊界條件失效通量在邊界處明顯偏離現(xiàn)象內(nèi)部區(qū)域擬合得不錯(cuò)但邊界上通量值明顯不對(duì)甚至出現(xiàn)負(fù)值。原因邊界點(diǎn)采樣太少或者邊界損失權(quán)重太低。解決把邊界點(diǎn)數(shù)量增加到內(nèi)部點(diǎn)的 10%邊界損失權(quán)重調(diào)到 5 到 10。另外檢查邊界條件寫的是 Dirichlet 還是 Neumann中子學(xué)里常見的是通量在邊界為零對(duì)應(yīng) Dirichlet 條件直接讓phi_bc趨近于零就行。4.3 二階導(dǎo)數(shù)計(jì)算報(bào)錯(cuò)提示“element 0 of tensors does not require grad”現(xiàn)象運(yùn)行到自動(dòng)微分那一步直接報(bào)錯(cuò)說張量不需要梯度。原因輸入點(diǎn)沒有設(shè)置requires_gradTrue或者在前向傳播之前就把它轉(zhuǎn)成了 numpy 數(shù)組。解決確保x_int是 torch 張量且requires_grad_(True)并且整個(gè)計(jì)算圖里不要出現(xiàn).detach()或.numpy()。這個(gè)坑我踩過好幾次血淚經(jīng)驗(yàn)就是采樣完之后立刻轉(zhuǎn)張量別等到用的時(shí)候再轉(zhuǎn)。4.4 訓(xùn)練速度慢到無法接受GPU 利用率極低現(xiàn)象一個(gè) epoch 要好幾秒GPU 顯存占用很低。原因采樣點(diǎn)太多或者網(wǎng)絡(luò)太大但 batch 設(shè)置不合理導(dǎo)致大量時(shí)間花在數(shù)據(jù)搬運(yùn)上。解決把內(nèi)部點(diǎn)分批送入網(wǎng)絡(luò)每批 500 到 1000 個(gè)點(diǎn)用DataLoader做批處理。另外檢查是不是在 CPU 上跑torch.cuda.is_available()返回 True 的話把模型和輸入都.cuda()。4.5 換一個(gè)工況就要重新訓(xùn)練沒有泛化能力現(xiàn)象在某個(gè)截面參數(shù)下訓(xùn)練好的模型換一組截面參數(shù)后預(yù)測(cè)完全不對(duì)。原因PINN 把截面參數(shù)當(dāng)成了固定值沒有作為輸入。解決把截面參數(shù)也作為網(wǎng)絡(luò)輸入的一部分讓模型學(xué)習(xí)參數(shù)到通量的映射。資源里默認(rèn)只輸入空間坐標(biāo)你可以改成輸入(x, y, z, D, Sigma_a)這樣訓(xùn)練一次就能覆蓋多個(gè)工況。代價(jià)是訓(xùn)練數(shù)據(jù)要包含不同參數(shù)組合采樣點(diǎn)數(shù)量也要增加。5. 進(jìn)階技巧把 PINN 從中子學(xué)擴(kuò)散方程推到輸運(yùn)方程擴(kuò)散方程只是中子輸運(yùn)的近似真正要做高保真計(jì)算還得回到輸運(yùn)方程。輸運(yùn)方程多了角度變量輸入維度變成空間加能量加方向網(wǎng)絡(luò)容量和采樣點(diǎn)數(shù)量都要大幅增加。我一般會(huì)先把角度離散成幾個(gè)方向每個(gè)方向單獨(dú)訓(xùn)練一個(gè)子網(wǎng)絡(luò)最后加權(quán)求和。這樣比直接上高維網(wǎng)絡(luò)穩(wěn)定得多。另一個(gè)技巧是用自適應(yīng)采樣在殘差大的區(qū)域多撒點(diǎn)殘差小的區(qū)域少撒點(diǎn)能顯著提升收斂速度。資源里提供了一個(gè)簡單的自適應(yīng)采樣函數(shù)你可以基于它改def adaptive_sampling(model, x_int, D, Sigma_a, S, n_new100): x_int.requires_grad_(True) phi model(x_int) phi_x torch.autograd.grad(phi, x_int, torch.ones_like(phi), create_graphTrue)[0] phi_xx torch.autograd.grad(phi_x, x_int, torch.ones_like(phi_x), create_graphTrue)[0] residual torch.abs(-D * phi_xx Sigma_a * phi - S) # 選殘差最大的前 n_new 個(gè)點(diǎn)附近加密 _, idx torch.topk(residual.flatten(), n_new) x_new x_int[idx] torch.randn_like(x_int[idx]) * 0.01 return torch.cat([x_int, x_new], dim0)這個(gè)函數(shù)每次調(diào)用會(huì)在殘差大的地方加新點(diǎn)訓(xùn)練過程中每隔 500 輪調(diào)用一次能明顯改善高梯度區(qū)域的擬合精度。驗(yàn)證方法也簡單訓(xùn)練完之后看殘差分布圖如果最大殘差和平均殘差差距在 10 倍以內(nèi)說明采樣比較均勻了。從那以后我每次做 PINN 訓(xùn)練都強(qiáng)制走一遍自適應(yīng)采樣不然總有些局部區(qū)域誤差下不來。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取