[Taichi] version 1.6.0, llvm 15.0.4, commit f1c6fbbd, linux, python 3.11.3
[Taichi] Starting on arch=x64
[E 09/30/23 11:20:50.480 226305] [type.h:TypedConstant@605] Not supported.
[E 09/30/23 11:20:50.480 226305] [ir.cpp:~DelayedIRModifier@408] Assertion failure: to_insert_before_.empty()
terminate called after throwing an instance of ‘std::__cxx11::basic_string<char, std::char_traits, std::allocator >’
已放弃 (核心已转储)
(base) user@user-Syste
简而言之我在ti.func()里面用if（if）【两个嵌套，然后报错】，我也很迷茫。这个过程是在复现CD-MPM中出现的。求教，因为代码略长方便可以放出代码。
部分节选：
import taichi as ti
import numpy as np
import math
ti.init(arch=ti.cpu)
quality = 1
d = 2
n_particles, n_grid = 20000 * quality ** 2, 128 * quality
E_s, nu_s = 3.537e5, 0.3 # 沙子的杨氏模量和泊松比
mu_s, lambda_s = E_s / (2 * (1 + nu_s)), E_s * nu_s / ((1 + nu_s) * (1 - 2 * nu_s)) # 沙的拉梅常数
xi = 1
kappa = (2/3)*mu_s+lambda_s
logJp = ti.field(dtype = float, shape = n_particles) # 屈服表面尺寸
beta = 0.5
M = 2.36
state = ti.field(dtype = int, shape = n_particles) #状态
F_s = ti.Matrix.field(2, 2, dtype = float, shape = n_particles) # deformation gradient
########################################################################################
#NACC Plasticity
########################################################################################
@ti.func
def deviatoric(AMatrix):
AMatrix_hat = AMatrix - AMatrix.trace() / d * ti.Matrix.identity(float, 2)
return AMatrix_hat
@ti.func
def projectStrainNACC(F_s,p):
U, sig, V = ti.svd(F_s)
for_sinh = xi * ti.max(-logJp[p],0)
p0 = kappa * ((0.00001) + 0.5 * (ti.exp(for_sinh) - ti.exp(-for_sinh))) #0.5 * (ti.exp(J) - ti.exp(-J))
B_hat_trial = ti.Matrix.identity(float, d)
J = sig[1,1]*sig[0,0]
B_hat_trial[0,0] = sig[0,0] * sig[0,0]
B_hat_trial[1,1] = sig[1,1] * sig[1,1]
s_hat_trial = mu_s * ti.pow(J,-2/d) * deviatoric(B_hat_trial)
prime = kappa / 2 * (J - 1 / J)
p_trial = -prime * J
y_s_half_coeff = (6 - d) / 2 * (1 + 2 * beta)
y_p_half = M * M * (p_trial + beta * p0) * (p_trial - p0)
y = y_s_half_coeff * (s_hat_trial**2).sum() + y_p_half
p_min = beta * p0
F_new = U @ sig @ V.transpose()
#Case 1:Project to max tip of YS
if(p_trial > p0):
Je_new = ti.sqrt(-2 * p0 / kappa + 1)
sig = ti.Matrix([[1.0, 0.0],[0.0,1.0]]) * ti.pow(Je_new, 1 / d)
logJp[p] += ti.log(J / Je_new)
F_new = U @ sig @ V.transpose()
state[p] = 1
#Csat 2:Project to min tip of YS
elif (p_trial < -p_min):
Je_new = ti.sqrt(2 * p_min / kappa + 1)
sig = ti.Matrix([[1.0, 0.0],[0.0,1.0]]) * ti.pow(Je_new, 1 / d)
logJp[p] += ti.log(J / Je_new)
F_new = U @ sig @ V.transpose()
state[p] = 0
#弹性阶段
else : #(p0 > 1e-4 and p_trial < p0 - 1e-4 and p_trial > 1e-4 - p_min):
p_center = (p0 - p_min) * .5
q_trial = ti.sqrt((6 - d) / 2) * s_hat_trial.norm()
direction = ti.Vector([0.0,0.0])
direction[0] = p_center - p_trial
direction[1] = 0 - q_trial
direction = direction / direction.norm()
C = M * M * (p_center + beta * p0) * (p_center - p0)
B = M * M * direction[0] * (2 * p_center - p0 + beta * p0)
A = M * M * direction[0] * direction[0] + (1 + 2 * beta) * direction[1] * direction[1]
l1 = (-B + ti.sqrt(B * B - 4 * A * C)) / (2 * A)
l2 = (-B - ti.sqrt(B * B - 4 * A * C)) / (2 * A)
p1 = p_center + l1 * direction[0]
p2 = p_center + l2 * direction[0]
print(q_trial)
p_fake = p2
state[p] = 2
if (p_trial - p_center) * (p1 - p_center) > 0 :
p_fake = p1
else :
p_fake = p2
Je_new_fake = ti.sqrt(ti.abs(-2 * p_fake / kappa + 1))
#if(Je_new_fake > 1e-4):
#logJp[p] += ti.log(J / Je_new_fake)
################################################bug所在前后区域
B_hat_new = ti.pow(J, 2 / d) / mu_s * ti.sqrt(-y_p_half / y_s_half_coeff) * s_hat_trial / s_hat_trial.norm()
#B_hat_new += 1 / d * B_hat_trial.sum() * ti.Matrix([[1.0, 0.0],[0.0,1.0]])
#print(B_hat_new)
#for i in range(d):
# sig[i,i] = ti.sqrt(B_hat_new[i])
return F_new
@ti.kernel
def come():
for p in range(2):
F_s[p] = ti.Matrix([[1.0, -0.98880],[-0.7880,1.0]])
new_F = projectStrainNACC(F_s[p],p)
come()