22# Taichi implementation by Ye Kuang (k-ye)
33
44import math
5-
5+ import os
66import numpy as np
77
88import taichi as ti
@@ -26,6 +26,7 @@ def round_up(f, s):
2626grid_size = (round_up (boundary [0 ], 1 ), round_up (boundary [1 ], 1 ))
2727
2828dim = 2
29+ max_frames = 1800
2930bg_color = 0x112F41
3031particle_color = 0x068587
3132boundary_color = 0xEBACA2
@@ -64,6 +65,7 @@ def round_up(f, s):
6465position_deltas = ti .Vector .field (dim , float )
6566# 0: x-pos, 1: timestep in sin()
6667board_states = ti .Vector .field (2 , float )
68+ wall_force = ti .Vector .field (dim , float , shape = ())
6769
6870ti .root .dense (ti .i , num_particles ).place (old_positions , positions , velocities )
6971grid_snode = ti .root .dense (ti .ij , grid_size )
@@ -122,8 +124,10 @@ def confine_position_to_boundary(p):
122124 bmin = particle_radius_in_world
123125 bmax = ti .Vector ([board_states [None ][0 ], boundary [1 ]]) - particle_radius_in_world
124126 for i in ti .static (range (dim )):
125- # Use randomness to prevent particles from sticking into each other after clamping
126127 if p [i ] <= bmin :
128+ diff = bmin - p [i ]
129+ if i == 0 :
130+ ti .atomic_add (wall_force [None ][i ], 1000.0 * diff / time_delta )
127131 p [i ] = bmin + epsilon * ti .random ()
128132 elif bmax [i ] <= p [i ]:
129133 p [i ] = bmax [i ] - epsilon * ti .random ()
@@ -145,6 +149,7 @@ def move_board():
145149
146150@ti .kernel
147151def prologue ():
152+ wall_force [None ] = ti .Vector ([0.0 , 0.0 ]) # ← Clear previous step force
148153 # save old positions
149154 for i in positions :
150155 old_positions [i ] = positions [i ]
@@ -294,13 +299,36 @@ def main():
294299 init_particles ()
295300 print (f"boundary={ boundary } grid={ grid_size } cell_size={ cell_size } " )
296301 gui = ti .GUI ("PBF2D" , screen_res )
297- while gui .running and not gui .get_event (gui .ESCAPE ):
302+
303+ # Prepare force tracking
304+ force_values = []
305+ time_values = []
306+
307+ # Output file
308+ force_filename = "forces.evt"
309+ if os .path .exists (force_filename ):
310+ os .remove (force_filename ) # clean if it exists
311+
312+ while gui .running and gui .frame < max_frames and not gui .get_event (gui .ESCAPE ):
298313 move_board ()
299314 run_pbf ()
315+
316+ # Record time and force
317+ current_time = gui .frame * time_delta
318+ time_values .append (current_time )
319+ fx = wall_force .to_numpy ()[None ][0 ].item (0 ) # x-component of the wall force
320+ force_values .append (fx )
321+
300322 if gui .frame % 20 == 1 :
301323 print_stats ()
324+ print ("Left wall force:" , fx )
325+
302326 render (gui )
303327
328+ # Write to file at end
329+ with open (force_filename , "w" ) as f :
330+ f .write (" " .join ("0.0" for _ in force_values ) + "\n " )
331+ f .write (" " .join (f"{ fval :.5f} " for fval in force_values ) + "\n " )
304332
305333if __name__ == "__main__" :
306334 main ()
0 commit comments