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

活动介绍

近期不少同学提到对于 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 吗?

3 Likes

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

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

DEM (1)

视频:

1 Like

我也来抛个砖,做了个搅拌机的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

4 Likes

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

我这就是个玩具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 Like

2022.9.29更新:

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

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

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

《偷功》

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

效果:

v02

代码:GitHub - mrzhuzhe/taichi_dem: A minimal DEM simulation demo written in Taichi.

已实现:

  • 事实上我们可以将邻域搜索的范围从 3 x 3 进一步缩小一半左右,你能想想怎么构建这个搜索算法吗?
  • 其实文中提到的并行计算方式完全可以用 ti.atomic_add 去掉所有串行区域,进一步并行化,你觉得应该咋搞?

PBF 开发中

2 Likes

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

代码在这: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弹,我也想不出哪里出错了,是因为一定要在前一步的切向力上跟新吗?

1 Like

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

1 Like