这项工作的难点在于如何将蒙特卡罗粒子输运（下简记MC输运）模拟计算可扩展点分离出来，在可扩展点设计接口，支持用户扩展个性功能。因此，我们从MC输运求解的积分形式粒子输运方程入手，分析和分离MC输运中的可扩展点。

稳定的粒子输运方程的积分形式如式（1）所示，其中ϕ是未知数，表示相空间Π对应的通量，或者说单位时间内通过的粒子数。Π是相空间的简写，在这里，它表示粒子的位置ρ、能量Ε和飞行方向 Ω，是一个六维相空间。

式（1）有Neumann级数解

其中

m表示从粒子源出发经过μ次碰撞后到达Π处的通量（粒子数）。

在Neumann级数解收敛的情况下，即存在一个Μ值，使得m大于Μ时，都是零，就获得了$ \phi $的求解，这就是可以使用MC方法通过模拟粒子行为而获得了Neumann级数解的理论基础^[18-19]。

实际上，核物理领域关心的物理量不仅是通量，更多的是各种响应量（g为已知分布），即

ΜΧ方法是可以同时计算多个响应量的，所以式（4）数量不定。将式（3）代入式（4）的任意一个，都可以有

对于式（5）， ΜΧ方法抽取大量粒子样本，通过统计方式求平均值来代替积分，即

式中：$ I_{m}^{(j)} $表示第ϕ个粒子，从粒子源中抽样确定初始状态$ {P}_{0} $，根据分布K，抽样确定$ {P}_{1} $，直到第μ次确定的状态正好是确定P（这就是该粒子对$ {\phi }_{m} $的小贡献），再乘$ g(P) $。

由于粒子事件的独立性，在获得$ I_{m}^{(j)} $的同时，也可以获得$ I_{1}^{(j)},I_{2}^{(j)},\cdots I_{m-1}^{(j)} $，就是说粒子发生了μ次碰撞，每次都对不同的$ {\phi }_{m} $有贡献。

如此，我们就获得了MC方法求Neumann级数解的一个简单流程：

（1）从粒子源中抽样获得一个粒子的$ {P}_{0} $；

（2）根据粒子状态$ {P}_{i},i=0,1,2,\cdots $，从分布Κ抽样确定$ {P}_{i}{}_{+1} $，乘上$ {g}_{{\mathrm{A}}},{g}_{{\mathrm{B}}},\cdots $，计数；

（3）重复第2步，直到粒子死亡或离开系统；计数。回到第1步；

（4）所有粒子都模拟完，进行统计处理和输出结果。

一般中子光子输运的MC程序都是按照上述流程，但这只是一个非常简化的流程，真实情况要远比这复杂和多变。

说明一点，算符Κ展开为

式（7）在MC算法代表了从一个相空间状态抽取另一个相空间的分布，它可以分成两个分布：在保持位置不变的情况下，能量和方向发生变化的分布，物理上就是粒子发生了核反应碰撞；在保持能量和方向不变的情况下，位置发生变化的分布，物理上就是粒子沿飞行方向飞行了一段距离，改变了位置。

从式（4）可知，MC程序的响应量计算是可扩展点，对于相同的输运模拟，不同的物理量，响应函数$ g(P) $显然是不同的，这来自用户需求，例如临界计算需要进行k_eff本征值计数，而固定源计算则不需要。

另一方面，MC输运一般伴随着降方差技巧的使用，这些技巧会改变粒子的权重，在响应物理量的最终计算上有所体现。降方差技巧多种多样，也是MC输运方法的一类研究内容，属于MC输运程序的可扩展点。

剔除统计计数和降方差技巧后，一般的MC粒子输运过程都可以整理成如图1所示的计算流程，包含了从外到内的6层的循环，前三个分别是：（1）时间步循环，涉及到多个时间步的输运计算时适用，输运-燃耗-热工的耦合在此处实现；（2）粒子代间循环，例如临界计算问题多代粒子，或固定源问题多个批次粒子；（3）粒子数循环，对一代/批次的粒子依次进行循环，其中每个粒子都要从粒子库中获得，若是初始源粒子则从粒子源中抽取。

在前三个循环中有两个可扩展点：第一个是输运问题，或者说粒子类型。在一定程度上，输运问题是有限的，只有临界和固定源两种，其中临界问题只有中子一种粒子类型，分成k_eff本征值临界和α本征值临界两种情况，固定源问题则有多种粒子类型，如中子、光子、电子、质子等；第二个可扩展点是粒子源项，对应公式（5）中的$ S({P}_{0}) $，这是与用户或工程应用相关性很大的一部分，现有流行的大部分MC软件都支持用户自定义源项，如MCNP、JMCT等。

