Skip to content

addressing some taichi issues #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
peterctang wants to merge 4 commits into
base: master
Choose a base branch
from
Open

addressing some taichi issues #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1ca5df2
Add files via upload
peterctang Jan 6, 2025
5c6842f
Delete test2d.py
peterctang Jan 6, 2025
21cd606
Add files via upload
peterctang Jan 6, 2025
3eef3b8
Add files via upload
peterctang Jan 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 9 additions & 8 deletions projects/mpm/basic/fixed_corotated.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,18 +1,19 @@
import taichi as ti
from projects.mpm.basic.math_tools import *
from taichi.lang import impl


@ti.func
def elasticity_energy(sig: ti.template(), la, mu):
if ti.static(sig.n == 2):
sigma = ti.Matrix.zero(ti.get_runtime().default_fp, sig.n, 1)
sigma = ti.Vector.zero(impl.get_runtime().default_fp, sig.n)
for i in ti.static(range(sig.n)):
sigma[i] = sig[i, 0 if ti.static(sig.m == 1) else i]
sigmam12Sum = (sigma - ti.Vector([1, 1])).norm_sqr()
sigmaProdm1 = sigma[0] * sigma[1] - 1
return mu * sigmam12Sum + la / 2 * sigmaProdm1 * sigmaProdm1
else:
sigma = ti.Matrix.zero(ti.get_runtime().default_fp, sig.n, 1)
sigma = ti.Matrix.zero(impl.get_runtime().default_fp, sig.n, 1)
for i in ti.static(range(sig.n)):
sigma[i] = sig[i, 0 if ti.static(sig.m == 1) else i]
sigmam12Sum = (sigma - ti.Vector([1, 1, 1])).norm_sqr()
Expand All @@ -23,16 +24,16 @@ def elasticity_energy(sig: ti.template(), la, mu):
@ti.func
def elasticity_gradient(sig: ti.template(), la, mu):
if ti.static(sig.n == 2):
sigma = ti.Matrix.zero(ti.get_runtime().default_fp, sig.n, 1)
sigma = ti.Vector.zero(impl.get_runtime().default_fp, sig.n)
for i in ti.static(range(sig.n)):
sigma[i] = sig[i, 0 if ti.static(sig.m == 1) else i]
sigma[i] = sig[i]
sigmaProdm1lambda = la * (sigma[0] * sigma[1] - 1)
sigmaProd_noI = ti.Vector([sigma[1], sigma[0]])
_2u = mu * 2
return ti.Vector([_2u * (sigma[0] - 1) + sigmaProd_noI[0] * sigmaProdm1lambda,
_2u * (sigma[1] - 1) + sigmaProd_noI[1] * sigmaProdm1lambda])
else:
sigma = ti.Matrix.zero(ti.get_runtime().default_fp, sig.n, 1)
sigma = ti.Matrix.zero(impl.get_runtime().default_fp, sig.n, 1)
for i in ti.static(range(sig.n)):
sigma[i] = sig[i, 0 if ti.static(sig.m == 1) else i]
sigmaProdm1lambda = la * (sigma[0] * sigma[1] * sigma[2] - 1)
Expand All @@ -46,9 +47,9 @@ def elasticity_gradient(sig: ti.template(), la, mu):
@ti.func
def elasticity_hessian(sig: ti.template(), la, mu):
if ti.static(sig.n == 2):
sigma = ti.Matrix.zero(ti.get_runtime().default_fp, sig.n, 1)
sigma = ti.Vector.zero(impl.get_runtime().default_fp, sig.n)
for i in ti.static(range(sig.n)):
sigma[i] = sig[i, 0 if ti.static(sig.m == 1) else i]
sigma[i] = sig[i]
sigmaProd = sigma[0] * sigma[1]
sigmaProd_noI = ti.Vector([sigma[1], sigma[0]])
_2u = mu * 2
Expand All @@ -57,7 +58,7 @@ def elasticity_hessian(sig: ti.template(), la, mu):
[la * ((sigmaProd - 1) + sigmaProd_noI[0] * sigmaProd_noI[1]),
_2u + la * sigmaProd_noI[1] * sigmaProd_noI[1]]])
else:
sigma = ti.Matrix.zero(ti.get_runtime().default_fp, sig.n, 1)
sigma = ti.Matrix.zero(impl.get_runtime().default_fp, sig.n, 1)
for i in ti.static(range(sig.n)):
sigma[i] = sig[i, 0 if ti.static(sig.m == 1) else i]
sigmaProd = sigma[0] * sigma[1] * sigma[2]
Expand Down
113 changes: 97 additions & 16 deletions projects/mpm/engine/mpm_solver_implicit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,6 @@ def up_to_2_pow_n(x: int):

