CN104317655A - 基于集群式gpu加速的多源全路径蒙特卡罗模拟方法 - Google Patents

Abstract

本发明涉及一种基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，包括：步骤1：Host节点向各Client节点进行计算任务的一次分配；步骤2：基于被分配的计算任务，各Client节点向本节点内的各GPU设备进行计算任务的二次分配；步骤3：基于被分配的计算任务，各GPU设备向其各block进行计算任务的三次分配；步骤4：各Client节点对被分配的相应光子的传输过程进行追踪。本方法通过采用集群式GPU实现了加速的三级并行高速运算，并能够实现对光子在组织中传播时所经各体素的路径信息进行保存，从而获取丰富的光子路径信息。

Description

基于集群式GPU加速的多源全路径蒙特卡罗模拟方法

技术领域

本发明涉及计算机科学，数学仿真和生物医学工程领域，具体涉及一种基于集群式GPU加速的多源全路径蒙特卡罗模拟方法。

背景技术

蒙特卡罗方法的基本特征是对随机性问题进行仿真，用于有效地解决随机性问题，甚至对许多确定性方法所难以解决的随机性问题都能够进行较为方便的解决。在生物组织中的光子传输问题的研究中，蒙特卡罗方法也显示出其相对于其它方法(如扩散近似)的优势：严密且灵活，不受组织的光学特性的限制，易于使用，模拟精度高，因此，其被该领域研究者奉为金标准。

自Wilson和Adam于1983年首次将蒙特卡罗方法引入到激光与生物组织之间相互作用的研究领域以来，越来越多的研究者投入到该研究中。绝大多数生物组织具有复杂的解剖结构，属于三维非均匀介质。通过蒙特卡罗方法模拟光在生物组织中传输特性，是为了重建组织体某一个横截面甚至整个三维空间的光学参数图像而服务的，以便于对生物组织光学检测成像系统进行设计和优化，进而实现对疾病的诊断、监测及光学临床治疗。为了实现对组织体三维空间的光学参数图像的重建，必须获取光子在该组织体中传播时，所经各体素的路径信息。

到目前为止，被报道的光子传输三维蒙特卡罗代码有三份，一份是Boas小组于2002年研发的，名为“tMCimg”,可应用于基于MRI和CT获取的数字影像数据搭建的体素介质模型。一份是Margallo-Baldas等人于2007年研发的“trMC”程序，该程序可应用于基于microCT获取的三角面片介质模型。另一份是骆清铭小组于2010年研发的MCVM程序，该程序可同样应用于MRI和CT获取的数字影像数据搭建的体素介质模型。

但是，以上的三份代码均不能获取光子在组织体中传播时所经各体素的路径信息。并且，随着模拟光子数目的不断增加，模拟源数目的不断增加，以及组织越来越复杂的结构，导致计算机运算时间也在成倍地增加。

有鉴于此，本发明提出了一种优化的程序结构：基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，该方法通过采用集群式GPU实现对加速的三级并行高速运算，并能够实现对光子在组织中传播时所经各体素的路径信息进行保存。

发明内容

本发明提供的一种基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，通过采用集群式GPU实现加速的三级并行高速运算。

本发明提供的一种基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，包括：

步骤1：基于MPI消息传递通信协议，Host节点向各Client节点进行计算任务的一次分配，并完成相应的数据交互；

步骤2：基于被分配的计算任务，各Client节点向本节点内的各GPU设备进行计算任务的二次分配，并完成相应的数据交互；

步骤3：基于被分配的计算任务，各GPU设备向其各block进行计算任务的三次分配，并完成相应的数据交互；

步骤4：各Client节点对被分配的相应光子的传输过程进行追踪。

所述步骤1为：

步骤1.1：各Client节点获取本节点的GPU设备的数量信息；

步骤1.2：各Client节点向其他节点广播本节点GPU设备的数量信息；

步骤1.3：Host节点接收各Client节点的广播信息，并获得各节点GPU设备的数量信息；

步骤1.4：Host节点根据各Client节点的GPU设备的数量以及与光源的光子数量相关的计算任务量，向各Client节点分配平均到每个GPU的计算任务。

所述步骤2为：

步骤2.1：各Client节点CPU的主进程根据本节点GPU设备的数目开辟相应数目的子进程，每个子进程对应一个本节点的GPU设备；

步骤2.2：各Client节点CPU的主进程将计算任务分配给子线程进行执行。

在所述步骤2.2中，各Client节点的子进程在本节点多核CPU的不同的核上并行运行。

在所述步骤2.2中，各Client节点的子进程在本节点单核CPU的单核上串行运行。

所述步骤3为：

步骤3.1：各Client节点确定本节点各GPU设备SM上的可用资源；

步骤3.2：各Client节点根据本节点GPU设备中各block使用资源，确定active block数量和active warp数量；

步骤3.3：各Client节点根据本节点GPU设备中各block使用资源，确定各block的尺寸和维度；

步骤3.4：按照计算任务的规模，各Client节点确定本节点各GPU设备grid中的block数量；

步骤3.5：各Client节点中的各GPU设备将计算任务分配给本设备的每个block。

所述步骤4为：

步骤4.1：启动cuda，各Client节点CPU对内存空间、显存空间进行分配，并将要计算的数据从内存复制到显存上；

步骤4.2：将入射光源表征为设定数目光子的集合，将入射光源位置和入射光方向赋给每个光子作为其初始位置和方向；

步骤4.3：追踪每个光子的传输过程；

步骤4.4：CPU分配内存，用于存放GPU输出数据，并将计算后显存上的数据复制到内存上，并将本节点的运行信息返回给host节点。

所述步骤4.3追踪每个光子的传输过程，包括以下情况：

当光子移动至生物组织边界时，发生反射或折射；

当光子折射到生物组织外时，停止追踪；

当光子因散射而改变移动方向时，继续追踪光子下一步的移动直到光子能量衰减至阙值以下。

所述步骤4.3可以与步骤4.1中的CPU将要计算的数据从内存复制到显存上的操作同时进行。

光子的移动步长由随机数确定。

在步骤4.3中，采用俄罗斯轮回盘算法对光子进行追踪：如果光子在轮回盘中复活，则继续追踪光子下一步的移动；如果光子在轮回盘中死亡，则停止对光子的追踪。

对光吸收量数据的访问采用原子操作，当多个线程同时访问全局存储器的同一位置时，保证每个线程能够实现对共享可写数据的互斥操作：在一个操作完成之前，其它任何线程都无法访问此地址，从而保证记录的光吸收量数据的绝对正确性。

在显存不足的情况下，Client节点将部分数据从显存复制到内存，在GPU并行执行程序的同时，CPU串行执行程序。

所述方法还包括：

步骤5：光子追踪结束后，Host节点输出光吸收矩阵和所有逸出光子信息及光子传播所经各体素的路径信息。

所述光子传播所经各体素的路径信息，由追踪每个光子传输过程中记录的每个光子所经体素的索引值及经过该体素光子所走的路径长度得到。

所述方法还包括：

步骤6：对步骤5的结果进行统计，转换为其它所需获取的传输特性，包括：将光吸收量矩阵中各元素除以其对应体素的吸收系数，得到通量分布；将光吸收量矩阵所有元素的值相加，然后除以光子数目得到光吸收量；将所有在光源同侧逸出的光子当前能量相加，然后除以光子数目得到总透射量。

综上所述，本发明提出了一种基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，该方法通过采用集群式GPU实现了加速的三级并行高速运算，并能够实现对光子在组织中传播时所经各体素的路径信息进行保存，从而获取丰富的光子路径信息，为生物组织光学检测成像系统设计的优化提供了方便，并为光治疗方法提供了精确的指导信息。

附图说明

参照如下附图将更加易于理解本发明：

图1为本发明提供的一种基于集群式GPU加速的多源全路径蒙特卡罗模拟方法的整体系统框架图。

图2为本发明提供的一种基于集群式GPU加速的多源全路径蒙特卡罗模拟方法的基本流程图。

具体实施方式

首先，将目标生物组织的空间结构描述为一个三维数字矩阵，即组织模型。矩阵中的元素对应目标生物组织的体素，每个体素的数值为标识组织类型的数字。体素越小，描述组织结构的组织模型越逼近于真实生物组织结构，模拟得到的光子传输特性精度就越高。

然后，设定一个与组织模型大小相同的三维空矩阵，其中每个元素用来记录其对应生物组织中的一个体素对光的吸收量。并设定目标生物组织的各类组织特性参数，包括吸收系数，散射系数，折射系数和各向异性因子。这些参数值的精度也会影响模拟精度，可查阅相应文献来获取这些参数。值得注意是，这些参数的选定要求与目标光源波长一致。

获得上述编程模型后即可采用本发明方法进行计算。

