Field建立的问题

F = ti.field(dtype=ti.f32, shape=N_total[None])
这个N_total 每个循环不一样大。
我存在两个问题:
1、循环步第一步我的结果是正确的,第二步之后,这个F就全为0(给他赋不上值);
2、循环步超过512步,报错。

RuntimeError: [D:/a/taichi/taichi/taichi/system/snode_tree_buffer_manager.cpp:allocate@47] LLVM backend supports up to 512 snode trees

Hi @dudu, 目前field还不支持 shape 动态改变。可以考虑用 Manual field allocation and destruction

这个方法用了,会重新启动kernel。
image

请问:这是为什么呢?
MBuilder = ti.linalg.SparseMatrixBuilder(N_total[None], N_total[None], max_num_triplets)是因为这个也不能是动态的吗?

目前是这样的。需要提前定义好系数矩阵。

建立动态稀疏矩阵这个目前没有办法解决吗?

这是我的code。freedom()是计算每一步的N_total[None],有什么办法可以解决这个问题吗?

def solve(): 
    MBuilder = ti.linalg.SparseMatrixBuilder(N_total[None], N_total[None], max_num_triplets)
    KBuilder = ti.linalg.SparseMatrixBuilder(N_total[None], N_total[None], max_num_triplets)
    assemble_M(MBuilder)    #质量一致矩阵
    M = MBuilder.build()
    # print(M)
    compute_Sti_element(KBuilder)   #单元刚度
    K = KBuilder.build()
    # print(K)
    A = M + K
    # print(A)
    solver = ti.linalg.SparseSolver(solver_type="LU")  #求解器
    solver.analyze_pattern(A)                          #分析稀疏矩阵
    solver.factorize(A)                                #分解稀疏矩阵
    fb1 = ti.FieldsBuilder()
    F = ti.field(dtype=ti.f32)
    fb1.dense(ti.i,N_total[None]).place(F)
    fb1_snode_tree = fb1.finalize()
    print(F)
    compute_RHS(F) 
    # print(F)
    boundary()                         
    new_u = solver.solve(F)                         #得到位移结果
    print(new_u)
    print(solver.info())
    update_newton(new_u)
    fb1_snode_tree.destroy()
initialize()
for i in range(5):
    P2G()
    boundary()
    freedom()
    # print(N_total[None],'freedom')
    solve()
    boundary()

请问老师,现在在taichi里能否实现先定义好一个field的大小后,后面再增加或缩小field的大小。
类似于python的list.append()的功能呢?
谢谢老师!

你可以看一下 dynamic snode

感谢您的建议!但是请问老师是不是没有针对struct field的方法,在dynamic snode中,我关注到利用root关键词只能建立field的dynamic snode,谢谢老师!

我没看出来有什么问题,你能分享一下能跑的代码么。我看一下哪里出问题了。

其实struct field对应的功能,按说用advanced field定义方式也可以实现的。目前struct field 应该没有动态改变shape的功能。

这个我的代码,应该就是稀疏矩阵那个shape每一步不一样导致的问题。稀疏矩阵是不是也能像field一样每一步.destroy()?
期待您的解答。

import taichi as ti  #引入taichi模块的同时为该模块取名为ti 
ti.init(arch=ti.cpu,default_fp=ti.f32,default_ip=ti.i32)
quality = 1 
n_particles, n_grid= 200 * quality **2,  50 * quality 
dx = 1  
inv_dx =  1/dx
dt = 1e-3 / quality    
p_vol, p_rho = (dx)**2, 1000  
p_mass = p_vol * p_rho/4  
E, nu = 1.0e6, 0.0 
G=E/(2*(1+nu))    
K_modulus=E/(3*(1-2*nu))  #体积模量
gravity = ti.Vector([0,-9.8])  #每个粒子的重力
x = ti.Vector.field(2, dtype=ti.f32, shape=n_particles)           # position
ini_x = ti.Vector.field(2, dtype=ti.f32, shape=n_particles) 
grid_u = ti.Vector.field(2, dtype=ti.f32, shape=(n_grid, n_grid))
grid_fext = ti.Vector.field(2, dtype=ti.f32, shape=(n_grid, n_grid))
grid_m = ti.field(dtype=ti.f32, shape=(n_grid, n_grid))           #grid node mass
base = ti.Vector.field(2,dtype = ti.f32, shape = n_particles)     #粒子坐标整数值
fx = ti.Vector.field(2,dtype = ti.f32, shape = n_particles)       #粒子坐标小数值(偏移整数坐标的量)
w = ti.Matrix.field(2,2,dtype = ti.f32, shape = n_particles)      #形函数
de_w=ti.Matrix.field(2,2,dtype = ti.f32, shape = n_particles)     #形函数导数                     
m_rho = ti.field(dtype = ti.f32, shape = n_particles)             #粒子密度 
extract = ti.Vector.field(6, dtype=ti.f32,shape=n_particles) 
#-----------------------------------------------------------------------------------
max_num_triplets = 10000000
#-----------------------------------------------------------------------------------
N_total = ti.field(dtype=ti.i32,shape=())
number_v = ti.field(dtype=ti.i32,shape=(n_grid, n_grid))

