From 1ca5df2751cc72c1731f8aabad18af88b828dff6 Mon Sep 17 00:00:00 2001 From: peterctang <150650021+peterctang@users.noreply.github.com> Date: Mon, 6 Jan 2025 18:36:41 -0500 Subject: [PATCH 1/4] Add files via upload --- test2d.py | 84 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) create mode 100644 test2d.py diff --git a/test2d.py b/test2d.py new file mode 100644 index 0000000..1e14e78 --- /dev/null +++ b/test2d.py @@ -0,0 +1,84 @@ +import os +import numpy as np +import taichi as ti + +from projects.mpm.engine.mpm_solver_implicit import MPMSolverImplicit +from common.utils.logger import * + +ti.init(arch=ti.cpu, default_fp=ti.f64) +#ti.init(arch=ti.gpu, default_fp=ti.f64) + +test_case = 10 +fps = 24 +max_num_frame = 200 + +# Remove or comment out the output directory creation +directory = 'output/2d/' + str(test_case) + '/' +os.makedirs(directory, exist_ok=True) + +gui = ti.GUI("PNMPM-2D", res=512, background_color=0x112F41) +#""" +mpm = MPMSolverImplicit(res=(128, 128)) + +if test_case == 0: # jello drop on ground + mpm.symplectic = True + mpm.setDXandDT(DX=0.0078125, DT=0.00015625) + mpm.setGravity((0, -1.0)) + mpm.setLameParameter(E=40, nu=0.2) + + mpm.add_cube(min_corner=(0.5, 0.6), + max_corner=(0.6, 0.9), + num_particles=3333, + rho=2) + mpm.add_cube(min_corner=(0.4, 0.2), + max_corner=(0.7, 0.5), + num_particles=10000, + rho=2) + mpm.add_analytic_box(min_corner=(0.0, 0.0), max_corner=(1.0, 0.05)) + +if test_case == 10: # jello drop on mound + mpm.symplectic = False + mpm.cfl = 0.6 + #mpm.setDXandDT(DX=0.00502513, DT=0.004) + mpm.setDXandDT(DX=0.00502513, DT=0.004) + mpm.setGravity((0, -2.0)) + mpm.setLameParameter(E=40, nu=0.2) + + mpm.add_cube(min_corner=(0.3, 0.7), + max_corner=(0.7, 0.9), + num_particles=20000, + rho=2) + mpm.add_analytic_box(min_corner=(0.4, 0.3), + max_corner=(0.6, 0.5), + rotation=(3.1415926 / 4, 0.0, 0.0, 0.0)) + mpm.add_analytic_box(min_corner=(0.0, 0.0), max_corner=(1.0, 0.05)) + +if test_case == 1: # physbam snow + mpm.symplectic = False + mpm.cfl = 0.6 + mpm.setDXandDT(DX=0.00502513, DT=0.004) + # mpm.symplectic = True + # mpm.setDXandDT(DX=0.00502513, DT=0.0005) + mpm.setGravity((0, -2.0)) + mpm.setLameParameter(E=40, nu=0.2) + + mpm.load_state(filename='data/2d.txt') + mpm.add_analytic_box(min_corner=(0.4, 0.3), + max_corner=(0.6, 0.5), + rotation=(3.1415926 / 4, 0.0, 0.0, 0.0)) + mpm.add_analytic_box(min_corner=(0.0, 0.0), max_corner=(1.0, 0.05)) +#""" + +with Logger(directory + f'log.txt'): + for frame in range(max_num_frame): + print("================== frame", frame, "===================") + #mpm.step(8e-3) + #print(mpm.p_v[10])#test initial vel + mpm.step(1 / fps) + particles = mpm.particle_info() + gui.circles(particles['position'], radius=1.5, color=0x068587) + gui.show(directory + f'{frame:06d}.png') + #print(mpm.getMaxVelocity()) + #gui.show() + +ti.print_profile_info() From 5c6842fc2caea4446a66c2cc31e740aeca2659f1 Mon Sep 17 00:00:00 2001 From: peterctang <150650021+peterctang@users.noreply.github.com> Date: Mon, 6 Jan 2025 18:42:44 -0500 Subject: [PATCH 2/4] Delete test2d.py --- test2d.py | 84 ------------------------------------------------------- 1 file changed, 84 deletions(-) delete mode 100644 test2d.py diff --git a/test2d.py b/test2d.py deleted file mode 100644 index 1e14e78..0000000 --- a/test2d.py +++ /dev/null @@ -1,84 +0,0 @@ -import os -import numpy as np -import taichi as ti - -from projects.mpm.engine.mpm_solver_implicit import MPMSolverImplicit -from common.utils.logger import * - -ti.init(arch=ti.cpu, default_fp=ti.f64) -#ti.init(arch=ti.gpu, default_fp=ti.f64) - -test_case = 10 -fps = 24 -max_num_frame = 200 - -# Remove or comment out the output directory creation -directory = 'output/2d/' + str(test_case) + '/' -os.makedirs(directory, exist_ok=True) - -gui = ti.GUI("PNMPM-2D", res=512, background_color=0x112F41) -#""" -mpm = MPMSolverImplicit(res=(128, 128)) - -if test_case == 0: # jello drop on ground - mpm.symplectic = True - mpm.setDXandDT(DX=0.0078125, DT=0.00015625) - mpm.setGravity((0, -1.0)) - mpm.setLameParameter(E=40, nu=0.2) - - mpm.add_cube(min_corner=(0.5, 0.6), - max_corner=(0.6, 0.9), - num_particles=3333, - rho=2) - mpm.add_cube(min_corner=(0.4, 0.2), - max_corner=(0.7, 0.5), - num_particles=10000, - rho=2) - mpm.add_analytic_box(min_corner=(0.0, 0.0), max_corner=(1.0, 0.05)) - -if test_case == 10: # jello drop on mound - mpm.symplectic = False - mpm.cfl = 0.6 - #mpm.setDXandDT(DX=0.00502513, DT=0.004) - mpm.setDXandDT(DX=0.00502513, DT=0.004) - mpm.setGravity((0, -2.0)) - mpm.setLameParameter(E=40, nu=0.2) - - mpm.add_cube(min_corner=(0.3, 0.7), - max_corner=(0.7, 0.9), - num_particles=20000, - rho=2) - mpm.add_analytic_box(min_corner=(0.4, 0.3), - max_corner=(0.6, 0.5), - rotation=(3.1415926 / 4, 0.0, 0.0, 0.0)) - mpm.add_analytic_box(min_corner=(0.0, 0.0), max_corner=(1.0, 0.05)) - -if test_case == 1: # physbam snow - mpm.symplectic = False - mpm.cfl = 0.6 - mpm.setDXandDT(DX=0.00502513, DT=0.004) - # mpm.symplectic = True - # mpm.setDXandDT(DX=0.00502513, DT=0.0005) - mpm.setGravity((0, -2.0)) - mpm.setLameParameter(E=40, nu=0.2) - - mpm.load_state(filename='data/2d.txt') - mpm.add_analytic_box(min_corner=(0.4, 0.3), - max_corner=(0.6, 0.5), - rotation=(3.1415926 / 4, 0.0, 0.0, 0.0)) - mpm.add_analytic_box(min_corner=(0.0, 0.0), max_corner=(1.0, 0.05)) -#""" - -with Logger(directory + f'log.txt'): - for frame in range(max_num_frame): - print("================== frame", frame, "===================") - #mpm.step(8e-3) - #print(mpm.p_v[10])#test initial vel - mpm.step(1 / fps) - particles = mpm.particle_info() - gui.circles(particles['position'], radius=1.5, color=0x068587) - gui.show(directory + f'{frame:06d}.png') - #print(mpm.getMaxVelocity()) - #gui.show() - -ti.print_profile_info() From 21cd606c3bc675be7c9ea3486c199db225a58ba5 Mon Sep 17 00:00:00 2001 From: peterctang <150650021+peterctang@users.noreply.github.com> Date: Mon, 6 Jan 2025 18:45:01 -0500 Subject: [PATCH 3/4] Add files via upload --- projects/mpm/engine/mpm_solver_implicit.py | 113 ++++++++++++++++++--- 1 file changed, 97 insertions(+), 16 deletions(-) diff --git a/projects/mpm/engine/mpm_solver_implicit.py b/projects/mpm/engine/mpm_solver_implicit.py index dc4a69f..1dd5845 100644 --- a/projects/mpm/engine/mpm_solver_implicit.py +++ b/projects/mpm/engine/mpm_solver_implicit.py @@ -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) @@ -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) @@ -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, @@ -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)): @@ -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): @@ -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 @@ -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 @@ -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) @@ -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() @@ -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] @@ -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 From 3eef3b895c7f8958d8ab1a8e394a297c80ec54dc Mon Sep 17 00:00:00 2001 From: peterctang <150650021+peterctang@users.noreply.github.com> Date: Mon, 6 Jan 2025 18:54:05 -0500 Subject: [PATCH 4/4] Add files via upload --- projects/mpm/basic/fixed_corotated.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/projects/mpm/basic/fixed_corotated.py b/projects/mpm/basic/fixed_corotated.py index 13fecf6..8115b0d 100644 --- a/projects/mpm/basic/fixed_corotated.py +++ b/projects/mpm/basic/fixed_corotated.py @@ -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() @@ -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) @@ -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 @@ -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]