截止 10/28 共收到 9 份 DEM 代码优化的投稿。最佳人气奖的投票正式开启！请大家为自己喜欢的作品投上一票吧～
近期不少同学提到对于 DEM、PBF 等方法中的快速碰撞检测很感兴趣，于是我和包乾写了一篇知乎文章和示例代码，希望能够讲清楚一种常用算法：用 39 行 Taichi 代码加速 GPU 粒子碰撞检测 - 知乎 12
本文中的 DEM 算例代码将公布在：GitHub - taichi-dev/taichi_dem: A minimal DEM simulation demo written in Taichi 16
DEM 是非常有趣的算法，稍加修改就能做出各种效果。欢迎大家在我们的模板仓库的基础之上进一步开发，优化算法或做出有趣的效果！我们将在所有参赛作品中评选出优化效果最好的 3 个作品送上 Airpods、胶囊咖啡机和键盘，还有更多社区参与奖等你来拿~
提交作品的小伙伴都可以添加小助手的微信（微信号：yanqingdw）进入交流群，围观更多有趣的创意！
参赛时间：9/18-10/28
在 GitHub 上使用 "use this template"从给定的 Template repo 中复制出自己的新 repo 并公开 ，在此基础上形成作品 repo，需要包含源代码 dem.py 文件、README.md 文件和 requiremements.txt。
在 dem.py 中添加你的代码，加速或者优化 DEM，实现知乎文末提到的各种算法拓展。
提交作品：请在本帖评论区提交作品。提交内容为：
简单介绍自己的作品优化了 DEM 哪个部分，如果有原理说明更佳~
附上作品的 GitHub Repo 地址
附上作品的 GIF 动图
可参考 @llinus 这位同学的投稿：
注：如果作品没有对算法进行拓展和优化，而仅仅是调整粒子配置实现视觉效果展示的也欢迎在评论区提交，我们也会送出 Taichi 周边套装~纯视觉展示不计入最后代码优化评选。
本次的代码优化挑战赛规则很简单，没有代码行数限制。仅有以下两个要求：
仅可改动 dem.py 完成代码优化。
不可以 import 其他库
本次大赛设置一二三等奖和最佳人气奖。所有奖项都将在 10/31-11/4 期间开启评审&投票。最佳人气奖投票通道将在 10/31 日开启~
还没有思路的同学，可以看看以下几个优化方向哦~
事实上我们可以将邻域搜索的范围从 3 x 3 进一步缩小一半左右，你能想想怎么构建这个搜索算法吗？
现在的算例中为了简化计算细节，并没有考虑 DEM 中常用到的粒子的角动量和不规则形状，你能扩展本算例并加入这些新的 feature 吗？
其实 Taichi 对于 GPU 上的动态数据结构有很好的支持，如 taichi_elements 中就采用了 dynamic SNode 来维护每个网格 cell 中的粒子下标数组。你能尝试用 Taichi 的dynamic SNode 来实现本文的功能吗？
其实 Position-based fluid (PBF) 也可以用这种方法进行加速，你可以试一试吗？
其实文中提到的并行计算方式完全可以用 ti.atomic_add 去掉所有串行区域，进一步并行化，你觉得应该咋搞？
本文的代码比较简单，完全是 2D 的，你能把它扩展成 3D 吗？
4 个赞
llinus
2022 年9 月 18 日 16:40
#4
我也来抛个砖，做了个搅拌机的demo，这是内置于商业离散元软件PFC中的案例，GPU的计算效率相比之下是真的高
将来会慢慢实现更符合真实物理的接触本构模型和三角网格描述的边界，以求得和PFC相当的计算精度。
代码仓库：GitHub - Linus-Civil/GeoBlender: A 3D demo of blender in civil/geotechnical/geological engineering based on Taichi DEM
5 个赞
大佬请问有详细的算法描述吗，主要是颗粒和中间的叶片作用是怎么实现的？
llinus
2022 年9 月 19 日 14:18
#6
我这就是个玩具Demo，代码里只实现了法向力。你如果对工程级别的力学本构细节感兴趣，推荐你去看看《计算颗粒力学及工程应用（季顺迎）》。
请问为啥case2 输出的结果里粒子多了很多不存在的碰撞？
for i, j in ti.ndrange(grid_n, grid_n):
for k in range(grid_n):
ti.atomic_add(_column_sum_cur[i,j], grain_count[i, j, k])
linear_idx = i * grid_n * grid_n + j * grid_n + k
list_head[linear_idx] = _column_sum_cur[i,j]- grain_count[i, j, k]
list_cur[linear_idx] = list_head[linear_idx]
list_tail[linear_idx] = _column_sum_cur[i,j]
for i, j, k in ti.ndrange(grid_n, grid_n, grid_n):
ti.atomic_add(_column_sum_cur[i,j], grain_count[i, j, k])
linear_idx = i * grid_n * grid_n + j * grid_n + k
list_head[linear_idx] = _column_sum_cur[i,j]- grain_count[i, j, k]
list_cur[linear_idx] = list_head[linear_idx]
list_tail[linear_idx] = _column_sum_cur[i,j]
代码：
python3 dem-3d.py 即可复现
这样也可以 ， 看来我还是得看看 cpp11
# case 3 test okay
for i, j, k in ti.ndrange(grid_n, grid_n, grid_n):
# we cannot visit prefix_sum[i,j] in this loop
pre = ti.atomic_add(prefix_sum[i,j], grain_count[i, j, k])
linear_idx = i * grid_n * grid_n + j * grid_n + k
list_head[linear_idx] = pre
list_cur[linear_idx] = list_head[linear_idx]
# only pre pointer is useable
list_tail[linear_idx] = pre + grain_count[i, j, k]
Otis
2022 年9 月 29 日 13:11
#14
请问在你的代码中折半邻居搜索后计算效率会有所提升嘛，我试过折半搜索，同时去掉了相邻网格时id1要小于id2的条件，照理来说减少搜索次数速度应该提升，但发现效率反而下降了
《偷功》
实现了一个简单旋转小球的效果
效果：
PBF 版本： 增加了不可压缩性，粘度等， 这个实现参考了 Ten Minute Physics 中的实现
代码：GitHub - mrzhuzhe/taichi_dem: A minimal DEM simulation demo written in Taichi.
python3 dem-3d.py # 黄色 dem 旋转小球
python3 pbf.py # 蓝色 pbf 旋转小球
python3 dem-3d-snode.py # 黄色 dem 旋转小球，但是临域搜索的链表结构是用dynamic snode实现
已实现：
事实上我们可以将邻域搜索的范围从 3 x 3 进一步缩小一半左右，你能想想怎么构建这个搜索算法吗？
其实 Taichi 对于 GPU 上的动态数据结构有很好的支持，如 taichi_elements 中就采用了 dynamic SNode 来维护每个网格 cell 中的粒子下标数组。你能尝试用 Taichi 的dynamic SNode 来实现本文的功能吗？
其实 Position-based fluid (PBF) 也可以用这种方法进行加速，你可以试一试吗？
其实文中提到的并行计算方式完全可以用 ti.atomic_add 去掉所有串行区域，进一步并行化，你觉得应该咋搞？
后续问题：
在GPU3090上 黄色 dem 旋转小球，如果临域搜索的数据结构是用dynamic snode 只有 30-40 fps
3d的情况下粒子靠近边界时会卡在一起，会产生数值下溢，目前只能随机把小球往回移动一小段距离，这个方法有点过于trick，有没有通用的解决办法？
if (pos[i][0] >= maxX):
pos[i][0] = maxX - eps * ti.random();
fp32 在小步长下很容易数值溢出，是否有通用的方法解决？
3 个赞
Mugi
2022 年10 月 2 日 14:35
#17
岩土专业，代码小白。在实现过程中遇到一点问题，想要请教社区的大家
代码在这：GitHub - chenqing253/taichiDEM
尝试实现加入tangential force and rotation, 遇到了以下几点问题：
如果根据以下公式来进行 切向力的更新\Delta F^s = -k^s * \Delta U^s , 其中 斜向分量 \Delta U^s is given by relative velocity V^s ：
t_i is Tangential unit vector. 上述就是求两个颗粒在接触切平面上的相对滑移来求摩擦力
可以看出切向力可能需要在前一步的接触上的切向力进行更新，有什么方法或者数据结构可以记录 前一步的接触吗 ？
我忽略了需要在前一步的切向力的基础上更新这一条件，尝试写了一下，结果不正确（我觉得比较合理的结果是形成一个沙丘状）。如图~
我用红色和蓝色标注了不同旋转方向的颗粒，旋转了但是低于平均值的为白色。而且颗粒模拟过程中比较Q弹，我也想不出哪里出错了，是因为一定要在前一步的切向力上跟新吗？
2 个赞
Otis
2022 年10 月 3 日 11:33
#18
赵仕威老师有篇文章里有写A thread-block-wise computational framework for large-scale hierarchical continuum-discrete modeling of granular media，采用多个列表存储数据
1 个赞
Mugi
2022 年10 月 10 日 15:48
#19
我也来分享一个小的demo~，我在代码中加入了切向力，即摩擦力，实现方法比较简单粗暴。为了实现静摩擦，添加了两个integer类型的ti.field容器用于‘记忆’上一个时间步和现在的接触，接触用两个接触颗粒的id的数对表示，还有两个float类型的ti.field 容器用于储存相应的切向力。Field容器的形状为nx10, n 为颗粒的个数，10是考虑到这个颗粒周围的接触不会超过10个。求算接触 (i,j) 的力时，会遍历第i行的容器。（本来想直接在class里加一个list的…)
第一次尝试用taichi语言，实现过程中踩了不少坑，例如调参数总是把试样调爆炸 orz 。希望随着不断学习taichi的过程中能实现得更漂亮简洁一些，像用 n \times n 的稀疏数组来实现对接触的记录。
代码仓库：A small demo implement shear force based on Taichi DEM
代码效果：
没加入之前总是会变得比较平：
这里没有放gif，因为我调小了重力整个模拟过程很慢，一点也不炫酷orz，否则加入或不加入切向力的效果是一样的，我尝试调过stiffness，但是会‘爆炸’，以后还得再调调参数。
参考文献：
Bell, N., et al. (2005). Particle-based simulation of granular materials. Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation. Los Angeles, California, Association for Computing Machinery: 77–86.
Lee, J. and H. J. Herrmann (1993). “Angle of repose and angle of marginal stability: molecular dynamics of granular particles.” Journal of Physics A: Mathematical and General 26(2): 373-383.
3 个赞
从工程角度使用太极语言完整实现离散单元法。
使用BIMBase ™二次开发可视化。BIMBase ™是北京构力科技有限公司研发的用于BIM的图形平台。BIMBase开发平台-BIMBase Platform
使用Houdini ™可视化。Houdini ™是一款特效软件，工作流程灵活，对用户友好。Houdini | 3D Procedural Software for Film, TV & Gamedev | SideFX
因Taichi论坛上传图片大小限制，上述Demo GIF效果欠佳。请前往GitHub代码仓库查阅高清GIF (地址下述)。
Denver Pilphis (Di Peng) - 离散单元法理论与实现
MuGdxy (Xinyu Lu) - 性能优化
提交分支 – 按照比赛规则要求，去除对其他库的import引用，可能导致运行效率变低：
开发分支 – 持续迭代和优化：
本例提供一个完整的离散单元法 (DEM) 仿真的实现，考虑了复杂的DEM力学机制，仿真结果达到工程定量精度。
本例使用Taichi语言，加上适当的数据结构和算法，以保障计算效率。
相比Taichi DEM原始版本，本例增加了如下功能：
二维离散单元法→三维离散单元法；
完整考虑并实现颗粒方位和旋转，预留了建模和仿真非球形颗粒的空间；
实现了墙 (离散元中的边界) 单元，并实现了颗粒-墙接触解算；
实现了复杂离散元接触模型，包括一个胶结模型 (爱丁堡胶结颗粒模型，EBPM) 和一个颗粒接触模型 (Hertz-Mindlin接触模型)；
因胶结模型已实现，可以使用胶结团块模拟非球形颗粒；
因胶结模型已实现，可以仿真颗粒破碎过程；
将材料特性与颗粒/墙体关联，以降低空间成本；
将接触面交互属性与接触关联，以降低空间成本；
基于Morton编码实现空间哈希表，进行邻接搜索 (宽碰撞检测)；
分别对低工作负载和高工作负载采用位表和并行扫描算法，存储邻接对以减少内核内部的发散，从而提高并行计算效率；
接触存储在每个颗粒所关联的动态列表中，以降低空间成本，列表在每个时间步都会维护（处理接触的增减）。
原理将另行上传附件说明。
本例展示了开仑台球开局第一球。白球正对其他球运动并发生碰撞，然后球散开。虽然该过程中有能量损失，但是由于没有抗转动力学响应耗散转动动能，所有的球在进入纯滚动状态后将会一直滚动下去。本例可用于验证Hertz-Mindlin模型。
本例展示了一个立方体胶结团块撞击平整表面的过程。
撞击过程中，团块内部的胶结受力将发生破坏，然后团块碎成碎片，飞向周围空间。
本例可用于验证EBPM模型。
本例与上例相似，唯独团块所含颗粒数不同。本例可用于大规模体系仿真的性能测试。
本例展示了一个斯坦福兔子形状的胶结团块自由落体撞击水平表面，然后破碎的过程。本例可用于大规模体系仿真的性能测试。
感谢谢菲尔德大学化学与生物工程学院Xizhong Chen博士对本研究给予的帮助和支持。
9 个赞
附件补充原理说明：HPC_Part.pdf (410.2 KB)DEM_Part.pdf (307.9 KB)
此外，因为发现不同学科对于同一概念使用的术语不同，为避免歧义，我们约定统一使用以下术语，这里进行补充解释。
DEM一个时间步包含5个流程：
邻接搜索 - 接触判断 - 接触解析 - 接触评估 - 时间积分
邻接搜索 Neighboring search == 碰撞检测宽检测 Broad phase collision detection
接触判断 Contact detection == 碰撞检测窄检测 Narrow phase collision detection
接触解析 Contact resolution ∈ 碰撞后的物理计算
接触评估 Contact evaluation ∈ 碰撞后的物理计算
时间积分 Time integration
邻接对 Neighboring pair
邻接对列表 Neighboring pair list
接触列表 Contact list
2 个赞
注：如果作品没有对算法进行拓展和优化，而仅仅是调整粒子配置实现视觉效果展示的也欢迎在评论区提交，我们也会送出 Taichi 周边套装~纯视觉展示不计入最后代码优化评选。
。主要的工作在于三维邻居搜索和复杂几何边界的交互，但是基于三角网格的几何边界实现有亿点点麻烦，索性就和MatDEM一样，拿固定的颗粒描述边界，其位置不根据受力更新。
感谢