@ti.kernel
def initialize():   # 初始条件,位置、速度、密度
    for i in range(n_particles):      
        x[i] = [i%2*0.5+1.25,i//2*0.5+0.25]
        ini_x[i] = [i%2*0.5+1.25,i//2*0.5+0.25]
        m_rho[i] = 1000 
        
@ti.kernel
def P2G():
    for i, j in grid_m:         #清空网格上速度、质量信息
        grid_m[i, j] = 0.0        
        grid_u[i, j] = [0.0, 0.0]  
    for p in range(n_particles) :  
        base[p] = (x[p] * inv_dx).cast(ti.i32)  #base为基础下标  x的第p个粒子  每个数×inv_dx,然后把全转成整数
        fx[p] = x[p] * inv_dx - base[p].cast(ti.f32)              #fx表示粒子偏移网格位置 即该粒子距离格子左下角的局部坐标 
        w0=[1-fx[p],fx[p]]  #形函数
        for i,j in ti.static(ti.ndrange(2,2)): #组装形函数矩阵
            w[p][i,j]=w0[i][j]
        de_w = [-1,1]  #形函数导数
        for i,j in ti.static(ti.ndrange(2,2)):                  #遍历相邻四个格子  
            offset=ti.Vector([i,j])      
            weight = w[p][i,0] * w[p][j,1]                      #权重比例
            grid_m[base[p] + offset] += weight * p_mass         #网格点上的质量信息   
          
@ti.kernel
def assemble_M(A: ti.linalg.sparse_matrix_builder()): #组装稀疏质量矩阵
    n_free = 0
    for k in range(1):
        for i,j in ti.ndrange(n_grid, n_grid):
            if grid_m[i,j] > 0:
                n_free += 1
                ele_num = int(N_total[None]/2)
                A[(n_free-1), (n_free-1)] += 4*grid_m[i,j]/(dt**2)
                A[(n_free-1) + ele_num, (n_free-1) + ele_num] += 4*grid_m[i,j]/(dt**2)
            else:
                n_free += 0
                
@ti.kernel
def compute_RHS(v: ti.template()):
    n_free = 0
    for i,j in grid_m:
        if grid_m[i,j] > 0:
            grid_fext[i,j] += 1.0
    for k in range(1):
        for i,j in ti.ndrange(n_grid, n_grid):
            if grid_m[i,j] > 0:
                n_free += 1
                ele_num = int(N_total[None]/2)
                v[(n_free-1)] = 1.0
                v[(n_free-1) + ele_num] = 2.0               
            else:
                n_free += 0

@ti.kernel
def update_newton(new_u: ti.ext_arr()):
    n_free = 0
    for k in range(1):
        for i,j in ti.ndrange(n_grid, n_grid):
            if grid_m[i,j] > 0:
                n_free += 1
                ele_num = int(N_total[None]/2)
                grid_u[i,j][0] = new_u[n_free-1]
                grid_u[i,j][1] = new_u[(n_free-1) + ele_num]
            else:
                n_free += 0
                grid_u[i,j][0] = 0.0
                grid_u[i,j][1] = 0.0      
def solve(): 
    MBuilder = ti.linalg.SparseMatrixBuilder(N_total[None], N_total[None], max_num_triplets)
    KBuilder = ti.linalg.SparseMatrixBuilder(N_total[None], N_total[None], max_num_triplets)
    assemble_M(MBuilder)    #质量一致矩阵
    M = MBuilder.build()
    A = M 
    # print(A)
    solver = ti.linalg.SparseSolver(solver_type="LU")  #求解器
    solver.analyze_pattern(A)                          #分析稀疏矩阵
    solver.factorize(A)                                #分解稀疏矩阵
    fb2 = ti.FieldsBuilder()
    F = ti.field(dtype=ti.f32)
    fb2.dense(ti.i, N_total[None]).place(F)
    fb2_snode_tree = fb2.finalize()
    compute_RHS(F)          
    new_u = solver.solve(F)                         #得到位移结果
    # print(solver.info())
    update_newton(new_u)
    fb2_snode_tree.destroy()
initialize()
for i in range(500):
    P2G()
    N_total[None] = 2*i
    # print(N_total[None],'freedom')
    solve()