【活动已结束】Taichi DEM 代码优化挑战赛来了!Airpods 胶囊咖啡机还有更多社区周边等你来拿~

投票已关闭

截止 10/28 共收到 9 份 DEM 代码优化的投稿。最佳人气奖的投票正式开启!请大家为自己喜欢的作品投上一票吧~

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
代码评审&投票:10/31-11/4

参赛方式

作品提交

  1. 在 GitHub 上使用 "use this template"从给定的 Template repo 中复制出自己的新 repo 并公开 ,在此基础上形成作品 repo,需要包含源代码 dem.py 文件、README.md 文件和 requiremements.txt。

  2. 在 dem.py 中添加你的代码,加速或者优化 DEM,实现知乎文末提到的各种算法拓展。

  3. 提交作品:请在本帖评论区提交作品。提交内容为:

  4. 简单介绍自己的作品优化了 DEM 哪个部分,如果有原理说明更佳~

  5. 附上作品的 GitHub Repo 地址

  6. 附上作品的 GIF 动图

可参考 @llinus 这位同学的投稿:

:star: 注:如果作品没有对算法进行拓展和优化,而仅仅是调整粒子配置实现视觉效果展示的也欢迎在评论区提交,我们也会送出 Taichi 周边套装~纯视觉展示不计入最后代码优化评选。

参赛规则

  • 本次的代码优化挑战赛规则很简单,没有代码行数限制。仅有以下两个要求:
      1. 仅可改动 dem.py 完成代码优化。
      1. 不可以 import 其他库

奖项

本次大赛设置一二三等奖和最佳人气奖。所有奖项都将在 10/31-11/4 期间开启评审&投票。最佳人气奖投票通道将在 10/31 日开启~

:star: 注:所有参与投稿的小伙伴都可以获得 Taichi 周边套装哦~

脑洞时刻

还没有思路的同学,可以看看以下几个优化方向哦~

  • 事实上我们可以将邻域搜索的范围从 3 x 3 进一步缩小一半左右,你能想想怎么构建这个搜索算法吗?

  • 现在的算例中为了简化计算细节,并没有考虑 DEM 中常用到的粒子的角动量和不规则形状,你能扩展本算例并加入这些新的 feature 吗?

  • 其实 Taichi 对于 GPU 上的动态数据结构有很好的支持,如 taichi_elements 中就采用了 dynamic SNode 来维护每个网格 cell 中的粒子下标数组。你能尝试用 Taichi 的dynamic SNode 来实现本文的功能吗?

  • 其实 Position-based fluid (PBF) 也可以用这种方法进行加速,你可以试一试吗?

  • 其实文中提到的并行计算方式完全可以用 ti.atomic_add 去掉所有串行区域,进一步并行化,你觉得应该咋搞?

  • 本文的代码比较简单,完全是 2D 的,你能把它扩展成 3D 吗?

5 个赞

我先来抛砖引玉,祝大家中秋快乐!

代码仓库:GitHub - yuanming-hu/taichi_dem_mid_autumn
代码效果:

DEM (1)

视频:

3 个赞

我也来抛个砖,做了个搅拌机的demo,这是内置于商业离散元软件PFC中的案例,GPU的计算效率相比之下是真的高 :grin:。主要的工作在于三维邻居搜索和复杂几何边界的交互,但是基于三角网格的几何边界实现有亿点点麻烦,索性就和MatDEM一样,拿固定的颗粒描述边界,其位置不根据受力更新。

将来会慢慢实现更符合真实物理的接触本构模型和三角网格描述的边界,以求得和PFC相当的计算精度。

代码仓库:GitHub - Linus-Civil/GeoBlender: A 3D demo of blender in civil/geotechnical/geological engineering based on Taichi DEM
代码效果:

ezgif.com-gif-maker

5 个赞

大佬请问有详细的算法描述吗,主要是颗粒和中间的叶片作用是怎么实现的?

我这就是个玩具Demo,代码里只实现了法向力。你如果对工程级别的力学本构细节感兴趣,推荐你去看看《计算颗粒力学及工程应用(季顺迎)》。

:pray:感谢

