CN103324780B - 颗粒流动仿真系统和方法 - Google Patents

颗粒流动仿真系统和方法 Download PDF

Info

Publication number
CN103324780B
CN103324780B CN201310195595.6A CN201310195595A CN103324780B CN 103324780 B CN103324780 B CN 103324780B CN 201310195595 A CN201310195595 A CN 201310195595A CN 103324780 B CN103324780 B CN 103324780B
Authority
CN
China
Prior art keywords
particle
gpu
computing node
node
grid
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201310195595.6A
Other languages
English (en)
Other versions
CN103324780A (zh
Inventor
杨磊
齐记
田园
高笑菲
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Modern Physics of CAS
Original Assignee
Institute of Modern Physics of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Modern Physics of CAS filed Critical Institute of Modern Physics of CAS
Priority to CN201310195595.6A priority Critical patent/CN103324780B/zh
Publication of CN103324780A publication Critical patent/CN103324780A/zh
Application granted granted Critical
Publication of CN103324780B publication Critical patent/CN103324780B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2210/00—Indexing scheme for image generation or computer graphics
    • G06T2210/52—Parallel processing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2210/00—Indexing scheme for image generation or computer graphics
    • G06T2210/56—Particle system, point based geometry or rendering

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了基于GPU的颗粒流动仿真系统及方法，包括：根据客户端输入的颗粒建模信息生成颗粒信息，并且生成几何体信息；接收颗粒信息和几何体信息，根据颗粒的数目及各个计算节点中空闲的GPU数目，确定使用哪些计算节点中的哪些GPU，然后根据确定的GPU的数目及颗粒在空间中的分布情况确定哪些颗粒由哪个计算节点的哪个GPU进行计算，并根据确定结果进行分配；在多个GPU中并行计算颗粒碰撞导致的每个颗粒的受力，进而计算出加速度，以仿真颗粒流动；展示仿真结果。根据本发明实施例，能够实现高密度颗粒虚拟实验仿真，并在降低能耗的同时提高运算效率。

Description

颗粒流动仿真系统和方法
技术领域
本发明涉及颗粒流动仿真技术领域，具体地，涉及一种可应用于颗粒物质或固体结构研究的基于GPU的颗粒仿真系统和方法。
背景技术
颗粒系统一直是人们关注的研究内容。在工业领域如：食品控制、化学、土木工程、油气、采矿、制药、粉末冶金、能源等领域都有大量应用；在理论研究方面，如何堆积能达到最密集堆积、沙堆在什么情况下坍塌以研究雪崩等问题。为研究相关颗粒系统，人们需要搭建大型实验用颗粒系统，费时费力。并且某些颗粒系统由于成本高，需要在极端的条件下运行，不可能通过实验搭建完成，而基于虚拟实验的仿真系统则不存在类似的问题。
目前颗粒系统模拟的计算方法主要以DEM(离散单元算法)方法为主。DEM方法是继有限元法、计算流体力学(CFD)之后，用于分析物质系统问题的又一种数值计算方法。DEM方法通过建立微元体系的参数化模型，进行颗粒行为模拟和分析，为解决众多涉及颗粒、结构、流体与电磁及其耦合等综合问题提供了一个平台，已成为科学过程分析、产品设计优化和研发的一种强有力的工具。目前DEM方法除了在科学研究中运用，也在科技应用和工业领域逐渐成熟，并已从颗粒物质研究、岩土工程和地质工程等科学及应用拓展至工业过程与工业产品的设计、研发领域，在诸多工业领域取得了重要成果。
DEM方法的特点是仿真精度高，但计算量大。当前，DEM方法的实现以CPU为主，这些方法由于CPU计算能力不足导致计算规模不足，在可以接受的机时内只能计算很小的空间尺寸和时间尺寸；或者因建设成本高，需要建设大规模甚至超大规模的CPU计算机集群，而且耗电量过大，使用和维护成本极高。此外，目前用CPU实现的DEM方法要么颗粒数较少，要么是低密集颗粒碰撞，无法完全实现高密度大量颗粒碰撞模拟。
使用GPU(图形处理器，GraphicsProcessingUnit)进行通用计算的技术越来越成熟，如当今两大显卡制造商nVIDIA和AMD均支持GPU通用计算。鉴于上述问题，本申请发明人提出了基于GPU的颗粒流动仿真系统和方法。
发明内容
根据本发明，提出一种基于GPU的颗粒流动仿真系统及方法，能够实现高密度颗粒虚拟实验仿真，并降低能耗的同时提高运算效率。
根据本发明一方面，提出了一种基于GPU的颗粒流动仿真方法，其在并行的多个GPU上执行离散单元算法(DEM)方法来进行颗粒流动仿真，该方法包括以下步骤：
a、使用DEM方法对颗粒建模，并将建立的DEM模型分配为多个粒子，将该多个粒子分配给多个计算节点进行处理，每个计算节点的CPU和GPU分别分配有内存空间，并在CPU中进行初始化数据，将初始化后的数据从CPU内存空间拷贝入所述GPU的内存空间中；
b、上述每个计算节点的GPU对每个粒子进行处理，其中每个计算节点的GPU的每个流处理器负责处理一个粒子，并更新粒子存储在GPU内存空间中的坐标及粒子速度；
c、在步骤b的处理过程中，确定每个计算节点控制的粒子，将每个计算节点控制的粒子个数拷贝至CPU内存空间中，并根据GPU内存空间中的粒子数进行动态划分，以按照均衡负载原则动态地确定出每个计算节点计算哪些粒子；
d、利用MPI(MessagePassingInterface)接口协议把上述进行数据动态划分后的粒子在每个计算节点间迁移；
e、根据步骤c得到的每个计算节点控制的粒子，在GPU中计算重叠区域，并将数据拷入CPU内存，而后通过MPI接口协议进行数据交换；
f、每个计算节点的GPU中的每个流处理器根据每个粒子的坐标，计算每个粒子所在GPU内存空间中网格的编号；
g、每个计算节点的GPU中的每个流处理器处理计算每个粒子运动中的受力和加速度；
h、每个计算节点的GPU中的每个流处理器处理每个粒子速度；
i、回到步骤b直至达到指定步数；
j、释放主控节点和计算节点的内存空间。
在一个实施例中，步骤b、步骤f、步骤g和步骤h利用GPU对每个粒子进行并行数据处理，即每个GPU对粒子的处理是同步进行的。
在一个实施例中，步骤d中所述粒子在各节点间迁移利用粒子在节点间传输迁移的方法，即用MPI接口发送和接收函数，实现粒子各物理量的发送和接收，进而实现了粒子在节点间的传输迁移。
在一个实施例中，在步骤e中，所述在GPU中计算重叠区域(Overlap区)包括：利用在GPU计算Overlap区，GPU的一个流处理器处理一个网格。在三维情况下，每个网格有26个网格与之相邻，判断相邻网格是否在当前计算节点中，如果不是，则算作overlap区，从其他节点迁移获得。
根据本发明另一方面，提出了一种基于GPU的颗粒流动仿真方法，包括：
建模步骤，确定颗粒材料、颗粒参数、边界条件和几何体形状，以及颗粒初始分布的区域，并根据预定的颗粒分布区域和数量生成颗粒；
任务管理步骤，根据颗粒总数和多个计算节点上空闲GPU的数目，确定最优GPU数目，并结合最优GPU数目和当前空闲GPU数目确定参与计算的GPU，并将参与计算的GPU状态设置为非空闲；以及
计算步骤，包括
初始化各计算节点的参与计算的GPU，并向各GPU发送计算所需的颗粒信息，
各GPU并行更新预定速度，对接收的颗粒信息排序以生成各自的排序元胞列表，
各GPU并行计算当前各自进程中非零网格编号及网格中颗粒数目，并发送至所述多个计算节点之中的主控节点，由主控节点依照每个GPU最优颗粒数目进行动态划分网格，确定每个GPU并行计算的网格数目和编号，
根据主控节点的确定结果，各GPU并行发送和接收颗粒信息，并在各GPU中重新生成各自的排序元胞列表，
在各GPU中生成当前时刻的碰撞列表，
根据当前时刻的碰撞列表以及前一时刻的碰撞列表和切向相对位移，在各GPU中并行调整切向相对位移的位置，使其与当前碰撞列表保持一致，
根据HM(Hertz-Mindlin)接触力学模型，在各GPU中并行计算每个颗粒的受力和加速度，
存储当前计算结果；
如果计算没有完成，就返回至各GPU并行更新预定速度的步骤，否则计算步骤结束。
在一个实施例中，所述方法还包括展示步骤，包括：确定边界条件，将几何体边界用透明曲面做出；根据颗粒位置和颗粒直径，将颗粒用同色或不同色的小球画出；以及用灰度图显示标量场，并通过将颗粒信息加权映射到网格上，用流线绘制方法画出向量场。
在一个实施例中，保存所有的颗粒的物理信息至外部存储装置。
在一个实施例中，各GPU并行计算相关物理统计量。
在一个实施例中，根据预定的颗粒分布区域和数量生成颗粒包括：在较小的空间内生成若干颗粒，再将这些颗粒平移复制填充至其它空间，直至满足颗粒数量要求。
在一个实施例中，排序元胞列表对所有颗粒按照颗粒所在网格进行排序。
在一个实施例中，采用动态划分方法，由GPU并行计算非零网格编号及网格中粒子数目。
在一个实施例中，在每个GPU中采用一个线程对应一个颗粒的方式来进行计算。
在一个实施例中，计算切向相对位移包括：记录上一时刻的切向相对位移，并根据当前时刻的碰撞列表对其进行更新。
在一个实施例中，使用拷贝或指针交换技术将当前计算结果存储至数组中。
根据本发明又一方面，提出了一种基于GPU的颗粒流动仿真系统，包括：
建模模块，配置为确定颗粒材料、颗粒参数、边界条件和几何体形状，以及颗粒初始分布的区域，并根据预定的颗粒分布区域和数量生成颗粒；
任务管理模块，配置为根据颗粒总数和多个计算节点上空闲GPU的数目，确定最优GPU数目，并结合最优GPU数目和当前空闲GPU数目确定参与计算的GPU，并将参与计算的GPU状态设置为非空闲；以及
计算模块，配置为
初始化各计算节点的参与计算的GPU，并向各GPU发送计算所需的颗粒信息，
各GPU并行更新预定的速度和坐标，对接收的颗粒信息排序以生成各自的排序元胞列表，
各GPU并行计算当前各自进程中非零网格编号及网格中颗粒数目，并发送至主控节点，由主控节点依照每个GPU最优颗粒数目进行动态划分网格，确定每个GPU并行计算的网格数目和编号，
根据主控节点的确定结果，各GPU并行发送和接收颗粒信息，并在各GPU中重新生成各自的排序元胞列表，
在各GPU中生成当前时刻的碰撞列表，根据当前时刻的碰撞列表以及前一时刻的碰撞列表和切向相对位移，在各GPU中并行调整切向相对位移的位置，
根据接触力学模型，在各GPU中并行计算每个颗粒的受力和加速度，
存储当前计算结果，和
如果计算没有完成，就返回至各GPU并行更新预定速度和坐标的步骤，否则计算结束。
在一个实施例中，所述系统还包括展示模块，配置为：确定边界条件，将几何体边界用透明曲面做出；根据颗粒位置和颗粒直径，将颗粒用同色或不同色的小球画出；以及用灰度图显示标量场，并通过将颗粒信息加权映射到网格上，用流线绘制方法画出向量场。
根据本发明再一方面，提出了一种基于GPU的颗粒流动仿真系统，包括：
前端服务器，配置为根据客户端输入的颗粒建模信息生成颗粒信息，并且生成几何体信息；
管理节点，配置为从前端服务器接收颗粒信息和几何体信息，根据颗粒的数目及各个计算节点中空闲的GPU数目，确定使用哪些计算节点中的哪些GPU，然后根据确定的GPU的数目及颗粒在空间中的分布情况确定哪些颗粒由哪个计算节点的哪个GPU进行计算，并根据确定结果进行分配；
多个计算节点，各自包括多个GPU，配置为在多个GPU中并行计算颗粒碰撞导致的每个颗粒的受力，进而计算出加速度，以仿真颗粒流动；
后端服务器，配置为展示仿真结果。
在一个实施例中，前端服务器通过将几何体分解成有限的曲面，并对这些曲面编号，来生成几何体信息。
在一个实施例中，后端服务器在展示的仿真结果中将几何体边界用透明曲面做出，根据颗粒位置和颗粒直径，将颗粒用同色或不同色的小球画出，以及用灰度图显示标量场，并通过将颗粒信息加权映射到网格上，用流线绘制方法画出向量场。
在一个实施例中，前端服务器、管理节点、计算节点和后端服务器通过IB(InfiniBand)网路通信。
根据本发明，实现了基于多个GPU的从建模到结果展示的仿真系统，并结合多GPU的硬件特点，实现了多GPU的颗粒流动仿真方法。根据本发明实施例，利用GPU的超强浮点运算能力、高带宽及多轻量计算核心的特点，充分利用GPU内众多流处理器，将分子动力学中的加速算法合理引入DEM算法，使得DEM算法更适合GPU硬件架构。在多GPU实现时，该算法采用动态划分数据以实现负载均衡的方法，减小了Overlap区及通讯量，极大地提高了GPU和CPU的利用率及运算效率。在可以接受的能耗和时间条件下，取得了非常好的计算效果，达到了在能耗小、维修成本低的同时提高运算效率的效果。
附图说明
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。附图中：
图1是根据本发明实施例的基于GPU的颗粒流动仿真系统的结构示意图；
图2是根据本发明一个实施例的基于GPU的颗粒流动仿真方法的流程图；
图3是根据本发明另一个实施例的基于GPU的颗粒流动仿真系统的模块结构示意图；以及
图4是根据本发明实施例的计算模块的操作流程图。
具体实施方式
以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。
图1是根据本发明实施例的基于GPU的颗粒流动仿真系统的结构示意图。如图1所示，该系统包括前端服务器10、后端服务器20、管理节点30、多个计算节点40-1，...，40N(N是大于1的整数)、IB交换机50和以太网交换机60。此外，图1还示出了该系统包括客户端和存储设备。客户端可以经由互联网与前端服务器10通信，这使得现场实验人员能够远程进行颗粒流动仿真实验。例如，用户可以在客户端上输入建模所需信息或参数，例如颗粒的数目、大小、材料等信息(杨氏模量、泊松(Poisson)比、密度、恢复系数等)及颗粒的分布范围、摩擦系数、边界条件等参数，并给出与颗粒球接触的几何体的材料信息，将这些信息或参数传入前端服务器中。该外部的存储设备可以存储例如各个计算节点的计算结果，以防止死机、断电等意外情况发生导致数据丢失。这里，客户端和外部存储设备是可选的，例如用户可以直接在前端服务器上进行输入，或者计算节点的计算结果可以存储到前端或后端服务器等中。
在图1中，前端服务器10、后端服务器20和计算节点40之间经由IB交换机50连接，并且前端服务器10、管理节点30和计算节点40之间经由以太网交换机60连接。然而，本发明实施例也可以采用其他任何适当的连接方式。在一个实施例中，计算节点40可以是有GPU加速卡的高性能集群。在一个实施例中，每个计算节点都具有GF110核心以上的NVIDIA通用计算卡。在一个实施例中，计算节点使用40Gb的IB(InfiniBand)网路连接。在一个实施例中，前后端服务器各自为一台具有Quadrio6000显卡的图形工作站。例如，工作站内存大于32G，并具有IB网卡。
在本发明实施例的基于GPU的颗粒流动仿真系统中，前端服务器10根据客户端输入的颗粒建模信息生成颗粒信息，并且生成几何体信息。例如，前端服务器10可以接收有关颗粒的尺寸、材料及几何构型的输入，也可通过交互方式增、删颗粒，并移动颗粒的位置。前端服务器10可以通过将几何体分解成有限的曲面，并对这些曲面编号，来生成几何体信息。管理节点30可以任意查看当前各计算节点的运行状态、GPU的工作状态、存储情况等，并可中止提交的任务，保证各任务之间不产生冲突。例如，管理节点30从前端服务器10接收颗粒信息和几何体信息，根据颗粒的数目及各个计算节点中空闲的GPU数目，确定使用哪些计算节点中的哪些GPU，然后根据确定的GPU的数目及颗粒在空间中的分布情况确定哪些颗粒由哪个计算节点的哪个GPU进行计算，并根据确定结果进行分配。各个计算节点40构成整个计算模块，能够处理复杂的边界问题，并行运行多个GPU，具有中断功能(如断电)，可接着中断前状态继续运算，该计算模块使用数据的动态划分方法及指针交换技术，以保证数据的动态平衡。例如，各个计算节点40在各自GPU中并行计算颗粒碰撞导致的每个颗粒的受力，进而计算出加速度，以仿真颗粒流动。后端服务器20展示仿真结果，例如动态展示颗粒当前的构型、温度场、流场、压强场等参数，也可通过交互方式调整观察角度，任意缩放颗粒群组。例如，后端服务器20可以包括显示器等输出设备。后端服务器20可以将几何体边界用透明曲面做出，根据颗粒位置和颗粒直径，将颗粒用同色或不同色的小球画出，以及用灰度图显示温度场等标量场，并通过将颗粒信息加权映射到网格上，用流线绘制方法画出流场、压强场等向量场。
以上系统仅仅是本发明基本构思的一种实现。本领域技术人员可以理解，上述各个部件的功能可以进行再分配或组合，以形成其他的系统构架。此外，如果功能足够强大，上述各个部件的功能可以集成到单个计算机或工作站中。
图2是可以在本发明实施例的仿真系统上执行的基于GPU的颗粒流动仿真方法的流程图。如图2所示，该仿真方法包括如下步骤。
201、使用DEM方法对颗粒建模，并将建立的DEM模型分配为多个粒子，将该多个粒子分配给多个计算节点进行处理，每个计算节点的CPU和GPU分别分配有内存空间，并在CPU中进行初始化数据，将初始化后的数据从CPU内存空间拷贝入所述GPU的内存空间中。
202、上述每个计算节点的GPU对每个粒子进行处理，其中每个计算节点的GPU的每个流处理器负责处理一个粒子，并更新粒子存储在GPU内存空间中的坐标及粒子速度。
203、由于存储在GPU内存空间中的粒子的坐标会发生改变，为保证负载均衡，因此每次计算中各个节点计算的粒子不相同，首先，每个计算节点的GPU计算出该节点控制的粒子，并将每个GPU控制的粒子个数拷贝至CPU内存空间中，并根据GPU内存空间中的网格中粒子数进行数据的动态划分，即按照均衡负载的原则，计算出每个节点计算哪些粒子。
204、利用MPI接口协议把上述进行数据动态划分后的粒子在每个计算节点间迁移。
205、根据步骤203得到的每个计算节点控制的粒子，在GPU中计算重叠区域，并将数据拷入CPU内存，而后通过MPI接口协议进行数据交换。
206、每个计算节点的GPU中的每个流处理器根据每个粒子的坐标，计算出每个粒子所在GPU内存空间中网格的编号。
207、每个计算节点的GPU中的每个流处理器处理计算每个粒子运动中的受力和加速度。
208、每个计算节点的GPU中的每个流处理器处理每个粒子速度。
209、回到步骤202直至达到指定步数，完成DEM方法。
210、释放主控节点和计算节点的内存空间。
其中所述步骤202、步骤206、步骤207和步骤208利用GPU对每个粒子进行并行数据处理，即每个GPU对粒子的处理是同步进行的。
步骤204中所述粒子在各节点间迁移利用粒子在节点间传输迁移的方法，即用MPI接口发送和接收函数，实现粒子各物理量的发送和接收，进而实现了粒子在节点间的传输迁移。接受函数为：MPI_Send()和MPI_Recv()函数。
步骤205中，所述在GPU中计算重叠区域(Overlap区)，利用在GPU计算Overlap区的方法，即GPU的一个流处理器处理一个网格，在三维情况下，每个网格有26个网格与之相邻，接着判断相邻网格是否在当前计算节点中，如果不是，则算作overlap区，则从其他节点迁移获得。
其具体如下：
步骤1，每个计算节点在CPU和GPU中分别开辟内存空间，于CPU中初始化数据，并拷贝入GPU中；
步骤2，
计算节点GPU的每个流处理器处理一个粒子，并行更新一步粒子坐标及半步长粒子速度。包含CUDA的Kernel函数：
在调用的时候，根据CUDA的语法要求，采用如下方式调用：
这两个函数的block和grid均采用一维方式，针对不同粒子数，block和grid的值可以调整，对计算机时有一定影响。
步骤3，在每个计算节点的GPU中，计算出该节点控制的粒子，并拷贝至CPU，根据网格内粒子数进行数据的动态划分。
在计算过程中，由于粒子会在不同节点之间迁移，为避免负载不均的情况，本发明使用动态划分数据的方式，以平衡每个节点的计算量。
初始状态下，假设有M个网格，每个网格中粒子数X相同，M个网格(G0～GM-1)被平均划分为N段，分别交给N个节点(P0～PN-1)处理，则每个节点计算的粒子数为(M/N)*X。经过迭代计算后，各个节点Pi所计算的网格范围内的粒子总数发生改变，可通过调整每个节点计算网格的范围来改变计算粒子的总数。动态数据划分的实现如下：
(1)每个节点维护一个与全局网格数量M相同的int型数组iCellCount，调用CUDA内核函数calcParticleNumPerCel()l计算出每个网格中粒子个数，将其保存到iCellCount中。此时iCellCount中的粒子个数只是局部的，仅仅记录了当前节点计算粒子所在各个网格中的个数。
(2)令PID＝0的节点为ROOT节点，调用MPI归约函数MPI_Reduce()，以加法操作将所有节点iCellCount中的信息收集到ROOT节点下的iGlobalCellCount数组中。此时iGlobalCellCount数组中记录的各个网格中粒子个数是全局的，是所有粒子所在各个网格中的个数。
(3)使用iGlobalCellCount数组对各个节点计算网格范围进行划分，划分采用CPU+GPU方式。划分步骤如下：
根据节点个数N将数组iGlobalCellCount平均分为N段并假设每个节点计算的网格范围相同，将各个节点计算网格范围保存在数组iDividedResult中，初始状态下iDividedResult中各元素值为节点i的范围可根据iDividedResult[i*2]和iDividedResult[i*2+1]得到。
调用CUDA内核函数dReducePerSeg()函数分别求出每个段中的粒子个数，并保存在数组iParticlesCountPerSeg＝{X0，X1，...，XN-1}中。
由CPU根据iDividedResult、iParticlesCountPerSeg和iGlobalCellCount来确定最终的划分结果。首先给定理想情况下每个节点计算粒子的个数iParticlesPerNodeIdeal，并读取iParticlesCountPerSeg[0]的值，若iParticlesCountPerSeg[0]＞iParticlesPerNodeIdeal，说明节点0处理的范围过大，则
重复上述过程直到iParticlesCountPerSeg[0]与iParticlesPerNodeIdeal相等或接近。若iParticlesCountPerSeg[0]＜iParticlesPerNodeIdeal说明节点0处理的范围过小，则对上述过程向相反方向处理。当iParticlesCountPerSeg[0]与iParticlesPerNodeIdeal相等或接近后，则iDividedResult[0]，iDividedResult[0*2+1]为节点0的计算范围。
重复(3)中的过程对所有分段进行处理后，便可得到各个节点处理网格的范围。
(4)ROOT节点调用MPI_BCast()函数将划分结果广播到所有节点。
步骤4：利用MPI接口协议把划分好数据后的粒子在各节点间迁移。
各节点根据网格划分结果iDividedResult确定iSendGridInfo数组和iSendParticlesOffset数组。数组iSendGridInfo和iSendParticlesOffset大小与全局网格数相同，其中iSendGridInfo记录各网格位于哪个节点，而iSendParticlesOffset则记录各网格中第一个粒子所在粒子数组中的位置。
根据链表gridInfo长度确定当前节点要向iSendNodeCount个节点发送粒子，并将发送信息写入数组iSendInfo中。数组iSendInfo长度为iSendNodeCount*3，其中iSendInfo[i*3]为接收粒子节点编号PIDR，iSendInfo[i*3+1]为发送粒子的个数，iSendInfo[i*3+2]为发送节点编号PIDS。
ROOT节点调用MPI_Gatherv()函数将所有节点的iSendInfo数组收集到iGlobalSendInfo数组中，根据iGlobalSendInfo[i*3]的值进行从小到大排序，再调用MPI_Scatterv()函数根据iGlobalSendInfo[i*3]的值将三元组发送到对应的节点。
各节点接收ROOT发送的三元组，将其保存在数组iRecvInfo中并开始发送与接收粒子。
步骤5，根据步骤3得到的每个节点控制的粒子，在GPU中计算Overlap区，并将数据拷入CPU内存，而后通过MPI接口协议进行数据交换。
三维DEM在计算过程中，每个网格需要相邻的26个网格(overlap网格)中的粒子数据，因此在重新划分每个节点的网格计算范围和传送粒子后，各节点必须获取overlap网格以确保计算能正确进行。Overlap交换过程实现如下：
将接收到的粒子保存到粒子数组中，同时将发送走的粒子从粒子数组中移除。对新的粒子数组根据所在网格号从小到大排序，并重新计算iCellCount和iSendParticlesOffset数组。
根据iDividedResult数组所记录的当前节点处理网格范围，计算每一个范围内网格的相邻网格，并确定不在当前节点的相邻网格编号以及所在节点编号。
ROOT节点调用MPI_Gatherv将各节点的iSendInfo数组收集到ROOT节点下的iGlobalSendInfo数组中。根据iGlobalSendInfo[i*3]从小到大排序后调用MPI_Scatterv()函数根据iGlobalSendInfo[i*3]的值将三元组发送到相应的节点。
各节点将ROOT发送来的三元组保存在数组iRecvInfo中，根据iCellCount[iRecvInfo[i*3+1]]确定要向编号为iRecvInfo[i*3+2]的节点发送多少个粒子，同时令iSendGridInfo[iRecvInfo[i*3+1]]＝iRecvInfo[i*3+2]。
采用步骤2中的方法将overlap网格中的粒子发送到指定的节点。
步骤6，在每个计算节点GPU中的一个流处理器处理一个粒子，根据粒子的坐标，计算每个粒子所在网格的编号。
为节省内存空间，网格编号采用按行一维存储，并调用CUDA核函数：
calcHash<<<gridsize，blocksize>>>(ParticleHash，ParticleIndex，
x1，x2，x3，
NumParticles)；
得到粒子所在网格编号ParticleHash。对于计算区域外的粒子，计算其所在网格时，可以将其人为地放入计算区域内某网格中，而不影响计算。
然后对根据Cell-list的要求，用下列kernel根据ParticleHash生成cell-list：
根据上述结果，调用下列kernel函数：
生成每个粒子的邻居表NbrLst。根据新生成的NbrLst，计算出新的切向相对位移U。
步骤7，在每个计算节点的GPU中的一个流处理器处理一个粒子，计算受力和加速度。
根据步骤6中得到的NbrLst以及U，结合粒子的坐标、速度、角速度，根据DEM方法的要求，计算出每个粒子的受力和力矩。根据牛顿第二定律，计算出每个粒子的加速度和角加速度。
步骤8，根据步骤7计算的加速度以及角加速度，更新半步粒子的速度，具体方式与步骤2相同。
步骤9，循环，回到步骤2继续计算，直至满足要求；
步骤10，将GPU设备内存中的需要的数据复制到CPU内存，释放主控节点和计算节点内存空间。
下表一示出了上述仿真方法的执行结果。程序在nVIDIA的GPU上运行不同步数，分别采用不同的block和Thread数量执行。
图3是根据本发明另一个实施例的基于GPU的颗粒流动仿真系统的模块结构示意图。如图3所示，该模块化的仿真系统包括建模模块302、任务管理模块304、计算模块306和展示模块308。结合图1，例如建模模块302可以在前端服务器10中实现，任务管理模块304可以在管理节点30中实现，计算模块306可以由计算节点40的集群实现，展示模块308可以在后端服务器20中实现。然而，这些模块也可以适当方式实现在例如单个或多个计算机中。
建模模块302接收颗粒生成所需的信息，例如即颗粒的数目、大小、材料等信息(杨氏模量、泊松(Poisson)比、密度、恢复系数等)及颗粒的分布范围、摩擦系数、边界条件等参数，并给出与颗粒接触的几何体的材料信息。
建模模块302根据接收的信息，生成需要的颗粒模型(也可简称为颗粒)。为保证所生成的颗粒间没有交叠，或者交叠较小，可使用以下几种方法生成颗粒模型：(1)规则生成法，即在规定范围内生成规则的颗粒，但须加入相对于颗粒半径0.1％～1％的变动；(2)每新生成一个颗粒都与前面所有颗粒进行比较，检测与其是否有交叠，若有，则重新生成，否则视为生成成功；(3)先在较小的空间内使用方法(2)生成若干粒子，再将这些颗粒平移复制填充至其它空间，直至满足颗粒数要求。这提高了颗粒分布的随机性并节省了计算时间。除上述三种方法外，对于颗粒数较少的情况，在确定空间范围后，也可通过交互式方法，使用鼠标点击生成。
颗粒生成后，建模模块302处理几何体信息。将几何体分解成有限的曲面，并对这些曲面编号。接着，将生成的颗粒、几何体及其它材料信息提供给任务管理模块304。
任务管理模块304首先根据传输的颗粒数目及空闲的GPU数目来给当前任务分配节点和GPU，若资源不足，则向用户通知或让用户选择等待还是放弃。在确定好GPU后，将初始的颗粒位置信息存入管理节点30的GPU中，并根据GPU的数目及颗粒在空间中的分布情况确定哪些颗粒由哪个计算节点40的哪个GPU卡进行计算。任务管理模块304向计算模块306传送确定结果以分配给各个计算节点40。
当每个计算节点40得到自己需要的颗粒后，先根据当前的加速度积分1/2步，得到1/2步后的速度，并根据这个速度和当前颗粒坐标值来更新所有颗粒的位置。
更新位置后进行碰撞检测。此时需要将空间划分为若干个网格，在计算任意一个颗粒的受力情况时，只需计算与该颗粒相邻网格内的颗粒是否与之碰撞，若发生碰撞，则将碰撞颗粒放入碰撞列表中，并将碰撞颗粒个数加1。
在计算颗粒球受力时，先将该碰撞颗粒的坐标、速度、角速度信息提取出来，计算受力，之后对所有碰撞颗粒求合力并计算颗粒加速度。对于颗粒周围几何体的受力，先计算颗粒与几何体之间的距离，若该距离小于颗粒半径，则认为该颗粒与几何体发生碰撞，将几何体当成质量无穷大但速度和角速度场为零的颗粒，同样可计算颗粒受到几何体的力
为保证能够中断后重新接着计算，根据实际需要，可以每个小时保存一步计算数据。该计算模块306也可根据需要计算密堆率、平均堆积密度、温度粘滞系数等物理量并存储。计算完成后用户若想对结果进行可视化展示，可将数据发送到展示模块308。
下面结合图4进一步描述计算模块306的操作流程。该实施例中，计算模块306的计算过程可以采用“排序的元胞列表法”。该方法对所有颗粒按照颗粒所在网格进行排序，并充分利用cellStart和CellEnd这两个数组的优势。该方法具有结构简单，实现容易，而不失效率的特点。由此，该方法适合各种高密集度的颗粒碰撞，能解决因颗粒速度大而引起跨节点传输的问题。
描述颗粒的物理量有：坐标pos，速度vel，角速度w，加速度a，角加速度beta，颗粒的切向相对位移U。这些变量均为三维变量。此外，还有颗粒所在网格编号hash、颗粒的永久全局编号pid及临时局部编号index，颗粒的碰撞列表CollideList，以及碰撞的颗粒数量CollideListCnt。
元胞是指上述划分所得的网格，在本文中“元胞”与“网格”含义相同，两者可以互换使用。描述元胞i的变量有：cellStart[i]，cellEnd[i]，cellCount[i]，其中i表示元胞的编号，cellStart[i]表示元胞i的开始颗粒的编号，cellEnd[i]表示元胞i的结束颗粒的编号，cellCount[i]表示元胞i中的颗粒总数。
用于描述进程通讯的二维数组可以称为ParticlesSendToEachNode，第i行第j列的元素[i][j]表示第i个节点向第j个节点发送的颗粒总数。
本发明采用的时间积分算法是速度verlet算法(现有的积分算法，参见例如http://en.wikipedia.org/wiki/Verlet_integration)。
如图4所示，在步骤401，进行初始化，包括开辟GPU和CPU的存储空间，将计算的颗粒信息发送至各计算节点的GPU中。
在步骤402，更新预定的速度和坐标，例如根据加速度(或角加速度)更新1/2步的速度(或角速度)，紧接着根据速度更新颗粒的坐标，如下公式所示。
v ( t + &Delta;t 2 ) &RightArrow; = v ( t ) &RightArrow; + a ( t ) &RightArrow; &Delta;t 2 , r ( t + &Delta;t ) &RightArrow; = r ( t ) &RightArrow; + v ( t + &Delta;t 2 ) &RightArrow; &Delta;t
以上两步都是在每个计算节点的GPU中并行完成的。GPU中每个线程(thread)对应着每个颗粒，达到了GPU的最高效率。
这样便得到了新的坐标，需要计算在新的坐标和新的速度(角速度)下的加速度(角加速度)。
由于颗粒的坐标已经发生改变，有可能原本应由A进程(或GPU)计算的颗粒此时应由B进程进行计算，这样就需要将A进程的该颗粒所有信息全部发送至B进程。
首先在每个计算节点的GPU中计算每个颗粒所在网格的编号Hash。以每个颗粒所在网格的编号Hash和粒子的局部自然编号index做key-value排序。此步用thrust库(现有成熟的库，已集成到cuda中，参见例如http://code.google.com/p/thrust/)完成。根据排序之后的hash，在GPU中并行计算得到每个网格i的cellStart[i]、cellEnd[i]和cellCount[i]，即执行步骤403。
根据排序的index，对颗粒所有的物理量进行排序。
至此，根据颗粒所在网格编号对颗粒所有的物理量进行重新排序，结合每个网格i的cellStart[i]，cellEnd[i]，cellCount[i]称为“排序的元胞列表”。
然后在步骤404进行动态划分。具体地，各计算节点将自己已有颗粒网格及颗粒数目发送至多个计算节点之中的主控节点，即在各计算节点中，若cellCount[i]！＝0，则将i与cellCount[i]发送至主动节点。由主控节点对每个计算节点发送的cellCount[i]进行累加，得到整个空间的cellCount[i]。主控节点根据整个空间的cellCount[i]重新划分每个GPU计算的颗粒。划分原则为：以网格为单位，每个GPU都计算连续的网格，并且网格中的颗粒总数接近每个GPU的平均颗粒数。这样，每个GPU都得到了由于颗粒坐标变化引起的颗粒计算范围。
根据新的计算范围和当前每个GPU的计算范围，发送和接收相关的颗粒信息。为确定每个GPU需要发送和接收的颗粒总数，创建二维数组：ParticlesSendToEachNode。该数组的每一维大小均为进程数目(或GPU数目)。ParticlesSendToEachNode[i][j]的含义是第i个GPU需要向第j个GPU发送的颗粒总数，也就是第j个GPU从第i个GPU接收的颗粒总数，该数组的对角线元素全部为零。对该数组第i行求和，就是第i个GPU发送的颗粒总数；对第j列求和，就是第j个GPU接收的颗粒总数。可以利用cellStart和cellCount作为输入来计算数组ParticlesSendToEachNode。同时计算出SendStart。SendStart也是二维数组，SendStart[i][j]是第i个GPU向第j个GPU发送的首个颗粒在数组中的位置。这样，就可以将要发送的颗粒信息从GPU中取出传入至发送颗粒缓冲区，以备发送。接着，通过对数组的列求和，可确定每个GPU接收的粒子总数，开辟相应的缓冲区，以备接收。通过MPI标准函数中的例如异步发送接收方式MPI_Irecv函数和MPI_Isend函数等发送接收相应颗粒的物理信息，直至等待所有发送接收完成。
将接收的数组直接通过cudaMemcpyHostToDevice函数(已知函数，用于GPU中存储与主存交换数据)，拷贝至GPU中各数组的尾端，并释放发送和接收缓冲区。
此时，对每个GPU需要计算的新颗粒的信息均已得到，但需要考虑新加入的颗粒和发送的颗粒，重新计算“排序的元胞列表”，可得到排序好的物理量数组。
由于每个GPU计算的颗粒不是独立的，即GPU之间有重叠(Overlap)区域，因此在步骤405，根据每个GPU计算的网格号，可计算出该GPU需要的Overlap区域。采用与动态划分类似的方法，每个GPU得到需要的Overlap区域的颗粒物理信息，存放至每个数组的尾部。这样，加上Overlap区域后的物理信息数组虽然不是完全排序，但同一个网格中的颗粒是连续存放的，同时计算每个网格的cellStart和cellEnd。
在步骤406，根据颗粒信息以及cellStart，cellEnd算出当前所有粒子的碰撞列表。其方法是：对于任意一个颗粒i，通过texturememory(纹理内存)先取出其坐标，计算其所在网格编号，扫描其周围包括自己在内的共27个网格中其余所有颗粒，如果其余颗粒与该颗粒的质心距离小于二者半径之和，则将这个颗粒标记为该颗粒的碰撞列表CollideList[i][CollideListCnt[i]]中，并将碰撞列表数目CollideListCnt[i]+1。
切向相对位移仅存在于两个颗粒发生接触之时。为计算当前时刻任意一个颗粒i的受力，需要上一时刻的切向相对位移。存储该切向相对位移的数组U维度大小和CollideList的维度大小一样。U[i][j]存储的是颗粒i与颗粒CollideList[i][j]的切向相对位移。因此，为保证计算结果的正确性，必须在计算颗粒的当前时刻的受力之前，根据当前时刻的CollideList和上一时刻的CollideListOld及UOld对数组U进行重新排列。这一排列过程在GPU中实现。具体地，利用输入上一时刻的碰撞列表CollideListOld、CollideListCnt，与CollideListOld对应的UOld，当前时刻的碰撞列表CoolideList[CollideListCnt]作为输入，调整UOld的顺序，得到当前时刻的数组U。
这样，就得到了正确的计算力需要的所有数组，在步骤407根据HM接触力学模型，计算出每个颗粒的受力。具体地，可以利用坐标pos，速度vel，角速度w，颗粒的切向相对位移U，碰撞列表CollideList[CollideListCnt]，来根据HM接触力学公式计算每个颗粒的加速度a和角加速度beta。
得到新的加速度a(角加速度b)之后，在步骤408再次根据以上速度更新1/2步速度。
至此，完成了计算模块中完整的一步的运算。
保存当前所有颗粒的物理信息的数组，为下一次数组做准备。这里在步骤409可以采用拷贝或指针交换技术。指针交换技术将当前数组与下一次计算的数组的首址进行交换，可以减少数据拷贝所需的较长时间。
在步骤410，判断是否进行外部存储。如果需要，在步骤411可以将所有颗粒的所有的物理信息至外部存储设备，防止断电后需要重新计算的危险。在步骤412，判断是否要进行统计。如果需要，在步骤413计算出相关统计物理量，例如均值、方差等。在步骤414，判断是否满足计算结束条件。例如是否执行了预定次数的计算。如果计算没有完成，就返回至步骤402，否则计算结束，保存结果，释放存储空间。
与国际著名软件lammps(广泛应用的开源软件，可参见http://lammps.sandia.gov/)基于8核CPU的实施相比，本发明的基于GPU(例如TelsaM2090)的仿真方法的运算速度提高了10倍左右。
本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (19)

1.一种基于GPU的颗粒流动仿真方法，其在并行的多个GPU上执行离散单元算法来进行颗粒流动仿真，该方法包括以下步骤：
a、使用离散单元算法对颗粒建模，并将建立的离散单元算法模型分配为多个粒子，将该多个粒子分配给多个计算节点进行处理，每个计算节点的CPU和GPU分别分配有内存空间，并在CPU中进行初始化数据，将初始化后的数据从CPU内存空间拷贝入所述GPU的内存空间中；
b、上述每个计算节点的GPU对每个粒子进行处理，其中每个计算节点的GPU的每个流处理器负责处理一个粒子，并更新粒子存储在GPU内存空间中的坐标及粒子速度；
c、在步骤b的处理过程中，确定每个计算节点控制的粒子，将每个计算节点控制的粒子个数拷贝至CPU内存空间中，并根据GPU内存空间中的粒子数进行动态划分，以按照均衡负载原则动态地确定出每个计算节点计算哪些粒子；
d、利用MPI接口协议把上述进行数据动态划分后的粒子在每个计算节点间迁移；
e、根据步骤c得到的每个计算节点控制的粒子，在GPU中计算重叠区域，并将数据拷入CPU内存，而后通过MPI接口协议进行数据交换；
f、每个计算节点的GPU中的每个流处理器根据每个粒子的坐标，计算每个粒子所在GPU内存空间中网格的编号；
g、每个计算节点的GPU中的每个流处理器处理计算每个粒子运动中的受力和加速度；
h、每个计算节点的GPU中的每个流处理器处理每个粒子速度；
i、回到步骤b直至达到指定步数；
j、释放主控节点和计算节点的内存空间。
2.根据权利要求1所述的方法，其中，步骤b、步骤f、步骤g和步骤h利用GPU对每个粒子进行并行数据处理。
3.根据权利要求1所述的方法，其中，步骤d中所述粒子在每个计算节点间迁移包括：利用MPI接口发送和接收函数实现粒子各物理量的发送和接收。
4.根据权利要求1所述的方法，其中，在步骤e中，所述在GPU中计算重叠区域包括：利用GPU的一个流处理器处理一个网格；
在三维情况下，每个网格有26个网格与之相邻，判断相邻网格是否在当前计算节点中，如果不是，则算作重叠区域，从其他节点迁移获得。
5.一种基于GPU的颗粒流动仿真方法，包括：
建模步骤，确定颗粒材料、颗粒参数、边界条件和几何体形状，以及颗粒初始分布的区域，并根据预定的颗粒分布区域和数量生成颗粒；
任务管理步骤，根据颗粒总数和多个计算节点上空闲GPU的数目，确定最优GPU数目，并结合最优GPU数目和当前空闲GPU数目确定参与计算的GPU，并将参与计算的GPU状态设置为非空闲；以及
计算步骤，包括
初始化各计算节点的参与计算的GPU，并向各GPU发送计算所需的颗粒信息，
各GPU并行更新预定的速度和坐标，对接收的颗粒信息排序以生成各自的排序元胞列表，
各GPU并行计算当前各自进程中非零网格编号及网格中颗粒数目，并发送至主控节点，由主控节点依照每个GPU最优颗粒数目进行动态划分网格，确定每个GPU并行计算的网格数目和编号，
根据主控节点的确定结果，各GPU并行发送和接收颗粒信息，并在各GPU中重新生成各自的排序元胞列表，
在各GPU中生成当前时刻的碰撞列表，根据当前时刻的碰撞列表以及前一时刻的碰撞列表和切向相对位移，在各GPU中并行调整切向相对位移的位置，
根据接触力学模型，在各GPU中并行计算每个颗粒的受力和加速度，
存储当前计算结果，和
如果计算没有完成，就返回至各GPU并行更新预定速度和坐标的步骤，否则计算步骤结束。
6.根据权利要求5所述的方法，还包括展示步骤，包括：确定边界条件，将几何体边界用透明曲面做出；根据颗粒位置和颗粒直径，将颗粒用同色或不同色的小球画出；以及用灰度图显示标量场，并通过将颗粒信息加权映射到网格上，用流线绘制方法画出向量场。
7.根据权利要求5所述的方法，还包括：保存作为计算结果的所有的颗粒信息至外部存储装置。
8.根据权利要求5所述的方法，还包括：各GPU并行计算颗粒相关的物理统计量。
9.根据权利要求5所述的方法，其中根据预定的颗粒分布区域和数量生成颗粒包括：在较小的空间内生成若干颗粒，再将这些颗粒平移复制填充至其它空间，直至满足颗粒数量要求。
10.根据权利要求5所述的方法，还包括：排序元胞列表是对所有颗粒按照颗粒所在网格进行排序的列表。
11.根据权利要求5所述的方法，其中采用动态划分方法，由GPU并行计算非零网格编号及网格中粒子数目。
12.根据权利要求5所述的方法，其中在每个GPU中采用一个线程对应一个颗粒的方式来进行计算。
13.根据权利要求5所述的方法，其中计算切向相对位移包括：记录上一时刻的切向相对位移，并根据当前时刻的碰撞列表对其进行更新。
14.根据权利要求5所述的方法，其中使用拷贝或指针交换技术将当前计算结果存储至数组中。
15.一种基于GPU的颗粒流动仿真系统，包括：
建模模块，配置为确定颗粒材料、颗粒参数、边界条件和几何体形状，以及颗粒初始分布的区域，并根据预定的颗粒分布区域和数量生成颗粒；
任务管理模块，配置为根据颗粒总数和多个计算节点上空闲GPU的数目，确定最优GPU数目，并结合最优GPU数目和当前空闲GPU数目确定参与计算的GPU，并将参与计算的GPU状态设置为非空闲；以及
计算模块，配置为
初始化各计算节点的参与计算的GPU，并向各GPU发送计算所需的颗粒信息，
各GPU并行更新预定的速度和坐标，对接收的颗粒信息排序以生成各自的排序元胞列表，
各GPU并行计算当前各自进程中非零网格编号及网格中颗粒数目，并发送至主控节点，由主控节点依照每个GPU最优颗粒数目进行动态划分网格，确定每个GPU并行计算的网格数目和编号，
根据主控节点的确定结果，各GPU并行发送和接收颗粒信息，并在各GPU中重新生成各自的排序元胞列表，
在各GPU中生成当前时刻的碰撞列表，根据当前时刻的碰撞列表以及前一时刻的碰撞列表和切向相对位移，在各GPU中并行调整切向相对位移的位置，
根据接触力学模型，在各GPU中并行计算每个颗粒的受力和加速度，
存储当前计算结果，和
如果计算没有完成，就返回至各GPU并行更新预定速度和坐标的步骤，否则计算结束。
16.根据权利要求15所述的系统，还包括展示模块，配置为：确定边界条件，将几何体边界用透明曲面做出；根据颗粒位置和颗粒直径，将颗粒用同色或不同色的小球画出；以及用灰度图显示标量场，并通过将颗粒信息加权映射到网格上，用流线绘制方法画出向量场。
17.一种基于GPU的颗粒流动仿真系统，包括：
前端服务器，配置为根据客户端输入的颗粒建模信息生成颗粒信息，并且生成几何体信息；
管理节点，配置为从前端服务器接收颗粒信息和几何体信息，根据颗粒的数目及各个计算节点中空闲的GPU数目，确定使用哪些计算节点中的哪些GPU，然后根据确定的GPU的数目及颗粒在空间中的分布情况确定哪些颗粒由哪个计算节点的哪个GPU进行计算，并根据确定结果进行分配；
多个计算节点，各自包括多个GPU，配置为在多个GPU中并行计算颗粒碰撞导致的每个颗粒的受力，进而计算出加速度，以仿真颗粒流动；
后端服务器，配置为展示仿真结果；
其中后端服务器在展示的仿真结果中将几何体边界用透明曲面做出，根据颗粒位置和颗粒直径，将颗粒用同色或不同色的小球画出，以及用灰度图显示标量场，并通过将颗粒信息加权映射到网格上，用流线绘制方法画出向量场。
18.根据权利要求17所述的系统，其中前端服务器通过将几何体分解成有限的曲面，并对这些曲面编号，来生成几何体信息。
19.根据权利要求17所述的系统，其中前端服务器、管理节点、计算节点和后端服务器通过InfiniBand网路通信。
CN201310195595.6A 2012-12-20 2013-05-22 颗粒流动仿真系统和方法 Active CN103324780B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201310195595.6A CN103324780B (zh) 2012-12-20 2013-05-22 颗粒流动仿真系统和方法

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
CN201210558134.6 2012-12-20
CN201210558134 2012-12-20
CN2012105581346 2012-12-20
CN201310195595.6A CN103324780B (zh) 2012-12-20 2013-05-22 颗粒流动仿真系统和方法

Publications (2)

Publication Number Publication Date
CN103324780A CN103324780A (zh) 2013-09-25
CN103324780B true CN103324780B (zh) 2016-03-16

Family

ID=49193522

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201310195595.6A Active CN103324780B (zh) 2012-12-20 2013-05-22 颗粒流动仿真系统和方法

Country Status (5)

Country Link
US (1) US10007742B2 (zh)
JP (1) JP6009075B2 (zh)
CN (1) CN103324780B (zh)
GB (1) GB2523640B (zh)
WO (1) WO2014094410A1 (zh)

Families Citing this family (52)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103617088B (zh) * 2013-11-29 2018-07-24 深圳中微电科技有限公司 在不同类型线程中分配内核资源的方法、装置及其处理器
US10210287B2 (en) 2013-12-31 2019-02-19 Disney Enterprises, Inc. Augmented material point method for simulating phase changes and varied materials
US10019827B2 (en) * 2013-12-31 2018-07-10 Disney Enterprises, Inc. Material point method for simulation of granular materials
WO2015108980A1 (en) * 2014-01-17 2015-07-23 Conocophillips Company Advanced parallel "many-core" framework for reservoir simulation
CN105389855B (zh) * 2014-08-26 2019-11-01 三星电子株式会社 对对象进行建模的方法和设备
JP6547547B2 (ja) * 2015-09-25 2019-07-24 富士通株式会社 粒子シミュレーションプログラム、計算機資源配分方法、および粒子シミュレーション装置
JP6543557B2 (ja) * 2015-11-11 2019-07-10 富士通株式会社 粒子シミュレーションプログラム、粒子シミュレーション装置、及び計算機資源配分方法
CN107016180A (zh) * 2017-03-30 2017-08-04 中国石油大学(北京) 一种颗粒流动仿真方法
US10262390B1 (en) * 2017-04-14 2019-04-16 EMC IP Holding Company LLC Managing access to a resource pool of graphics processing units under fine grain control
US10275851B1 (en) 2017-04-25 2019-04-30 EMC IP Holding Company LLC Checkpointing for GPU-as-a-service in cloud computing environment
CN107230242B (zh) * 2017-06-07 2020-09-25 广州酷狗计算机科技有限公司 粒子映射方法和装置
US10325343B1 (en) 2017-08-04 2019-06-18 EMC IP Holding Company LLC Topology aware grouping and provisioning of GPU resources in GPU-as-a-Service platform
CN107766148B (zh) * 2017-08-31 2021-02-19 北京百度网讯科技有限公司 一种异构集群及任务处理方法和装置
KR101980699B1 (ko) * 2017-10-31 2019-05-22 한국과학기술원 공간 데이터를 분산 처리하는 시스템 및 방법
US10698766B2 (en) 2018-04-18 2020-06-30 EMC IP Holding Company LLC Optimization of checkpoint operations for deep learning computing
CN110795228B (zh) 2018-08-03 2023-08-25 伊姆西Ip控股有限责任公司 用于训练深度学习模型的方法和制品、以及计算系统
KR102181989B1 (ko) 2018-08-27 2020-11-23 이에이트 주식회사 복수의 프로세서를 이용한 입자 기반의 유체 시뮬레이션 방법 및 유체 시뮬레이션 장치
US10776164B2 (en) 2018-11-30 2020-09-15 EMC IP Holding Company LLC Dynamic composition of data pipeline in accelerator-as-a-service computing environment
CN109753726A (zh) * 2019-01-04 2019-05-14 东南大学 一种基于边界盒搜索方法和gpu的球磨机介质运动仿真方法
CN109977505B (zh) * 2019-03-13 2024-02-09 南京大学 基于gpu矩阵计算的离散元孔隙系统快速搜索方法
CN110275732B (zh) * 2019-05-28 2023-02-21 上海交通大学 在ARMv8处理器上质点网格法的并行实现方法
WO2020254066A1 (en) * 2019-06-20 2020-12-24 Asml Netherlands B.V. Method for patterning process modelling
JP7244757B2 (ja) 2019-07-01 2023-03-23 富士通株式会社 情報処理装置、粒子シミュレーション方法および粒子シミュレーションシステム
CN110348156B (zh) * 2019-07-18 2022-10-14 河南理工大学 一种炉料颗粒在高炉旋转溜槽内部运动的仿真方法
CN112035995A (zh) * 2019-07-19 2020-12-04 交通运输部天津水运工程科学研究所 基于gpu计算技术的非结构网格潮汐潮流数值模拟方法
CN112528456B (zh) * 2019-09-18 2024-05-07 曙光信息产业(北京)有限公司 一种异构节点计算系统及方法
CN110766709B (zh) * 2019-10-30 2021-03-30 成都理工大学 基于图像识别的滑坡颗粒堆积特征识别方法
CN111222262B (zh) * 2019-10-30 2023-09-29 中国中元国际工程有限公司 基于质量比影响的气浮隔振平台性能优化设计方法
CN110992456B (zh) * 2019-11-19 2021-09-07 浙江大学 一种基于位置动力学的雪崩模拟方法
CN111027244B (zh) * 2019-12-03 2024-03-12 天津大学 一种百亿级颗粒模型的构建方法
KR102371345B1 (ko) * 2020-05-07 2022-03-07 국방과학연구소 도심지 유동해석 방법 및 장치
CN112001108B (zh) * 2020-07-08 2024-02-02 中国人民解放军战略支援部队信息工程大学 锥束ct蒙特卡洛仿真集群并行加速方法及系统
CN112100939B (zh) * 2020-09-14 2023-06-16 福建天晴在线互动科技有限公司 一种基于Compute Shader的实时流体仿真方法及其系统
CN112380788B (zh) * 2020-11-06 2022-03-01 天津大学 一种超椭球颗粒与流场双向耦合的半解析计算方法
CN113177346B (zh) * 2020-11-17 2022-06-10 西北工业大学 一种锅炉煤粉运输安全性判断方法及系统
CN112380793B (zh) * 2020-11-18 2024-02-13 上海交通大学 基于gpu的湍流燃烧数值模拟并行加速实现方法
CN112464543B (zh) * 2021-01-28 2021-04-06 北京科技大学 一种计算vim冶炼过程中的夹杂物运动的方法
CN113283048A (zh) * 2021-03-12 2021-08-20 南京航空航天大学 一种多枝晶运动相场法并行模拟的碰撞检测和合并方法
CN113221200B (zh) * 2021-04-15 2022-10-25 哈尔滨工程大学 一种适用于堆芯颗粒分布不确定性分析的三维高效随机排布方法
CN113239522B (zh) * 2021-04-20 2022-06-28 四川大学 一种基于计算机集群的大气污染物扩散模拟方法
CN113282976B (zh) * 2021-04-30 2023-04-11 重庆大学 基于comsol的粉床构建方法
CN112948643B (zh) * 2021-05-13 2021-08-06 中国空气动力研究与发展中心计算空气动力研究所 一种基于线程并行的结构化网格流线积分方法
CN113343605B (zh) * 2021-06-29 2022-03-25 华南理工大学 一种汽车涂层抗石击性标准实验的仿真对标方法
CN114254572B (zh) * 2021-12-16 2024-01-02 西北工业大学太仓长三角研究院 考虑污染物沉积的航发压气机流场性能预测方法及系统
CN114186468B (zh) * 2021-12-16 2024-03-29 西安交通大学 一种基于gpu加速的二次电子发射模拟方法及系统及设备
CN114491864B (zh) * 2022-01-26 2022-12-13 哈尔滨工程大学 一种具有参数化、可重构特征的核动力管网模型预处理方法
CN115600449A (zh) * 2022-07-05 2023-01-13 浙江大学（Cn） 基于细长柔性物件组成的大规模颗粒系统的模拟预测方法
CN115719047B (zh) * 2022-11-14 2023-06-06 沐曦集成电路(上海)有限公司 基于波形gpu联合仿真系统
CN116738806A (zh) * 2023-04-18 2023-09-12 西北农林科技大学 含丁坝水流交汇区微塑料输运规律的模拟方法
CN116562065B (zh) * 2023-07-12 2023-09-12 北京凌云智擎软件有限公司 一种网格拓扑的转换方法、设备及装置
CN117272769A (zh) * 2023-09-27 2023-12-22 长江大学 一种颗粒混合过程仿真数据的处理方法
CN117217132B (zh) * 2023-11-08 2024-02-09 四川新能源汽车创新中心有限公司 一种基于cfd-dem耦合的颗粒运动模拟方法

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7952583B2 (en) * 2000-06-19 2011-05-31 Mental Images Gmbh Quasi-monte carlo light transport simulation by efficient ray tracing
JP5371221B2 (ja) * 2007-09-11 2013-12-18 プロメテック・ソフトウェア株式会社 粒子法シミュレーションのためのスライスデータ構造、およびスライスデータ構造を利用した粒子法シミュレーションのｇｐｕへの実装方法
JP5124845B2 (ja) 2007-09-14 2013-01-23 本田技研工業株式会社 金属リング周長補正方法
CN101727653B (zh) * 2008-10-31 2012-03-07 中国科学院过程工程研究所 一种基于图形处理器的多组分系统离散模拟计算方法
JP5408611B2 (ja) * 2009-03-31 2014-02-05 独立行政法人海洋研究開発機構 粒子シミュレーション装置及び粒子シミュレーション方法
US8531463B2 (en) 2009-08-10 2013-09-10 Dem Solutions Limited Method and apparatus for discrete element modeling with a virtual geometry object
CA2848319C (en) * 2010-09-15 2020-12-29 Commonwealth Scientific And Industrial Research Organisation Discrete element method
US9038088B2 (en) * 2011-03-10 2015-05-19 Nec Laboratories America, Inc. Load balancing on hetrogenous processing cluster based on exceeded load imbalance factor threshold determined by total completion time of multiple processing phases
US8554527B2 (en) * 2011-04-08 2013-10-08 Japan Agency For Marine-Earth Science And Technology Particle simulator and method of simulating particles
US9760376B1 (en) * 2016-02-01 2017-09-12 Sas Institute Inc. Compilation for node device GPU-based parallel processing

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
GPU中的流体场景实时模拟算法;陈曦 等;《计算机辅助设计与图形学学报》;20100331;第22卷(第3期);396-405 *
Modelling the solids in-flow and solids conveying of single-screw extruders using the discrete element method;MOYSEY P A 等;《Powder Technology》;20050516;第153卷(第2期);95-107 *
基于CUDA的离散粒子系统模拟仿真及其实现;江顺亮 等;《南昌大学学报(工科版)》;20110930;第33卷(第3期);290-294 *

Also Published As

Publication number Publication date
WO2014094410A1 (zh) 2014-06-26
CN103324780A (zh) 2013-09-25
GB2523640A (en) 2015-09-02
US10007742B2 (en) 2018-06-26
US20150213163A1 (en) 2015-07-30
GB2523640B (en) 2020-05-27
GB201500658D0 (en) 2015-03-04
JP6009075B2 (ja) 2016-10-19
JP2015530636A (ja) 2015-10-15

Similar Documents

Publication Publication Date Title
CN103324780B (zh) 颗粒流动仿真系统和方法
Karloff et al. A model of computation for MapReduce
US9529590B2 (en) Processor for large graph algorithm computations and matrix operations
CN107016180A (zh) 一种颗粒流动仿真方法
CN110929456B (zh) 移动粒子法并行计算等效粒子负载均衡加速方法
Chaudhary et al. Cloudy GSA for load scheduling in cloud computing
MacFarland et al. A new parallel P3M code for very large-scale cosmological simulations
CN103049241A (zh) 一种提高cpu+gpu异构装置计算性能的方法
CN103440163A (zh) 使用gpu并行实现的基于pic模型的加速器仿真方法
CN102411658A (zh) 一种基于cpu和gpu协作的分子动力学加速方法
CN114490011A (zh) N体模拟在异构架构的并行加速实现方法
Guan et al. pRPL 2.0: Improving the parallel raster processing library
Choquette Nvidia hopper gpu: Scaling performance
Zhang et al. An fast simulation tool for fluid animation in VR application based on GPUs
Chien et al. Dynamic load balancing on a network of workstations for solving computational fluid dynamics problems
CN100492371C (zh) 分布型cad装置
CN113703955A (zh) 计算系统中数据同步的方法及计算节点
CN105531602A (zh) 利用多个加速处理部件(apc)实现时域有限差分模型的系统和方法
EP2141594A2 (en) Vectorized parallel collision detection pipeline
Cambier et al. A task-based distributed parallel sparsified nested dissection algorithm
Plimpton et al. Rendezvous algorithms for large-scale modeling and simulation
Murphy Computational social science and high performance computing: A case study of a simple model at large scales
Yong et al. GVM based intuitive simulation web application for collision detection
Boyle et al. Performance and Energy Simulation of Spiking Neuromorphic Architectures for Fast Exploration
CN105389220B (zh) 在交互式r语言平台中进行并行线性代数计算的方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant