Add in additional example numerical simulations for TaichiEvent (FEM,…

… MPM, Euler, etc.) and fix runtime pathing

modules/createEVENT/TaichiEvent/TaichiEvent.py

@@ -186,9 +186,8 @@ def writeEVENT(forces, eventFilePath='EVENT.json', floorsCount=1):  # noqa: N802
 
     # subtype = "StochasticWindModel-KwonKareem2006"
     eventClassification = 'Hydro'  # noqa: N806
-    eventType = 'StochasticWave'  # noqa: N806
-    eventSubtype = 'StochasticWaveJonswap'  # noqa: N806, F841
-    # subtype = "StochasticWaveJonswap" # ?
+    eventType = 'TaichiEvent'  # noqa: N806
+    eventSubtype = 'TaichiEvent-OSU-LWF'  # noqa: N806, F841
     # timeSeriesName = "HydroForceSeries_1X"
     # patternName = "HydroForcePattern_1X"
 
@@ -223,39 +222,29 @@ def GetFloorsCount(BIMFilePath):  # noqa: N802, N803, D103
     return int(bim['GeneralInformation']['stories'])
 
 
-def main():  # noqa: D103
-    return 0
-    # """
-    # Entry point to generate event file using Stochastic Waves
-    # """
-    # #CLI parser
-    # parser = argparse.ArgumentParser(description="Get sample EVENT file produced by StochasticWave")
-    # parser.add_argument('-b', '--filenameAIM', help="BIM File", required=True)
-    # parser.add_argument('-e', '--filenameEVENT', help= "Event File", required=True)
-    # parser.add_argument('--getRV', help= "getRV", required=False, action='store_true')
-
-    # #parsing arguments
-    # arguments, unknowns = parser.parse_known_args()
-
-    # exec(open("Ex1_WaveKinematics.py").read())
-    # exec(open("Ex2_Jonswap_Spectrum.py").read())
-    # exec(open("Ex3_WaveTimeSeries.py").read())
-    # # exec(open("Ex4_WaveLoads.py").read())
-
-    # # Run Ex4_WaveLoads.py with the given parameters
-    # # result = Ex4_WaveLoads.main(arguments.water_depth, arguments.peak_period, arguments.significant_wave_height, arguments.pile_diameter, arguments.drag_coefficient, arguments.mass_coefficient, arguments.number_of_recorders_z, arguments.time)
-    # import subprocess
-    # result = subprocess.run(["python", "Ex4_WaveLoads.py", "-hw", 30.0, "-Tp", 12.7, "-Hs", 5.0, "-Dp", 1.0, "-Cd", 2.1, "-Cm", 2.1, "-nz", GetFloorsCount(arguments.filenameAIM), "-t", 10.0], stdout=subprocess.PIPE)
+def GetTaichiScript(BIMFilePath):  # noqa: N802, N803, D103
+    filePath = BIMFilePath  # noqa: N806
+    with open(filePath, encoding='utf-8') as file:  # noqa: PTH123
+        evt = json.load(file)
+    file.close  # noqa: B018
 
-    # if arguments.getRV == True:
-    #     #Read the number of floors
-    #     floorsCount = GetFloorsCount(arguments.filenameAIM)
-    #     forces = []
-    #     for i in range(floorsCount):
-    #         forces.append(FloorForces())
+    fileNameKey = 'interfaceSurface' 
+    filePathKey = fileNameKey + 'Path'
+
+    for event in evt['Events']:
+        fileName = event[fileNameKey]
+        filePath = event[filePathKey]
+        return os.path.join(filePath, fileName)
+
+    defaultScriptPath = f'{os.path.realpath(os.path.dirname(__file__))}'  # noqa: ISC003, PTH120
+    defaultScriptName = 'taichi_script.py'
+    return defaultScriptPath + defaultScriptName 
 