图1显示的流程中后三层的循环：（4）历史事件循环、（5）碰撞核素循环和（6）核素反应道循环，他们共同构成了一个粒子的完整模拟过程，也称为完整历史，对应着式（5）中的算符K。每个粒子开始模拟后，第4层历史事件循环负责处理粒子所发生的历史事件处理，每个循环步都会移动粒子的空间位置，对应着算符K展开式即式（7）中的T；当粒子发生碰撞事件时，则进入第5层的碰撞核素循环，进行碰撞核素的抽样，这是粒子碰撞处理的数值部分，可以认为对应着算符K展开式即式（7）中的C的数值部分；当确定碰撞核素后，则就进入第6层核素反应道循环，在核素内部的具体反应道中进行抽样，这部分与核素物理性质及物理模型密切相关，对应着算符K展开式即式（7）中的$ \mathit{\Sigma } $和C的物理部分。

实质上，这三层循环对应着就是三个可扩展点，例如核素反应道循环，可以采用连续能量模式也可以采用多群模式；碰撞核素循环可以采用一些数值优化方法来提交效率；历史事件也可采用单碰撞方法也可采用浓缩历史方法等，这都是可变的。

而且，这三者合到一起成为一个大的可扩展点，就是粒子模拟历史，例如可以实现不发生碰撞的权窗参数计算历史或者某种特殊需求的模拟历史等。

把上述三层变化分析出来的可扩展点合到一起，就可以看出式（5）中的被积分项都是可变的，剩下的就只有积分算符了，这对应的软件代码就是软件架构，六层循环架构以及不同类之间的组织关系、各种容器管理器、数据结构、接口及I/O。在JMCT软件中，这些软件架构代码占据10%左右。

为了最大限度的实现可扩展性，基于上述分析，我们设计出新的MC输运软件架构，主要包括微内核、信息池、特征组件接口和特征组件，以下逐一介绍：

（1）微内核：MC输运中稳定的、不可扩展的流程架构组成，主要对应图2所示的多层循环架构和粒子管理器。

（2）信息池：用于存储输运中特征组件所需要的各种信息，分成两类：第一类是与当前MC粒子状态极为相关的，直接影响输运模拟流程的信息，包括能量、位置、飞行方向、权重、类型等；第二类是为了满足各种特征组件的功能需求而存储的信息，包括粒子发生过的事件信息、过去的粒子状态、粒子来源、产地、父亲或用户标记等。其中第二类具有可扩展性，允许用户自定义信息存储其中。

（3）特征组件接口：用于支持用户研发特征组件的接口类，规范了特征组件的成员函数的输入输出参数和函数调用位置。目前具有源粒子抽样、粒子历史、输运、碰撞器、物理模型、统计计数和降方差技巧共计7个接口。这7个接口在微内核输运方程上都有对应的注册函数，可以将基于该接口开发的特征组件安装到微内核上。

（4）特征组件：基于特征组件接口开发的具有某种单一功能的函数类，该功能体现在修改或读取信息池的某个或某些数据，例如输运特征组件会修改当前粒子的能量、位置、方向等，物理模型特征组件会修改粒子的能量、统计计数特征组件会读取信息池的粒子过去和现在状态来计算统计量并计数等；对应特征组件接口，特征组件也分为7类。

四者关系如图所示，整个软件架构只有一个微内核，但可以有多个特征组件，特征组件通过继承特征组件接口而研发而来，组装到微内核上形成完整的输运方程。微内核按照其内部流程架构，通过回调特征组件接口的函数，实现特征组件的功能调用，完成粒子模拟。模拟其间，微内核与各个特征组件之间通过信息池来互通信息。

从JMCT软件出发进行架构重构，我们实现了上述MC输运软件架构设计的落地，研制出MCPT输运平台，如图3所示。

MCPT主要分成微内核和特征组件两层。MCPT的微内核包含了输运信息池，就是上节中提到的第一类信息池，由MCPT负责维护。微内核对用户提供了MC输运方程组装接口，它含有多个注册函数，用来组装各种特征组件；初始化函数，用来读取输入文件对方程进行初始化；求解函数，对MC输运方程进行求解；将计算结果进行规约输出的后处理函数。

MCPT微内核提供了历史类等7类特征组件接口，对应扩展出的7类特征组件形成了MCPT的特征组件层。这些特征组件来自JMCT软件对应的功能模块，经过了各种算例的反复验证。

（1）物理模型接口：对应粒子的核反应物理模型，相当于输运方程的碰撞核（式（7）中的C）

（2）输运算法接口：对应粒子的空间输运算法，相当于输运方程的碰撞核（式（7）中的T）

