.397395619494580:e589aa2f1246dcb64182b291ddc4a559_69e6210c180bcf8a442e5ffe.69e62c42180bcf8a442e60b6.69e62c42966af563bdeec98f:Trae CN.T(2026/4/20 21:38:10)

demo/demo_3d.py

@@ -1,12 +1,16 @@
+import sys
+import os
+sys.path.insert(0, os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), 'local_packages'))
+
 import taichi as ti
 import numpy as np
 import utils
 from engine.mpm_solver import MPMSolver
 
 write_to_disk = False
 
-# Try to run on GPU
-ti.init(arch=ti.cuda, device_memory_GB=4.0)
+# Try to run on CPU (for testing)
+ti.init(arch=ti.cpu)
 
 gui = ti.GUI("Taichi Elements", res=512, background_color=0x112F41)
 
@@ -25,6 +29,12 @@
 
 mpm.set_gravity((0, -50, 0))
 
+mpm.add_wind_field(
+    lower_bound=[0, 0, 0],
+    upper_bound=[3, 10, 10],
+    force=[200, 0, 0]
+)
+
 for frame in range(1500):
     mpm.step(4e-3)
     colors = np.array([0x068587, 0xED553B, 0xEEEEF0, 0xFFFF00],

demo/test_wind_field.py

@@ -0,0 +1,60 @@
+import sys
+import os
+import warnings
+
+warnings.filterwarnings("ignore")
+os.environ["TI_LOG_LEVEL"] = "error"
+
+project_root = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+sys.path.insert(0, project_root)
+sys.path.insert(0, os.path.join(project_root, 'local_packages'))
+
+import taichi as ti
+import numpy as np
+from engine.mpm_solver import MPMSolver
+
+ti.init(arch=ti.cpu, log_level=ti.ERROR)
+
+print("=== Testing Wind Field Feature ===")
+print()
+
+mpm = MPMSolver(res=(32, 32, 32), size=10)
+
+mpm.add_cube(lower_corner=[1, 4, 4],
+             cube_size=[1, 1, 1],
+             material=MPMSolver.material_water,
+             velocity=[0, 0, 0])
+
+mpm.set_gravity((0, -20, 0))
+
+mpm.add_wind_field(
+    lower_bound=[0, 0, 0],
+    upper_bound=[3, 10, 10],
+    force=[100, 0, 0]
+)
+
+print("Initial particles:")
+particles = mpm.particle_info()
+print(f"  Number of particles: {len(particles['position'])}")
+print(f"  Initial X positions range: [{particles['position'][:, 0].min():.3f}, {particles['position'][:, 0].max():.3f}]")
+print(f"  Initial mean X position: {particles['position'][:, 0].mean():.3f}")
+print()
+
+print("Running simulation for 100 steps...")
+for frame in range(100):
+    mpm.step(4e-3)
+    if frame % 20 == 0:
+        particles = mpm.particle_info()
+        mean_x = particles['position'][:, 0].mean()
+        mean_vx = particles['velocity'][:, 0].mean()
+        print(f"  Step {frame:3d}: mean X = {mean_x:.3f}, mean Vx = {mean_vx:.3f}")
+
+print()
+print("Final state:")
+particles = mpm.particle_info()
+print(f"  Final X positions range: [{particles['position'][:, 0].min():.3f}, {particles['position'][:, 0].max():.3f}]")
+print(f"  Final mean X position: {particles['position'][:, 0].mean():.3f}")
+print(f"  Final mean X velocity: {particles['velocity'][:, 0].mean():.3f}")
+print()
+print("=== Wind Field Test Passed! ===")
+print("Particles should have moved to the right due to the wind field.")

demo/wind_field_visual_test.py

@@ -0,0 +1,97 @@
+import sys
+import os
+import warnings
+
+warnings.filterwarnings("ignore")
+os.environ["TI_LOG_LEVEL"] = "error"
+
+project_root = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+sys.path.insert(0, project_root)
+sys.path.insert(0, os.path.join(project_root, 'local_packages'))
+
+import taichi as ti
+import numpy as np
+from engine.mpm_solver import MPMSolver
+
+ti.init(arch=ti.cpu)
+
+print("=" * 60)
+print("Wind Field Visual Test")
+print("=" * 60)
+print()
+print("Scene setup:")
+print("  - Simulation size: 10 units, 64x64x64 resolution")
+print("  - Wind field: X from 0 to 3, Y from 0 to 10, Z from 0 to 10")
+print("  - Wind force: (150, 0, 0) - pushing particles to the right")
+print("  - Particles: water cube at (1, 4, 4) with size 1x1x1")
+print("  - Gravity: (0, -30, 0)")
+print()
+
+mpm = MPMSolver(res=(64, 64, 64), size=10)
+
+mpm.add_cube(lower_corner=[1, 4, 4],
+             cube_size=[1, 1, 1],
+             material=MPMSolver.material_water,
+             velocity=[0, 0, 0])
+
+mpm.set_gravity((0, -30, 0))
+
+mpm.add_wind_field(
+    lower_bound=[0, 0, 0],
+    upper_bound=[3, 10, 10],
+    force=[150, 0, 0]
+)
+
+gui = ti.GUI("Wind Field Test", res=800, background_color=0x112F41)
+
+colors = np.array([0x068587, 0xED553B, 0xEEEEF0, 0xFFFF00], dtype=np.uint32)
+
+print("Press ESC to exit. Press SPACE to pause/resume.")
+print()
+
+paused = False
+frame = 0
+
+while gui.running:
+    for e in gui.get_events(ti.GUI.PRESS):
+        if e.key == ti.GUI.ESCAPE:
+            gui.running = False
+        elif e.key == ti.GUI.SPACE:
+            paused = not paused
+            print(f"Paused: {paused}")
+
+    if not paused:
+        mpm.step(4e-3)
+        frame += 1
+
+    particles = mpm.particle_info()
+    np_x = particles['position'] / 10.0
+
+    screen_x = ((np_x[:, 0] + np_x[:, 2]) / 2**0.5) - 0.2
+    screen_y = np_x[:, 1]
+    screen_pos = np.stack([screen_x, screen_y], axis=-1)
+
+    gui.circles(screen_pos, radius=2.0, color=colors[particles['material']])
+
+    gui.line((0.0, 0.0), (0.3, 0.0), color=0xFF0000, radius=3)
+    gui.line((0.3, 0.0), (0.3, 1.0), color=0xFF0000, radius=3)
+
+    mean_x = particles['position'][:, 0].mean()
+    mean_vx = particles['velocity'][:, 0].mean()
+
+    gui.text(f"Frame: {frame}", pos=(0.02, 0.98), color=0xFFFFFF, font_size=20)
+    gui.text(f"Mean X: {mean_x:.2f}", pos=(0.02, 0.93), color=0xFFFFFF, font_size=20)
+    gui.text(f"Mean Vx: {mean_vx:.2f}", pos=(0.02, 0.88), color=0xFFFFFF, font_size=20)
+
+    wind_zone_text = "Wind Zone (X < 3)"
+    gui.text(wind_zone_text, pos=(0.02, 0.83), color=0xFF6666, font_size=16)
+
+    if frame % 50 == 0 and not paused:
+        print(f"Frame {frame:4d}: Mean X = {mean_x:6.3f}, Mean Vx = {mean_vx:6.3f}")
+
+    gui.show()
+
+print()
+print("=" * 60)
+print("Test completed!")
+print("=" * 60)

engine/mpm_solver.py

@@ -318,6 +318,26 @@ def set_gravity(self, g):
         assert len(g) == self.dim
         self.gravity[None] = g
 
+    def add_wind_field(self, lower_bound, upper_bound, force):
+        lower_bound = list(lower_bound)
+        upper_bound = list(upper_bound)
+        force = list(force)
+        assert len(lower_bound) == self.dim
+        assert len(upper_bound) == self.dim
+        assert len(force) == self.dim
+
+        @ti.kernel
+        def apply_wind(t: ti.f32, dt: ti.f32, grid_v: ti.template()):
+            for I in ti.grouped(grid_v):
+                grid_pos = I * self.dx
+                inside = True
+                for d in ti.static(range(self.dim)):
+                    inside = inside and (lower_bound[d] <= grid_pos[d] < upper_bound[d])
+                if inside:
+                    grid_v[I] += dt * ti.Vector(force)
+
+        self.grid_postprocess.append(apply_wind)
+
     @ti.func
     def sand_projection(self, sigma, p):
         sigma_out = ti.Matrix.zero(ti.f32, self.dim, self.dim)

local_packages/bin/get_gprof

@@ -0,0 +1,75 @@
+#!E:\ProgramFiles\Anaconda\python.exe
+#
+# Author: Mike McKerns (mmckerns @caltech and @uqfoundation)
+# Copyright (c) 2008-2016 California Institute of Technology.
+# Copyright (c) 2016-2026 The Uncertainty Quantification Foundation.
+# License: 3-clause BSD.  The full license text is available at:
+#  - https://github.com/uqfoundation/dill/blob/master/LICENSE
+'''
+build profile graph for the given instance
+
+running:
+  $ get_gprof <args> <instance>
+
+executes:
+  gprof2dot -f pstats <args> <type>.prof | dot -Tpng -o <type>.call.png
+
+where:
+  <args> are arguments for gprof2dot, such as "-n 5 -e 5"
+  <instance> is code to create the instance to profile
+  <type> is the class of the instance (i.e. type(instance))
+
+For example:
+  $ get_gprof -n 5 -e 1 "import numpy; numpy.array([1,2])"
+
+will create 'ndarray.call.png' with the profile graph for numpy.array([1,2]),
+where '-n 5' eliminates nodes below 5% threshold, similarly '-e 1' eliminates
+edges below 1% threshold
+'''
+
+if __name__ == "__main__":
+    import sys
+    if len(sys.argv) < 2:
+        print ("Please provide an object instance (e.g. 'import math; math.pi')")
+        sys.exit()
+    # grab args for gprof2dot
+    args = sys.argv[1:-1]
+    args = ' '.join(args)
+    # last arg builds the object
+    obj = sys.argv[-1]
+    obj = obj.split(';')
+    # multi-line prep for generating an instance
+    for line in obj[:-1]:
+        exec(line)
+    # one-line generation of an instance
+    try:
+        obj = eval(obj[-1])
+    except Exception:
+        print ("Error processing object instance")
+        sys.exit()
+
+    # get object 'name'
+    objtype = type(obj)
+    name = getattr(objtype, '__name__', getattr(objtype, '__class__', objtype))
+
+    # profile dumping an object
+    import dill
+    import os
+    import cProfile
+    #name = os.path.splitext(os.path.basename(__file__))[0]
+    cProfile.run("dill.dumps(obj)", filename="%s.prof" % name)
+    msg = "gprof2dot -f pstats %s %s.prof | dot -Tpng -o %s.call.png" % (args, name, name)
+    try:
+        res = os.system(msg)
+    except Exception:
+        print ("Please verify install of 'gprof2dot' to view profile graphs")
+    if res:
+        print ("Please verify install of 'gprof2dot' to view profile graphs")
+
+    # get stats
+    f_prof = "%s.prof" % name
+    import pstats
+    stats = pstats.Stats(f_prof, stream=sys.stdout)
+    stats.strip_dirs().sort_stats('cumtime')
+    stats.print_stats(20) #XXX: save to file instead of print top 20?
+    os.remove(f_prof)

local_packages/bin/get_objgraph

@@ -0,0 +1,54 @@
+#!E:\ProgramFiles\Anaconda\python.exe
+#
+# Author: Mike McKerns (mmckerns @caltech and @uqfoundation)
+# Copyright (c) 2008-2016 California Institute of Technology.
+# Copyright (c) 2016-2026 The Uncertainty Quantification Foundation.
+# License: 3-clause BSD.  The full license text is available at:
+#  - https://github.com/uqfoundation/dill/blob/master/LICENSE
+"""
+display the reference paths for objects in ``dill.types`` or a .pkl file
+
+Notes:
+    the generated image is useful in showing the pointer references in
+    objects that are or can be pickled.  Any object in ``dill.objects``
+    listed in ``dill.load_types(picklable=True, unpicklable=True)`` works.
+
+Examples::
+
+    $ get_objgraph ArrayType
+    Image generated as ArrayType.png
+"""
+
+import dill as pickle
+#pickle.debug.trace(True)
+#import pickle
+
+# get all objects for testing
+from dill import load_types
+load_types(pickleable=True,unpickleable=True)
+from dill import objects
+
+if __name__ == "__main__":
+    import sys
+    if len(sys.argv) != 2:
+        print ("Please provide exactly one file or type name (e.g. 'IntType')")
+        msg = "\n"
+        for objtype in list(objects.keys())[:40]:
+            msg += objtype + ', '
+        print (msg + "...")
+    else:
+        objtype = str(sys.argv[-1])
+        try:
+            obj = objects[objtype]
+        except KeyError:
+            obj = pickle.load(open(objtype,'rb'))
+            import os
+            objtype = os.path.splitext(objtype)[0]
+        try:
+            import objgraph
+            objgraph.show_refs(obj, filename=objtype+'.png')
+        except ImportError:
+            print ("Please install 'objgraph' to view object graphs")
+
+
+# EOF

local_packages/bin/undill

@@ -0,0 +1,22 @@
+#!E:\ProgramFiles\Anaconda\python.exe
+#
+# Author: Mike McKerns (mmckerns @caltech and @uqfoundation)
+# Copyright (c) 2008-2016 California Institute of Technology.
+# Copyright (c) 2016-2026 The Uncertainty Quantification Foundation.
+# License: 3-clause BSD.  The full license text is available at:
+#  - https://github.com/uqfoundation/dill/blob/master/LICENSE
+"""
+unpickle the contents of a pickled object file
+
+Examples::
+
+    $ undill hello.pkl
+    ['hello', 'world']
+"""
+
+if __name__ == '__main__':
+    import sys
+    import dill
+    for file in sys.argv[1:]:
+        print (dill.load(open(file,'rb')))
+

local_packages/colorama-0.4.6.dist-info/INSTALLER

@@ -0,0 +1 @@
+pip