如图1所示，本方法的整体系统框架为Host-client结构，包括用以分配和控制计算任务的Host端和至少两个用以接收计算任务进行运算的Client节点，每个Client节点可将运算结果上传给所述的Host端，相互间能进行通信。该模拟方法实现了三层次并行，分别是集群分布式并行，多核并行和GPU内多线程并行。

如图2所示，本发明的基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，包括：

步骤1：基于MPI消息传递通信协议，Host节点向各Client节点进行计算任务的一次分配，并完成相应的数据交互；

可选地，所述步骤1可以包括：

步骤1.1：各Client节点获取本节点的GPU设备的数量信息；

步骤1.2：各Client节点向其他节点广播本节点GPU设备的数量信息；

步骤1.3：Host节点接收各Client节点的广播信息，并获得各节点GPU设备的数量信息；

步骤1.4：Host节点根据各Client节点的GPU设备的数量以及与光源的光子数量相关的计算任务量，向各Client节点分配平均到每个GPU的计算任务。

步骤2：基于被分配的计算任务，各Client节点向本节点内的各GPU设备进行计算任务的二次分配，并完成相应的数据交互；

可选地，所述步骤2为：

步骤2.1：各Client节点CPU的主进程根据本节点GPU设备的数目开辟相应数目的子进程，每个子进程对应一个本节点的GPU设备；

步骤2.2：各Client节点CPU的主进程将计算任务分配给子线程进行执行。可选地，各Client节点的CPU可以是多核CPU或者单核CPU，其中，在多核CPU情况下，各Client节点的子进程在本节点多核CPU的不同的核上并行运行；在单核CPU情况下，各Client节点的子进程在本节点单核CPU的单核上串行运行。

步骤3：基于被分配的计算任务，各GPU设备向其各block进行计算任务的三次分配，并完成相应的数据交互；

可选地，所述步骤3为：

步骤3.1：各Client节点确定本节点各GPU设备SM上的可用资源；

步骤3.2：各Client节点根据本节点GPU设备中各block使用资源，确定active block数量和active warp数量；

步骤3.3：各Client节点根据本节点GPU设备中各block使用资源，确定各block的尺寸和维度；

步骤3.4：按照计算任务的规模，各Client节点确定本节点各GPU设备grid中的block数量；

步骤3.5：各Client节点中的各GPU设备将计算任务分配给本设备的每个block。

步骤4：各Client节点对被分配的相应光子的传输过程进行追踪。

可选地，所述步骤4为：

步骤4.1：启动cuda，各Client节点CPU对内存空间、显存空间进行分配：将x,y,z方向体素的索引值使用纹理存储器存储，先对其进行纹理参考声明，再对其进行三维纹理绑定；将组织模型在x,y,z方向的尺寸以及及组织的光学特性参数用常量存储器存储；将光子传播所经各体素的路径信息用页锁定内存存储；其他变量则用全局存储器存储。然后将要计算的数据从内存复制到显存上；

步骤4.2：将入射光源表征为设定数目光子的集合，将入射光源位置和入射光方向赋给每个光子作为其初始位置和方向；

步骤4.3：追踪每个光子的传输过程：当光子移动至生物组织边界时，发生反射或折射；当光子折射到生物组织外时，停止追踪；当光子因散射而改变移动方向时，继续追踪光子下一步的移动直到光子能量衰减至阙值以下。其中，光子的移动步长可以由随机数确定。可选地，在本步骤中，可以将其中耗费时间长的大规模数据并行、高计算密度的步骤安排在GPU上并行执行，其余步骤安排在CPU上串行执行。

可选地，在步骤4.3中可以采用俄罗斯轮回盘算法对光子进行追踪：如果光子在轮回盘中复活，则继续追踪光子下一步的移动；如果光子在轮回盘中死亡，则停止对光子的追踪。

步骤4.4：CPU分配内存，用于存放GPU输出数据，并将计算后显存上的数据复制到内存上，并将本节点的运行信息返回给host节点。

可选地，所述步骤4.3可以与步骤4.1中的CPU将要计算的数据从内存复制到显存上的操作同时进行。

可选地，对光吸收量数据的访问可以采用原子操作，当多个线程同时访问全局存储器的同一位置时，保证每个线程能够实现对共享可写数据的互斥操作：在一个操作完成之前，其它任何线程都无法访问此地址，从而保证记录的光吸收量数据的绝对正确性。

可选地，在显存不足的情况下，Client节点将部分数据从显存复制到内存，在GPU并行执行程序的同时，CPU串行执行程序。

可选地，所述方法还可以包括：

步骤5：光子追踪结束后，Host节点输出光吸收矩阵和所有逸出光子信息及光子传播所经各体素的路径信息。其中，所述光子传播所经各体素的路径信息，由追踪每个光子传输过程中记录的每个光子所经体素的索引值及经过该体素光子所走的路径长度得到。

可选地，所述方法还可以包括：

步骤6：对步骤5的结果进行统计，转换为其它所需获取的传输特性，包括：将光吸收量矩阵中各元素除以其对应体素的吸收系数，得到光通量分布；将光吸收量矩阵所有元素的值相加，然后除以光子数目得到光吸收量；将所有在光源同侧逸出的光子当前能量相加，然后除以光子数目得到总透射量。

Claims (10)

1.一种基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，其特征在于，所述方法包括：

步骤1：基于MPI消息传递通信协议，Host节点向各Client节点进行计算任务的一次分配，并完成相应的数据交互；

步骤2：基于被分配的计算任务，各Client节点向本节点内的各GPU设备进行计算任务的二次分配，并完成相应的数据交互；

步骤3：基于被分配的计算任务，各GPU设备向其各block进行计算任务的三次分配，并完成相应的数据交互；

步骤4：各Client节点对被分配的相应光子的传输过程进行追踪。

2.根据权利要求1所述的基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，其特征在于，所述步骤1为：

步骤1.1：各Client节点获取本节点的GPU设备的数量信息；

步骤1.2：各Client节点向其他节点广播本节点GPU设备的数量信息；

步骤1.3：Host节点接收各Client节点的广播信息，并获得各节点GPU设备的数量信息；

步骤1.4：Host节点根据各Client节点的GPU设备的数量以及与光源的光子数量相关的计算任务量，向各Client节点分配平均到每个GPU的计算任务。

3.根据权利要求2所述的基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，其特征在于，所述步骤2为：

步骤2.1：各Client节点CPU的主进程根据本节点GPU设备的数目开辟相应数目的子进程，每个子进程对应一个本节点的GPU设备；

步骤2.2：各Client节点CPU的主进程将计算任务分配给子线程进行执行。

4.根据权利要求3所述的基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，其特征在于，在所述步骤2.2中，各Client节点的子进程在本节点多核CPU的不同的核上并行运行。

5.根据权利要求3所述的基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，其特征在于，在所述步骤2.2中，各Client节点的子进程在本节点单核CPU的单核上串行运行。

6.根据权利要求3所述的基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，其特征在于，所述步骤3为：

步骤3.1：各Client节点确定本节点各GPU设备SM上的可用资源；

步骤3.2：各Client节点根据本节点GPU设备中各block使用资源，确定active block数量和active warp数量；

步骤3.3：各Client节点根据本节点GPU设备中各block使用资源，确定各block的尺寸和维度；

步骤3.4：按照计算任务的规模，各Client节点确定本节点各GPU设备grid中的block数量；

步骤3.5：各Client节点中的各GPU设备将计算任务分配给本设备的每个block。

7.根据权利要求6所述的基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，其特征在于，所述步骤4为：

步骤4.1：启动cuda，各Client节点CPU对内存空间、显存空间进行分配，并将要计算的数据从内存复制到显存上；

步骤4.2：将入射光源表征为设定数目光子的集合，将入射光源位置和入射光方向赋给每个光子作为其初始位置和方向；

步骤4.3：追踪每个光子的传输过程；

步骤4.4：CPU分配内存，用于存放GPU输出数据，并将计算后显存上的数据复制到内存上，并将本节点的运行信息返回给host节点。

8.根据权利要求7所述的基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，其特征在于，所述方法还包括：

步骤5：光子追踪结束后，Host节点输出光吸收矩阵和所有逸出光子信息及光子传播所经各体素的路径信息。

9.根据权利要求8所述的基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，其特征在于，所述光子传播所经各体素的路径信息，由追踪每个光子传输过程中记录的每个光子所经体素的索引值及经过该体素光子所走的路径长度得到。

10.根据权利要求8所述的基于集群式GPU加速的多源全路径蒙特卡罗模拟方法，其特征在于，所述方法还包括：

步骤6：对步骤5的结果进行统计，转换为其它所需获取的传输特性，包括：将光吸收量矩阵中各元素除以其对应体素的吸收系数，得到通量分布；将光吸收量矩阵所有元素的值相加，然后除以光子数目得到光吸收量；将所有在光源同侧逸出的光子当前能量相加，然后除以光子数目得到总透射量。