工程精度离散单元法计算程序优化¶
一、题目背景¶
离散单元法 (discrete element method, DEM) 由岩土工程学者 Cundall 和 Strack 于 1979 年提出，是一种数值仿真方法。它将世界抽象为若干细小颗粒 (particle) 的集合，将时间抽象为若干个细小的时间步 (time step)，在每个时间步内，分析颗粒之间的接触 (contact)，通过适当的接触模型 (contact model) 计算接触力 (矩)，利用牛顿第三定律 (Newton's third law) 赋予颗粒上。归集每个颗粒所受的合力 (矩)，然后利用牛顿第二定律 (Newton's second law) 求得颗粒 (角) 加速度，再对 (角) 加速度两次时间积分，得到颗粒位置 (方位) 的变动。重复以上过程，即可模拟时间的流逝和世界的演变。
DEM 因为基本假设简单，适用范围广，限制少。相比有限单元法 (finite element method, FEM)，DEM 可以自然模拟开裂、破碎等大变形问题，显式 DEM 也没有迭代收敛问题。DEM 被广泛用于土木工程、化学工程、工业工程等领域，通过数值仿真预测产品性能或工艺效果，从而指导工程实践。DEM 具有数值仿真共同的低成本优势，能以相比实验更低的成本，创造显著的经济效益。
然而，DEM 有一个显著的劣势——计算量大，计算速度慢。因此，高性能计算 (high performance computing, HPC) 对 DEM 仿真的加速必不可少，以保证 DEM 的计算时间处于可接受的范围内，确保生产效率。
2022 年 09 月，Taichi 社区举办了活动“Taichi DEM 代码优化挑战赛”。参赛作品使用 Taichi 高性能计算语言，编写 DEM 程序，丰富和完善其 DEM 各项功能。
您将获得上述比赛优胜作品代码，作为本赛题的内容之一。该作品实现了一个复杂的 DEM 接触模型——胶结模型 (bond model)，胶结模型的特点是，在仿真一开始，通过外部命令控制颗粒之间按照一定判断准则形成胶结接触 (bond contact)。在此之后的仿真过程中，胶结可以进行受压、受拉、受剪等一系列力学响应，但当其正应力或剪应力满足破坏条件时，胶结将破坏消失，不再恢复。除此之外，颗粒之间的天然碰撞，采用不同于胶结接触的 Hertz-Mindlin 接触模型。
作品的一个代表 Demo 如下，图中，一个斯坦福兔子形状的胶结团块自由下落撞击地面并摔碎。团块摔碎，伴随着其中颗粒间的胶结的破坏 （图片较大，请耐心等待。下载附件 D 以获得较好播放体验。
二、赛题描述¶
您将通过以下链接获取上述作品代码和 Demo：
https://github.com/Denver-Pilphis/taichi_dem/tree/question/PKU_HPC_Competition
其中，dem.py 文件包含了所有的程序代码，包含了详细的注释。运行 dem.py，将会在同目录下写出 output.p4p 和 output.p4c 两个文件，作为输出数据文件，分别存储每个存档时间点的颗粒信息和颗粒接触信息。
对于使用 Taichi 实验的团队，具体任务是：
- 请设计测试方法，测试程序性能瓶颈位置，并将其找出。
- 设法对至少一处性能瓶颈进行优化。
注意：
- 只优化上述链接内提供的程序的代码，不优化 Taichi Runtime。
- 本次比赛测试环境限定为 Nvidia CUDA。
选手可以使用 CUDA 重新实现本程序，或作为组件调用，但需提交完整的依赖、编译脚本和运行说明。DEM 的实现需为选手原创。
三、提示¶
-
第二部分给出的 GitHub 仓库中，程序代码对应的文档请见附件 A: DEM 理论部分文档 和 附件 B: DEM HPC 实现文档
-
可以使用 Taichi kernel profiling 测试，输出各个 kernel 的耗时。
四、测试方法¶
第二部分给出的 GitHub 代码仓库中的 README.md 提供了 6 个 Example。6 个 Example 的输入数据文件在 Resource 目录下。dem.py 开头的注释详细说明了每个 Example 对应的输入数据文件和对应的物理参数 (注释中标注的物理参数，为代码中给定的默认参数列表中需要改变的部分; 未标注的参数，应保持原始值)。除 Example 1 外，所有的 Example 都可以用于赛中测试使用，但需要注意的是，DEM 的物理参数对计算结果有显著影响 (评分要求计算结果与原有代码结果一致，详见第五部分)，因此，在变更输入数据文件时，需要同步准确变更对应的物理参数，以确保计算结果正确。Example 2 为包含胶结的最小规模算例。
最终的测试采用第 6 个 Example: 斯坦福兔子自由落体撞击水平表面 (快速计算版)。程序运行结束后，抽取 output.p4p 和 output.p4c 中任意若干个存档时间点，校验颗粒与接触信息 (均为双精度浮点数) 与原始代码运行结果是否一致。双精度浮点数一致的判据如下：
abs(x′−x)<tol=10e−12
式中 x′为优化后代码输出的双精度浮点数，x 为原始代码输出的双精度浮点数，tol 为双精度浮点容差。
五、提交说明¶
选手需提交一个压缩包，包括代码和解题文档。
代码部分¶
如使用 Taichi 实现赛题，可保留原目录架构。如使用 CUDA 实现，需提供编译说明。
解题文档¶
一个文字版 pdf 文件，命名为writeup.pdf。开篇需声明本文为（团队名）第零届 PKU HPCGame DEM 赛题 writeup，此文由本团队原创。
文档中应简要说明思路和优化策略，文档长度不作为评分依据，请尽量精简。
六、评分说明¶
选手得分为性能评测和综合评测两部分。综合评测主要考察选手的代码质量和解题文档质量，除非质量不可接受，否则不作为区分手段。
性能评测部分，以时间作为评分标准。得分不设上限。最终奖金的瓜分规则与比赛通用规则一致。注意，最终测试时，测试输出结果与原始代码运行结果不一致，得分为零分。
Taichi¶
对于使用Taichi 的团队，dem.py 文件运行时，默认会输出计算至存档点的耗时、存档 IO 耗时；在计算完成后，会输出总耗时。第四部分最终测试时，以最后输出的总耗时为结果。组委会将重点核查选手程序的计时模块和正确性核查模块，请勿随意更改。
CUDA¶
对于使用CUDA 的团队，程序需要从time.txt 读取时间并输出结果，我们将核验结果的正确性。
输入格式：¶
10 个浮点数，代表着需要输出结果的时刻。例：
（仅作为参考，实际测试中并不相同）
1.2345678
1.2345678
1.2345678
1.2345678
1.2345678
1.2345678
1.2345678
1.2345678
1.2345678
1.2345678
输出格式：¶
输出到文件，命名为outputX.p4p 和outputX.p4p,X 测试点，编号是 1～10。如output8.p4p 和output8.p4p 将采用程序运行时间−CUDA 初始化时间平均值−IO 时间最优值为评分的时间标准。
附件¶
附件 A: DEM 理论部分文档
下载DEM.pdf
附件 B: DEM HPC 实现文档
下载HPC.pdf