请问为啥case2 输出的结果里粒子多了很多不存在的碰撞?
case1 却是正常的?
正在查c++ 原子操作的说明

case 1 ok

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]

case 2 wrong

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]   
  • 其实文中提到的并行计算方式完全可以用 ti.atomic_add 去掉所有串行区域,进一步并行化,你觉得应该咋搞?

对应这个问题 提了一个PR给你

1 个赞

2022.9.29更新:

  1. 原子操作,由 @mrzhuzhe 贡献
  2. 折半邻居网格数枚举
1 个赞

请问在你的代码中折半邻居搜索后计算效率会有所提升嘛,我试过折半搜索,同时去掉了相邻网格时id1要小于id2的条件,照理来说减少搜索次数速度应该提升,但发现效率反而下降了

会提升的,i<j的条件需要选择性保留

《偷功》

实现了一个简单旋转小球的效果

效果:

v02

PBF 版本: 增加了不可压缩性,粘度等, 这个实现参考了 Ten Minute Physics 中的实现
ezgif.com-gif-maker

代码: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 去掉所有串行区域,进一步并行化,你觉得应该咋搞?

后续问题:

  1. 在GPU3090上 黄色 dem 旋转小球,如果临域搜索的数据结构是用dynamic snode 只有 30-40 fps
    但是用blog里的紧凑内存实现 60-70 fps, 瓶颈在哪里?
  2. 3d的情况下粒子靠近边界时会卡在一起,会产生数值下溢,目前只能随机把小球往回移动一小段距离,这个方法有点过于trick,有没有通用的解决办法?
 if (pos[i][0] >= maxX):
            pos[i][0] = maxX - eps * ti.random();
  1. fp32 在小步长下很容易数值溢出,是否有通用的方法解决?
3 个赞

岩土专业,代码小白。在实现过程中遇到一点问题,想要请教社区的大家

代码在这:GitHub - chenqing253/taichiDEM

尝试实现加入tangential force and rotation, 遇到了以下几点问题:

  1. 如果根据以下公式来进行 切向力的更新
    image
    where F and M are force and moment fro entity i and j,and F^s is given by

    and the increment tangential force is \Delta F^s = -k^s * \Delta U^s , 其中 斜向分量 \Delta U^s is given by relative velocity V^s
    image

t_i is Tangential unit vector. 上述就是求两个颗粒在接触切平面上的相对滑移来求摩擦力

可以看出切向力可能需要在前一步的接触上的切向力进行更新,有什么方法或者数据结构可以记录 前一步的接触吗

  1. 我忽略了需要在前一步的切向力的基础上更新这一条件,尝试写了一下,结果不正确(我觉得比较合理的结果是形成一个沙丘状)。如图~

ezgif.com-gif-maker (1)

我用红色和蓝色标注了不同旋转方向的颗粒,旋转了但是低于平均值的为白色。而且颗粒模拟过程中比较Q弹,我也想不出哪里出错了,是因为一定要在前一步的切向力上跟新吗?

2 个赞

赵仕威老师有篇文章里有写A thread-block-wise computational framework for large-scale hierarchical continuum-discrete modeling of granular media,采用多个列表存储数据

1 个赞

我也来分享一个小的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 个赞

使用 Taichi DEM 进行工程定量离散单元法仿真

从工程角度使用太极语言完整实现离散单元法。

carom

cube_911_particles_impact

cube_18112_particles_impact

使用BIMBase ™二次开发可视化。BIMBase ™是北京构力科技有限公司研发的用于BIM的图形平台。BIMBase开发平台-BIMBase Platform

bunny