（3）源抽样接口：对应模型粒子源，相当于输运方程的源项（式（5）中的S）

（4）碰撞器接口：对应粒子碰撞反应的第一层处理，包括输运方程的总截面/分截面计算和抽样（式（7）中的$ \mathit{\Sigma } $）

（5）历史类接口：对应一个粒子出生后的完整模拟，相当于式（5）中的算符K

（6）统计计数接口：对应着各种物理响应量计数，对应着输运方程的响应函数（式（5）中的G）

（7）降方差技巧接口：对应着输运过程的降方差技巧，输运方程上没有对应项

限于篇幅，此处仅以统计计数接口及特征组件来展开介绍。

统计计数接口(Effector)是MCPT的一个典型设计，它可以说粒子输运MC模拟中除了不可或缺的粒子输运流程外，其他所有可附加的物理效应的统一接口，不仅支持开发各种统计计数，也可以开发降方差技巧，甚至个性物理模型和其他特殊用户需求。

如图4所示，Effector接口提供了很多成员函数，这些成员函数在粒子输运模拟的流程中对应的位置被调用，例如粒子碰撞前后的操作、或者粒子跨越几何界面前后的操作等，都是粒子属性状态发生改变，或者可能触发统计计数的位置。成员函数本身具有默认的空实现，用户在实现派生子类时，只需重载自己需要的成员函数即可。

Effector接口的设计，实现了粒子输运过程本身的代码和统计计数的解耦分离，可以由各自的相关人员开发；将一个物理效应的碎片化代码归并在一个类文件中管理，对开发人员而言不仅利于代码维护，更是一种权益保护；只需在main函数中一句注册语句就可以将特征组件组装在MCPT上，简单易用；用户基于MCPT开发定制软件，只需编写个性的Effector特征组件，不用编写粒子输运模拟的代码，减少了编码工作量，可节省时间。Effector接口设计完整地诠释了MCPT的目标和意义。

多个统计计数特征组件可以一起组装到MC输运方程上，其相同的成员函数执行顺序依据注册顺序。

7类特征组件之间需要配合使用，彼此之间是有适配性的，图5给出了中子固定源、中子临界本征值、中子时间本征值、光子固定源等四类输运问题的特征组件配置图。

2020年左右，MCPT推出1.0版，提供了多个历史类、碰撞器类和物理模型类的特征组件，特别是统计计数和降方差技巧特征组件近20多个。这些特征组件都来自JMCT的代码重构，复现了JMCT大约70%的功能，通过了近百个JMCT已经成功验证应用过的算例^[1]的正确性测试，计算结果与JMCT完全一致。

基于MCPT已有的特征组件，可直接组装出中子、光子类型的临界k_eff，临界α及固定源问题的专用功能并行输运软件（图5）。针对某项目辐射效应模拟需求，我们使用MCPT开发了专用的中子辐射效应、高能伽玛辐射效应、可见光辐射效应、辐射电磁耦合等4个轻量级程序。新开发代码（main函数）仅有负责组装特征组件和主流程控制部分，实现行数都在200行以内，很大程度上减少了代码开发量。专用功能软件用途单一，结构简单，操作简单，改动频率低，易于维护，更适合成熟的工程应用。

2020年，MCPT在中国工程物理研究院内各所开始应用。

中物院高性能数值模拟软件中心基于MCPT研制了光子在大气的输运软件。该软件需要模拟高能光子与大气分子的相互作用，而MCPT本身提供的光子物理模型特征组件仅有光子与原子的相互作用模型，不满足需求，因此需要开发专用特征组件。研发人员采用Effector类研制了光子分子相互作用物理模型，将其组装到MCPT上，实现了可以同时模拟光子在大气中输运时能同时与原子和大气分子相互作用的输运过程。

简单介绍如下：

（1）从Effector接口上派生一个组件类UserEffector；

（2）在UserEffector内实现actionAfterTotalCX成员函数，在MCPT计算的总截面上增加粒子与大气分子相互作用的截面，这样总截面就同时考虑光子与原子和大气分子碰撞的两种情况。该函数调用位置是在MCPT输运模拟中计算材料的总截面之后；

（3）在UserEffector内实现actionBeforeCollisionToMaterial成员函数，抽样确定光子是否与大气分子发生碰撞，若是则根据相应的物理模型进行处理，获得出射粒子的能量和方向，确定反应类型，返回。调用位置是对整个物质材料进行碰撞处理前；

（4）编写main函数，组装光子固定源历史、光子碰撞器、连续光子物理、标准源抽样等必要的特征组件和UserEffector组件，如下：

编译链接生成可执行程序。