#### Declare taichi fields for simulation data
self.gravity = ti.Vector.field(self.dim, dtype=real, shape=())
self.pid = ti.field(ti.i32)

# position
self.p_x = ti.Vector.field(self.dim, dtype=real)
Expand Down Expand Up @@ -109,8 +108,22 @@ def up_to_2_pow_n(x: int):

offset = tuple(-self.grid_size // 2 for _ in range(self.dim))
self.offset = offset
"""
self.pid = ti.field(ti.i32)
# grid node momentum/velocity
self.grid_v = ti.Vector.field(self.dim, dtype=real)
self.grid_f = ti.Vector.field(self.dim, dtype=real)
self.grid_idx = ti.field(dtype=ti.i32)
self.grid_dv = ti.Vector.field(self.dim, dtype=real)
# grid node mass
self.grid_m = ti.field(dtype=real)

self.num_active_grid = ti.field(dtype=ti.i32, shape=())

grid_block_size = 128
self.grid = ti.root.pointer(indices, self.grid_size // grid_block_size) # 32
"""
self.pid = ti.field(ti.i32)
# grid node momentum/velocity
self.grid_v = ti.Vector.field(self.dim, dtype=real)
self.grid_f = ti.Vector.field(self.dim, dtype=real)
Expand All @@ -136,19 +149,28 @@ def block_component(c):
block.dense(indices, self.leaf_block_size).place(c, offset=offset) # 16 (-2048, 2048)

block_component(self.grid_m)
for v in self.grid_v.entries:
block_component(v)
for dv in self.grid_dv.entries:
block_component(dv)
for f in self.grid_f.entries:
block_component(f)
#for v in self.grid_v.entries:
#block_component(v)
#for dv in self.grid_dv.entries:
#block_component(dv)
#for f in self.grid_f.entries:
#block_component(f)
for v in range(self.dim):
block_component(self.grid_v.get_scalar_field(v))
for dv in range(self.dim):
block_component(self.grid_dv.get_scalar_field(dv))
for f in range(self.dim):
block_component(self.grid_f.get_scalar_field(f))

block_component(self.grid_idx)

block_offset = tuple(o // self.leaf_block_size
for o in self.offset)

block.dynamic(ti.indices(self.dim),
block.dynamic(ti.axes(self.dim),
1024 * 1024,
chunk_size=self.leaf_block_size**self.dim * 8).place(
self.pid, offset=offset + (0, ))
self.pid, offset=block_offset + (0, ))

self.particle = ti.root.dynamic(ti.i, max_num_particles, 2**17) # 2**20 will trigger error in cuda + 3d
self.particle.place(self.p_x, self.p_v, self.p_C, self.p_F, self.p_Fp,
Expand Down Expand Up @@ -205,25 +227,60 @@ def block_component(c):

self.delta = ti.field(real, shape=()) # Store a double precision delta

@ti.func
def stencil_range(self):
return ti.ndrange(*((3, ) * self.dim))

"""
@ti.kernel
def build_pid(self):
ti.block_dim(64)
ti.loop_config(block_dim=64)
for p in self.p_x:
base = int(ti.floor(self.p_x[p] * self.inv_dx - 0.5))
ti.append(self.pid.parent(), base - ti.Vector(list(self.offset)),
ti.append(self.pid.parent(), base - ti.Vector(self.offset),
p)
#print(self.pid.parent())#for bugging
"""
@ti.kernel
def build_pid(self):
"""
grid has blocking (e.g. 4x4x4), we wish to put the particles from each block into a GPU block,
then used shared memory (ti.block_local) to accelerate
:param pid:
:param grid_m:
:param offset:
:return:
"""
ti.loop_config(block_dim=64)
for p in self.p_x:
base = int(ti.floor(self.p_x[p] * self.inv_dx - 0.5)) \
- ti.Vector(self.offset)
# Pid grandparent is `block`
base_pid = ti.rescale_index(self.grid_m, self.pid.parent(2), base)
#print("base is", base)
#print("base_pid is", base_pid)
#print("p is", p)
ti.append(self.pid.parent(), base_pid, p)

#for bugging pid
@ti.kernel
def print_pid(self):
print("Grid Cell")
for I in ti.grouped(self.pid): # Iterate over all grid indices in pid
print(f"Grid Cell: {I}")
print(self.pid[I])
#for i in range(len(self.pid[I])): # Access the particles in the cell
#print(f" Particle: {self.pid[I][i]}")


@ti.kernel
def p2g(self, dt: real):
ti.no_activate(self.particle)
ti.block_dim(256)
ti.loop_config(block_dim=256)#block_dim(256)
# ti.block_local(*self.grid_v.entries)
# ti.block_local(*self.grid_dv.entries) ##########
# ti.block_local(self.grid_m)
for I in ti.grouped(self.pid):
#print("I is", I)########
p = self.pid[I]
base = ti.floor(self.p_x[p] * self.inv_dx - 0.5).cast(int)
for D in ti.static(range(self.dim)):
Expand Down Expand Up @@ -251,6 +308,7 @@ def p2g(self, dt: real):
oidx = self.offset2idx(offset)
self.p_cached_w[p, oidx] = weight
self.p_cached_dw[p, oidx] = dN
#print(self.grid_v[base + offset])#for test

@ti.func
def applyPlasticity(self, dt):
Expand Down Expand Up @@ -282,7 +340,7 @@ def grid_normalization(self):
self.grid_v[I] = (1 / self.grid_m[I]
) * self.grid_v[I] # Momentum to velocity
self.grid_dv[I] = self.grid_v[I]
idx = self.num_active_grid[None].atomic_add(1) # Avoid error in parallel computing
idx = ti.atomic_add(self.num_active_grid[None],1) # Avoid error in parallel computing
# print(idx)
self.grid_idx[I] = idx
self.dof2idx[idx] = I
Expand Down Expand Up @@ -569,6 +627,9 @@ def TotalEnergy(self, dt: real) -> real:
F = self.p_F[p]
la, mu = self.p_la[p], self.p_mu[p]
U, sig, V = svd(F)
#############################
########### 需要再改##########
#######################################
psi = elasticity_energy(sig, la, mu)
ee += self.p_vol[p] * psi

Expand Down Expand Up @@ -982,7 +1043,7 @@ def implicit_newton2(self, dt):

@ti.kernel
def g2p(self, dt: real):
ti.block_dim(256)
ti.loop_config(block_dim=256)#ti.block_dim(256)
# ti.block_local(*self.grid_v.entries)
# ti.block_local(*self.grid_dv.entries)
ti.no_activate(self.particle)
Expand Down Expand Up @@ -1076,6 +1137,7 @@ def evaluatePerNodeCNTolerance(self, eps:real, dt:real):
def advanceOneStepExplicit(self, dt):
self.grid.deactivate_all()
self.build_pid()
#self.print_pid() #for bugging pid
self.p2g(dt)
self.grid_normalization()

Expand Down Expand Up @@ -1257,8 +1319,9 @@ def sample_cube2(self, lower_corner: ti.template(), cube_size: ti.template(), ne
self.p_mu[idx] = self.mu_0
self.n_particles[None] = self.n_particles[None] + new_particles

"""
@ti.kernel
def copy_dynamic_nd(self, np_x: ti.ext_arr(), input_x: ti.template()):
def copy_dynamic_nd(self, np_x: ti.template(), input_x: ti.template()):
for i in self.p_x:
for j in ti.static(range(self.dim)):
np_x[i, j] = input_x[i][j]
Expand All @@ -1273,6 +1336,24 @@ def particle_info(self):
'position': np_x,
'velocity': np_v,
}
"""

@ti.kernel
def copy_dynamic_nd(self, np_x: ti.types.ndarray(), input_x: ti.template()):
for i in self.p_x:
for j in ti.static(range(self.dim)):
np_x[i, j] = input_x[i][j]

def particle_info(self):
np_x = ti.ndarray(dtype=ti.f64, shape=(self.n_particles[None], self.dim))
self.copy_dynamic_nd(np_x, self.p_x)
np_v = ti.ndarray(dtype=ti.f64, shape=(self.n_particles[None], self.dim))
self.copy_dynamic_nd(np_v, self.p_v)

return {
'position': np_x.to_numpy(),
'velocity': np_v.to_numpy(),
}

def simulation_info(self):
i = 0
Expand Down