使用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原始版本,本例增加了如下功能:

  1. 二维离散单元法→三维离散单元法;

  2. 完整考虑并实现颗粒方位和旋转,预留了建模和仿真非球形颗粒的空间;

  3. 实现了墙 (离散元中的边界) 单元,并实现了颗粒-墙接触解算;

  4. 实现了复杂离散元接触模型,包括一个胶结模型 (爱丁堡胶结颗粒模型,EBPM) 和一个颗粒接触模型 (Hertz-Mindlin接触模型);

  5. 因胶结模型已实现,可以使用胶结团块模拟非球形颗粒;

  6. 因胶结模型已实现,可以仿真颗粒破碎过程;

  7. 将材料特性与颗粒/墙体关联,以降低空间成本;

  8. 将接触面交互属性与接触关联,以降低空间成本;

  9. 基于Morton编码实现空间哈希表,进行邻接搜索 (宽碰撞检测);

  10. 分别对低工作负载和高工作负载采用位表和并行扫描算法,存储邻接对以减少内核内部的发散,从而提高并行计算效率;

  11. 接触存储在每个颗粒所关联的动态列表中,以降低空间成本,列表在每个时间步都会维护(处理接触的增减)。

原理将另行上传附件说明。

示例

开仑台球

本例展示了开仑台球开局第一球。白球正对其他球运动并发生碰撞,然后球散开。虽然该过程中有能量损失,但是由于没有抗转动力学响应耗散转动动能,所有的球在进入纯滚动状态后将会一直滚动下去。本例可用于验证Hertz-Mindlin模型。

carom

911个颗粒胶结组成的立方体团块撞击平整表面

本例展示了一个立方体胶结团块撞击平整表面的过程。

撞击过程中,团块内部的胶结受力将发生破坏,然后团块碎成碎片,飞向周围空间。

本例可用于验证EBPM模型。

cube_911_particles_impact

18112个颗粒胶结组成的立方体团块撞击平整表面

本例与上例相似,唯独团块所含颗粒数不同。本例可用于大规模体系仿真的性能测试。

cube_18112_particles_impact

斯坦福兔子自由落体撞击水平表面

本例展示了一个斯坦福兔子形状的胶结团块自由落体撞击水平表面,然后破碎的过程。本例可用于大规模体系仿真的性能测试。

bunny

致谢

感谢谢菲尔德大学化学与生物工程学院Xizhong Chen博士对本研究给予的帮助和支持。

10 个赞

附件补充原理说明:
HPC_Part.pdf (410.2 KB)
DEM_Part.pdf (307.9 KB)

此外,因为发现不同学科对于同一概念使用的术语不同,为避免歧义,我们约定统一使用以下术语,这里进行补充解释。

DEM一个时间步包含5个流程:

邻接搜索 - 接触判断 - 接触解析 - 接触评估 - 时间积分
Neighboring search - Contact detection - Contact resolution - Contact evaluation - Time integration

  1. 邻接搜索 Neighboring search == 碰撞检测宽检测 Broad phase collision detection
    此过程解决两个颗粒 (==粒子,英文均为particle) 是否相邻,因为只有相邻的颗粒才有可能接触。虽然该过程对于DEM算法本身并非必要,但是此过程可以剪枝很大一部分冗余计算,大大提升计算效率。

  2. 接触判断 Contact detection == 碰撞检测窄检测 Narrow phase collision detection
    此过程解决两个颗粒是否真正有接触 (== 碰撞, contact == collision)。接触的意思是,二者之间存在一个相互作用 (interaction),可在颗粒间传导力 (矩)。

  3. 接触解析 Contact resolution ∈ 碰撞后的物理计算
    此过程解决两个颗粒的接触的几何特性,以作为接触评估时输入接触模型力-位移表达式的位移数据。

  4. 接触评估 Contact evaluation ∈ 碰撞后的物理计算
    此过程调用预设的接触模型 (contact model),接触模型中含有力-位移表达式来表达接触的力学响应。通过输入位移数据,计算后得到接触传导的力 (矩),然后作用于接触两端的颗粒上。

  5. 时间积分 Time integration
    此过程汇总颗粒所受的所有力 (矩),利用牛顿第二定律和时间积分得到颗粒的运动学数据,然后更新体系。

  6. 邻接对 Neighboring pair
    为两个经邻接搜索确认相邻的颗粒的编号 (ID) 所组成的数对 {ID_i, ID_j}。

  7. 邻接对列表 Neighboring pair list
    存储所有邻接对的列表。

  8. 接触列表 Contact list
    存储所有接触的列表。接触也是以 {ID_i, ID_j} 的模式进行存储,与邻接对类似。

2 个赞