搞流体模拟的隐式粘度求解器,直接求解线性系统是真算不动啊。网上没搜到很直观的关于matrix-free方法的介绍啊555(具体如 matrix-free conjugate gradient。有好一些的参考资料吗?)
Matrix-free 的求解器可以参考一下 Scipy 的文档:scipy.sparse.linalg.LinearOperator — SciPy v1.10.1 Manual
目前 Taichi 里也有了一个 matrix-free 的 CG 求解器:[lang] Implement experimental CG(Conjugate Gradient) solver in Taichi-lang by houkensjtu · Pull Request #7690 · taichi-dev/taichi · GitHub
使用方法暂时可以参考测试脚本:taichi/test_taichi_cg.py at master · taichi-dev/taichi · GitHub
简单说就是你可以定义一个 LinearOperator
来代替矩阵 A ,然后传入 CG 求解器来求解方程组。有问题的话欢迎继续讨论!
1 Like
非常感谢,我先学习学习~
这个测试脚本已经失效,请问还有别的参考嘛
谢谢指出!我来整理一下最新的版本稍后发在 repo 里~
最近的测试脚本在这里:https://github.com/taichi-dev/taichi/blob/master/tests/python/test_matrixfree_cg.py
下面是一个直接可运行的版本:
import taichi as ti
from taichi.linalg import taichi_cg_solver, LinearOperator # Use taichi_cg_solver in 1.6.0
# from taichi.linalg import MatrixFreeCG, LinearOperator # Use MatrixFreeCG from 1.7.0 (not released yet)
ti.init(arch=ti.cpu)
n = 4
x = ti.field(dtype=ti.f32, shape=(n, n))
Ax = ti.field(dtype=ti.f32, shape=(n, n))
b = ti.field(dtype=ti.f32, shape=(n, n))
@ti.kernel
def init():
for i in ti.grouped(b):
b[i] = 1.0
x[i] = 0.0
@ti.kernel
def compute(v:ti.template(), mv:ti.template()):
for i in ti.grouped(v):
mv[i] = 2 * v[i]
init()
A = LinearOperator(compute)
print(b)
taichi_cg_solver(A, b, x, maxiter=10*n*n, tol=1e-18, quiet=True)
# MatrixFreeCG(A, b, x, maxiter=10*n*n, tol=1e-18, quiet=True) # Use MatrixFreeCG if you're on master branch
print(x)
注意如果你是 <=1.6.0 的版本,那么需要使用 taichi_cg_solver
这个调用名。如果是源码编译的,调用名是 MatrixFreeCG
。另外,1.6.0 版本还有一个 bug 导致只能求解 2 维问题,新的 master 上已经修复了不过还没反映在 pip 版本里。