请教Implicit FEM demo中的关于stiffness matrix的问题

在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”

请问在这种情况下,如果想要做matrix preconditioner,如何能够在这种非矩阵的储存形式上完成呢?还是需要转化为矩阵再进行预处理?

对此有些疑惑,非常感谢

Hi @Jingzheng , 非常欢迎来到太极论坛。

针对你的问题,我觉得有两种方案:

  1. 如果 preconditioner 比较简单,可以考虑在目前这个结构上手写。vertex 上的矩阵对应主对角线,edge 上的矩阵对应其余位置。
  2. 存成一般的 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呢?

你可以暂时先不用管这个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]

请问在多线程情况下怎么处理会更好呢?(不太理解这里为什么会出现错误 :thinking:,在ti.init中取cpu_max_num_threads = 1结果就是正常的,所以是因为多线程在hes没有完全初始化为0就参与到后面的计算中了吗?)

在Taichi的 ti.kernel 函数里面,两个 并列的 for 语句之间是线性执行的。

@ti.kernel
def test_for():
    for i in range(10):  
         ....
    # 前面的for循环执行完之后,再执行下面的语句
    for j in range(10):
         ....

你这里是不是会出现不同线程的race condition?