该程序在模拟过程中，计算总截面后就会调用UserEffector的actionAfterTotalCX函数，总截面被修改；根据总截面抽样判断是否粒子发生碰撞；确定发生碰撞后，进行碰撞处理前，会调用UserEffector的actionAfterTotalCX函数。根据反应类型来判断光子是否已经完成碰撞，若没有（没有发生光子与大气分子的碰撞）则继续处理光子与原子的碰撞抽样（连续光子物理特征组件的函数），若有则不进行光子与原子的碰撞处理。如此实现了在MCPT提供的共性的核反应物理模型基础上，增加了个性个核反应物理模型。

战略研究中心基于MCPT研制了专用辐射探测模拟程序^[6]。研发人员取用了MCPT提供的中子固定源历史类、中子固定源历史、中子碰撞器、连续中子物理、光子固定源历史、光子碰撞器、连续光子物理、标准源抽样等特征组件，自行实现了基于Effector接口开发的NaI探测器探测组件，在main函数中，通过注册函数将上述特征组件进行注册后，编译即专用辐射探测模拟程序。为了实现 NaI 探测器探测效率模拟，研发人员在NaI 探测器探测组件中，通过统计进入 NaI 探测器有效探测区域内的光子触发探测信号的数量与总光子数目之比确定探测器探测效率。由于光子在探测区域内多次发生能量沉积经过的时间仅有几ns，而探测器内光子沉积能量，产生荧光，触发电子学信号的时间尺度达到 ms 量级，光子在探测区域内发生多次能量沉积时一般只能触发一个电子学信号，因此可以认为光子在探测区域内次数与测量结果无关。因此在模拟过程中，可以通过统计在有效探测区域内碰撞次数大于1的光子数目，给出探测器信号数。探测信号统计流程框图见文献[20]。

此外，还有多个专用程序已进入实际应用中：

特殊核材料建模与模拟程序：用户基于Effector接口，开发了针对特殊材料的详细中子裂变模拟的特征组件，实现计算不同能量中子引发的裂变产生的次级粒子个数分布，缓发中子的出射时间等。还开发了匹配的中子计数统计类，完成了特殊核材料的探测模拟。

质子输运程序：用户基于History接口，实现了质子历史的完整模拟，包括质子碰撞模块及质子反应截面模块。如果质子反应产生中子、光子等其他粒子，可以组装MCPT提供的粒子历史类来实现耦合输运。目前，该程序已支撑院内某所用户进行质子输运模拟。

ICF项目实验诊断仿真程序：该程序需要粒子源按照ICF辐射参数进行发射，输运中粒子同样按照ICF辐射参数发生吸收反应。用户分别基于Sampler和Effector接口，实现了自定义的粒子源模块和吸收反应模块，编码约1300行，组合到MCPT上，获得了中子、X光对ICF内爆靶进行成像的诊断后处理能力。

某闪光照相专用程序：该程序需要单独记录初级光子的直穿项贡献，为此用户基于Effector接口和粒子状态类ParticleState的粒子父母属性和碰撞次数信息，开发了能区分初级光子和次级光子的特征组件，实现了初级光子的直穿贡献计数。

针对JMCT的接口开发和可维护性问题，本文从粒子输运方程入手，分析了蒙特卡罗粒子输运模拟过程，从理论上找出了MC粒子输运模拟软件存在的可扩展点，设计了MC粒子输运软件的高可扩展的软件架构。该架构包括稳定的微内核、信息池、特征组件接口和特征组件四部分，实现了特征分离的解耦研发，且在开发过程中无需对方源码。编译后的二进制组件制品，可以通过规范的特征组件接口，组装到二进制微内核制品上，通过二次编译，产生一个并行蒙特卡罗粒子输运程序。这种研发方式不仅实现了不同功能模块的研发解耦，而且能够对微内核和特征组件两方面的源码进行保护。

我们在蒙特卡罗粒子输运解算器库MCPT上实现了上述设计。MCPT提供了一个凝练了蒙卡输运模拟中不变项的微内核，同时开放了多种物理特征组件的编程接口，支持用户柔性扩展个性的物理特征组件。同时也提供了多个通用的数值特征组件和共性的物理特征组件供用户选用；

MCPT已经通过实际应用的验证，成功支持用户开发了多款专用软件。基于MCPT，JMCT软件重构成平台型软件，可扩展性和可维护性都有所提高。经初步评估，通过MCPT和用户自编代码的联合编译而产生用户专用定制的蒙特卡罗粒子输运模拟软件，可节省用户30%～95%的编码工作量。

这种从数值方程分析入手进行软件架构设计模式，对于软件研发，尤其是多物理耦合软件、多人研发或维护的复杂架构软件，具有参考价值。