由于本人还是苦逼本科生+疫情把final推迟到了暑假末，所以没有在ddl前赶上hw2（其实那个时候也才刚啃完MPM）。最近趁复习摸鱼时间重新把hw2做了，复现了蒋老师的两篇和沙子有关的论文：
Multi-species simulation of porous sand and water mixtures 和 Drucker-Prager Elastoplasticity for Sand Animation，用MPM实现了沙子和水的耦合。
干沙子 (10000个沙子粒子，128x128网格)
湿沙子+水 （10000个沙子粒子，每次迭代+50个水粒子，128x128网格，约26fps）
湿沙子+水 显示出project操作的不同状态（和论文中使用相同的着色）
干沙子+水，根据我的理解，论文在耦合时默认沙子初始是湿的并没有考虑毛细作用，如果我代码没有写错的话，感觉这个结果并不是非常正确。
复现的过程可谓是非常痛苦，踩了许多的坑，据不完全统计有以下这些：
- hardening操作的参数 h0,h1,h2,h3 文本中提到的是弧度制而参考表中给的是角度制，虽然下面有一行说明但是一开始没看到。。
- 一开始分别实现了沙子和水，耦合的时候忽略修正沙子和水的密度，导致数量级差距过大动量项交换异常。
- 应力求错了，不得不说在simulator里调试NAN相关的问题非常痛苦。
- 如果沙子处在屈服面的范围内，project返回的时候不考虑论文里修正体积的项
- 在处理 cohesion 的时候论文里把不等式右端项直接换成了 cC 但是没有给出明确的推导。。本人数理基础较差，最后是参考 这里 的方法，类似于体积修正项在 project 操作前给 e 额外加上 \frac{c_C}{d*\alpha_p} \mathbf{I} 的项，也就是：
e += (c_C0[p] * (1.0 - phi_s[p])) / (d * alpha_s[p]) * ti.Matrix.identity(float, 2) # effects of cohesion
此外还有几点对17年论文的疑惑，希望能够得到解答：
- section 3.1.1的应力公式貌似写错了，应该是 \frac{\lambda}{2} \mathbf{tr^2(\epsilon)} ，原文是 \frac{\lambda}{2} \mathbf{tr(\epsilon)}
- 在实现中，感觉体积修正项的正负号和论文里的正好相反，按论文里写貌似不对？
- section 4.3.4 末尾关于更新 \mathbf{F}^{s E, n+1} 的公式貌似是错的？
代码托管在：
https://github.com/g1n0st/GAMES201/blob/master/hw2/sand_water_mpm.py
鸣谢：
非常感谢 archibate 学长热心帮忙解决基础问题：NAN check and throw exception on debug mode · Issue #1777 · taichi-dev/taichi · GitHub
蒋老师在 GAMESWebinar 上关于本论文的报告：GAMESWebinar2017-09æ-èéå¡å¤«_è
¾è®¯è§é¢
总结：
Taichi 复现论文真的是很好用的工具，希望以后有时间多写写hw2-extra！
大家战斗力太强了…
刚刚看到了仓库