在taichi AOT demo “implicit FEM”，将element K matrix组装成stiffness matrix过程中，并没有存成矩阵的形式，而是存成Vert+Edge的形式
z = 0
for u_i in ti.static(range(4)):
u = verts[u_i]
for v_i in ti.static(range(4)):
v = verts[v_i]
if u < v:
hes_edge[c2e[c][z]] += hes[u_i * 3, v_i * 3]
z += 1
for zz in ti.static(range(4)):
hes_vert[verts[zz]] += hes[zz * 3, zz * 3]
line155-164 “implicit_fem.py”
import argparse
import os
import pathlib
import shutil
import numpy as np
import taichi as ti
parser = argparse.ArgumentParser()
parser.add_argument('--dim', type=int, default=3)
parser.add_argument('--aot', default=False, action='store_true')
args = parser.parse_args()
# TODO: asserts cuda or vulkan backend
ti.init(arch=ti.vulkan)
SCRIPT_PATH = os.path.dirname(os.path.realpath(__file__))
def get_rel_path(*segs):
show original
请问在这种情况下，如果想要做matrix preconditioner，如何能够在这种非矩阵的储存形式上完成呢？还是需要转化为矩阵再进行预处理？
对此有些疑惑，非常感谢
1 Like
YuPeng
May 10, 2022, 10:02am
#2
Hi @Jingzheng , 非常欢迎来到太极论坛。
针对你的问题，我觉得有两种方案：
如果 preconditioner 比较简单，可以考虑在目前这个结构上手写。vertex 上的矩阵对应主对角线，edge 上的矩阵对应其余位置。
存成一般的 sparse matrix 的形式，之后用通用的方法进行处理，比如用 scipy 或是 gpu 上的 sparse matrix。
非常感谢！关于这个demo我还有个小问题，之前在调试代码的时候，get_matrix()函数中有一个初始化Hessian矩阵步骤
hes = ti.Matrix.zero(ti.f32, 12, 12)
然后给出的建议是替换成ti.field(ti.f32,12,12)
可是我在替换成ti.field后，使用ti.ndrange初始化，在并行的情况下会得到错误的结果（观察到的是Hessian没有完全清0），在开cpu单线程的情况下就是完全正确的。请问在这种情况下应该怎么来初始化这个ti.field呢？
YuPeng
May 11, 2022, 4:21am
#4
你可以暂时先不用管这个warning。
初始化ti.field 是有 addition 操作才造成结果错误的?atomic_add 或者 += .
我这里也是每一步初始化hes为0矩阵，就直接用ndrange定义为0
hes = ti.field(ti.f32, shape=(12,12))
@ti.kernel
def get_matrix(c2e: ti.types.ndarray(), vertices: ti.types.ndarray()):
for c in vertices:
verts = vertices[c]
W_c = W[c]
B_c = B[c]
#hes = ti.Matrix.zero(ti.f32, 12, 12)
for i,j in ti.ndrange(12,12):
hes[i,j] = 0
for u in ti.static(range(4)):
for d in ti.static(range(3)):
dD = ti.Matrix.zero(ti.f32, 3, 3)
if ti.static(u == 3):
for j in ti.static(range(3)):
dD[d, j] = -1
else:
dD[d, u] = 1
dF = dD @ B_c
dP = 2.0 * mu * dF
dH = -W_c * dP @ B_c.transpose()
for i in ti.static(range(3)):
for j in ti.static(range(3)):
hes[i * 3 + j, u * 3 + d] = -dt**2 * dH[j, i]
hes[3 * 3 + j, u * 3 + d] += dt**2 * dH[j, i]
请问在多线程情况下怎么处理会更好呢？（不太理解这里为什么会出现错误
YuPeng
May 13, 2022, 1:49pm
#6
在Taichi的 ti.kernel 函数里面，两个 并列的 for 语句之间是线性执行的。
@ti.kernel
def test_for():
for i in range(10):
....
# 前面的for循环执行完之后，再执行下面的语句
for j in range(10):
....
YuPeng
May 13, 2022, 1:52pm
#7
你这里是不是会出现不同线程的race condition?
，在ti.init中取cpu_max_num_threads = 1结果就是正常的，所以是因为多线程在hes没有完全初始化为0就参与到后面的计算中了吗？）