-    #     #write the event file
-    #     writeEVENT(forces, arguments.filenameEVENT)
+def main():  # noqa: D103
+    """
+    Entry point to generate event file using TaichiEvent.
+    """
+    return 0
 
 
 if __name__ == '__main__':
@@ -286,22 +275,24 @@ def main():  # noqa: D103
     # parsing arguments
     arguments, unknowns = parser.parse_known_args()
 
-    # Run Ex4_WaveLoads.py with the given parameters
-    # result = Ex4_WaveLoads.main(arguments.water_depth, arguments.peak_period, arguments.significant_wave_height, arguments.pile_diameter, arguments.drag_coefficient, arguments.mass_coefficient, arguments.number_of_recorders_z, arguments.time)
-
     # import subprocess
 
+    # Get json of filenameAIM
+    scriptName = GetTaichiScript(arguments.filenameAIM)  # noqa: N816        
+
     if arguments.getRV == True:  # noqa: E712
-        print('RVs requested in StochasticWave.py')  # noqa: T201
+        print('RVs requested')  # noqa: T201
         # Read the number of floors
         floorsCount = GetFloorsCount(arguments.filenameAIM)  # noqa: N816
         filenameEVENT = arguments.filenameEVENT  # noqa: N816
 
+
         result = subprocess.run(  # noqa: S603
             [  # noqa: S607
                 'ti',
-                f'{os.path.realpath(os.path.dirname(__file__))}'  # noqa: ISC003, PTH120
-                + '/taichi_script.py',
+                scriptName,
+                # f'{os.path.realpath(os.path.dirname(__file__))}'  # noqa: ISC003, PTH120
+                # + '/taichi_script.py',
             ],
             stdout=subprocess.PIPE,
             check=False,
@@ -315,13 +306,14 @@ def main():  # noqa: D103
         writeEVENT(forces, filenameEVENT, floorsCount)
 
     else:
-        print('No RVs requested in StochasticWave.py')  # noqa: T201
+        print('No RVs requested')  # noqa: T201
         filenameEVENT = arguments.filenameEVENT  # noqa: N816
         result = subprocess.run(  # noqa: S603
             [  # noqa: S607
                 'ti',
-                f'{os.path.realpath(os.path.dirname(__file__))}'  # noqa: ISC003, PTH120
-                + '/taichi_script.py',
+                scriptName,
+                # f'{os.path.realpath(os.path.dirname(__file__))}'  # noqa: ISC003, PTH120
+                # + '/taichi_script.py',
             ],
             stdout=subprocess.PIPE,
             check=False,

modules/createEVENT/TaichiEvent/ad_gravity.py

@@ -0,0 +1,55 @@
+import taichi as ti
+
+ti.init()
+
+N = 8
+dt = 1e-5
+
+x = ti.Vector.field(2, dtype=ti.f32, shape=N, needs_grad=True)  # particle positions
+v = ti.Vector.field(2, dtype=ti.f32, shape=N)  # particle velocities
+U = ti.field(dtype=ti.f32, shape=(), needs_grad=True)  # potential energy
+
+
+@ti.kernel
+def compute_U():
+    for i, j in ti.ndrange(N, N):
+        r = x[i] - x[j]
+        # r.norm(1e-3) is equivalent to ti.sqrt(r.norm()**2 + 1e-3)
+        # This is to prevent 1/0 error which can cause wrong derivative
+        U[None] += -1 / r.norm(1e-3)  # U += -1 / |r|
+
+
+@ti.kernel
+def advance():
+    for i in x:
+        v[i] += dt * -x.grad[i]  # dv/dt = -dU/dx
+    for i in x:
+        x[i] += dt * v[i]  # dx/dt = v
+
+
+def substep():
+    with ti.ad.Tape(loss=U):
+        # Kernel invocations in this scope will later contribute to partial derivatives of
+        # U with respect to input variables such as x.
+        compute_U()  # The tape will automatically compute dU/dx and save the results in x.grad
+    advance()
+
+
+@ti.kernel
+def init():
+    for i in x:
+        x[i] = [ti.random(), ti.random()]
+
+
+def main():
+    init()
+    gui = ti.GUI("Autodiff gravity")
+    while gui.running:
+        for i in range(50):
+            substep()
+        gui.circles(x.to_numpy(), radius=3)
+        gui.show()
+
+
+if __name__ == "__main__":
+    main()

modules/createEVENT/TaichiEvent/comet.py

@@ -0,0 +1,105 @@
+import math
+
+import taichi as ti
+
+ti.init(arch=ti.cuda)
+
+dim = 3
+N = 1024 * 8
+dt = 2e-4
+steps = 7
+sun = ti.Vector([0.5, 0.5] if dim == 2 else [0.5, 0.5, 0.0])
+gravity = 0.5
+pressure = 0.3
+tail_paticle_scale = 0.4
+color_init = 0.3
+color_decay = 1.6
+vel_init = 0.07
+res = 640
+
+inv_m = ti.field(ti.f32)
+color = ti.field(ti.f32)
+x = ti.Vector.field(dim, ti.f32)
+v = ti.Vector.field(dim, ti.f32)
+ti.root.bitmasked(ti.i, N).place(x, v, inv_m, color)
+count = ti.field(ti.i32, ())
+img = ti.field(ti.f32, (res, res))
+
+
+@ti.func
+def rand_unit_2d():
+    a = ti.random() * 2 * math.pi
+    return ti.Vector([ti.cos(a), ti.sin(a)])
+
+
+@ti.func
+def rand_unit_3d():
+    u = rand_unit_2d()
+    s = ti.random() * 2 - 1
+    c = ti.sqrt(1 - s**2)
+    return ti.Vector([c * u[0], c * u[1], s])
+
+
+@ti.kernel
+def substep():
+    ti.no_activate(x)
+    for i in x:
+        r = x[i] - sun
+        r_sq_inverse = r / r.norm(1e-3) ** 3
+        acceleration = (pressure * inv_m[i] - gravity) * r_sq_inverse
+        v[i] += acceleration * dt
+        x[i] += v[i] * dt
+        color[i] *= ti.exp(-dt * color_decay)
+
+        if not all(-0.1 <= x[i] <= 1.1):
+            ti.deactivate(x.snode.parent(), [i])
+
+
+@ti.kernel
+def generate():
+    r = x[0] - sun
+    n_tail_paticles = int(tail_paticle_scale / r.norm(1e-3) ** 2)
+    for _ in range(n_tail_paticles):
+        r = x[0]
+        if ti.static(dim == 3):
+            r = rand_unit_3d()
+        else:
+            r = rand_unit_2d()
+        xi = ti.atomic_add(count[None], 1) % (N - 1) + 1
+        x[xi] = x[0]
+        v[xi] = r * vel_init + v[0]
+        inv_m[xi] = 0.5 + ti.random()
+        color[xi] = color_init
+
+
+@ti.kernel
+def render():
+    for p in ti.grouped(img):
+        img[p] = 1e-6 / (p / res - ti.Vector([sun.x, sun.y])).norm(1e-4) ** 3
+    for i in x:
+        p = int(ti.Vector([x[i].x, x[i].y]) * res)
+        if 0 <= p[0] < res and 0 <= p[1] < res:
+            img[p] += color[i]
+
+
+def main():
+    inv_m[0] = 0
+    x[0].x = +0.5
+    x[0].y = -0.01
+    v[0].x = +0.6
+    v[0].y = +0.4
+    color[0] = 1
+
+    gui = ti.GUI("Comet", res)
+    while gui.running:
+        gui.running = not gui.get_event(gui.ESCAPE)
+        generate()
+        for s in range(steps):
+            substep()
+        render()
+        gui.set_image(img)
+        gui.show()
+
+
+if __name__ == "__main__":
+    main()