diff --git a/appveyor.yml b/appveyor.yml index 10253539b..00af4b0ac 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -1,4 +1,4 @@ -version: 24.07.{build} +version: 24.09.{build} image: - Visual Studio 2019 @@ -26,7 +26,7 @@ for: install: - uname - rm -rf SimCenterBackendApplications - - git clone https://github.com/NHERI-SimCenter/SimCenterBackendApplications.git + - git clone https://github.com/JustinBonus/SimCenterBackendApplications.git build_script: @@ -63,7 +63,7 @@ for: - conan user - conan remote add simcenter https://nherisimcenter.jfrog.io/artifactory/api/conan/simcenter - rm -rf SimCenterBackendApplications - - git clone https://github.com/NHERI-SimCenter/SimCenterBackendApplications.git + - git clone https://github.com/JustinBonus/SimCenterBackendApplications.git build_script: # build SimCenterBackendApplications @@ -104,7 +104,7 @@ for: install: - cmd: rm -rf SimCenterBackendApplications - - cmd: git clone https://github.com/NHERI-SimCenter/SimCenterBackendApplications.git + - cmd: git clone https://github.com/JustinBonus/SimCenterBackendApplications.git - cmd: dir build_script: diff --git a/makeEXE.sh b/makeEXE.sh index 50f9ebd54..9cb613d21 100755 --- a/makeEXE.sh +++ b/makeEXE.sh @@ -1,11 +1,11 @@ #!/bin/bash -echo "Building SimCenterBackendApplicationss ..." +echo "Building SimCenterBackendApplications ..." mkdir -p build cd build conan install .. --build missing status=$?; if [[ $status != 0 ]]; then echo "conan failed"; exit $status; fi -cmake .. +cmake -DCMAKE_CXX_STANDARD=17 -DCMAKE_CXX_STANDARD_REQUIRED=ON -DCMAKE_CXX_EXTENSIONS=OFF .. status=$?; if [[ $status != 0 ]]; then echo "cmake failed"; exit $status; fi cmake --build . --config Release status=$?; if [[ $status != 0 ]]; then echo "make failed"; exit $status; fi diff --git a/modules/createEVENT/TaichiEvent/TaichiEvent.py b/modules/createEVENT/TaichiEvent/TaichiEvent.py index 0357198a1..b81435fe7 100755 --- a/modules/createEVENT/TaichiEvent/TaichiEvent.py +++ b/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, diff --git a/modules/createEVENT/TaichiEvent/ad_gravity.py b/modules/createEVENT/TaichiEvent/ad_gravity.py new file mode 100644 index 000000000..1b22347d6 --- /dev/null +++ b/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() diff --git a/modules/createEVENT/TaichiEvent/comet.py b/modules/createEVENT/TaichiEvent/comet.py new file mode 100644 index 000000000..789711928 --- /dev/null +++ b/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() diff --git a/modules/createEVENT/TaichiEvent/euler.py b/modules/createEVENT/TaichiEvent/euler.py new file mode 100644 index 000000000..647293bfd --- /dev/null +++ b/modules/createEVENT/TaichiEvent/euler.py @@ -0,0 +1,461 @@ +from matplotlib import cm + +import taichi as ti + +# A compressible euler equation solver using two methods +# 1: 2nd order muscl +# 2: thinc BVD, ref: "Limiter-free discontinuity-capturing scheme +# for compressible gas dynamics with reactive fronts" + +real = ti.f32 +ti.init(arch=ti.gpu, default_fp=real) + +N = 512 # grid resolution +CFL = 0.9 # keep below 1 +method = 1 # 0:muscl, 1:thinc +IC_type = 0 # 0:sod +BC_type = 0 # 0:walls + +img_field = 0 # 0:density, 1: schlieren, 2:vorticity, 3: velocity mag +res = 512 # gui resolution +cmap_name = "magma_r" # python colormap +use_fixed_caxis = 0 # 1: use fixed caxis limits, 0: automatic caxis limits +fixed_caxis = [0.0, 5.0] # fixed caxis limits + +Q = ti.Vector.field(4, dtype=real, shape=(N, N)) # [rho, rho*u, rho*v, rho*e] consv vars +Q_old = ti.Vector.field(4, dtype=real, shape=(N, N)) +W = ti.Vector.field(4, dtype=real, shape=(N, N)) # [rho, u, v, p] cell avg +W_xl = ti.Vector.field(4, dtype=real, shape=(N, N, 3)) # left side of x-face +W_xr = ti.Vector.field(4, dtype=real, shape=(N, N, 3)) # right side of x-face +W_yl = ti.Vector.field(4, dtype=real, shape=(N, N, 3)) # left side of y-face +W_yr = ti.Vector.field(4, dtype=real, shape=(N, N, 3)) # right side of y-face +F_x = ti.Vector.field(4, dtype=real, shape=(N, N)) # x-face flux +F_y = ti.Vector.field(4, dtype=real, shape=(N, N)) # y-face flux +dt = ti.field(dtype=real, shape=()) +img = ti.field(dtype=ti.f32, shape=(res, res)) + +beta_smooth = 1.2 +beta_sharp = 2.0 +gamma = 1.4 # ratio of specific heats +h = 1.0 / (N - 2) # cell size +vol = h * h # cell volume + + +@ti.func +def is_interior_cell(i, j): + return 0 < i < N - 1 and 0 < j < N - 1 + + +@ti.func +def is_interior_x_face(i, j): + return 1 < i < N - 1 and 0 < j < N - 1 + + +@ti.func +def is_boundary_x_face(i, j): + return (i == 1 or i == N - 1) and 0 < j < N - 1 + + +@ti.func +def is_interior_y_face(i, j): + return 0 < i < N - 1 and 1 < j < N - 1 + + +@ti.func +def is_boundary_y_face(i, j): + return 0 < i < N - 1 and (j == 1 or j == N - 1) + + +@ti.func +def get_cell_pos(i, j): + return ti.Vector([i * h - h / 2.0, j * h - h / 2.0]) + + +@ti.kernel +def compute_W(): + # conversion from conservative variables to primitive variables + for i, j in Q: + W[i, j] = q_to_w(Q[i, j]) + + +@ti.kernel +def copy_to_old(): + for i, j in Q: + Q_old[i, j] = Q[i, j] + + +@ti.kernel +def set_ic(): + for i, j in Q: + if IC_type == 0: + # primitive variable initial conditions + w_in = ti.Vector([10.0, 0.0, 0.0, 10.0]) + w_out = ti.Vector([0.125, 0.0, 0.0, 0.1]) + + pos = get_cell_pos(i, j) + center = ti.Vector([0.5, 0.5]) + + if (pos - center).norm() < 0.25: + Q[i, j] = w_to_q(w_in) + else: + Q[i, j] = w_to_q(w_out) + + # implement more ic's later + + +@ti.kernel +def set_bc(): + # enforce boundary conditions by setting ghost cells + for i, j in Q: + if not is_interior_cell(i, j): + if BC_type == 0: # walls + # enforce neumann=0 and zero normal velocity on face + if i == 0: + Q[i, j] = Q[i + 1, j] + Q[i, j][1] = -Q[i + 1, j][1] + if i == N - 1: + Q[i, j] = Q[i - 1, j] # neumann 0 bc + Q[i, j][1] = -Q[i - 1, j][1] # enforce 0 normal velocty at face + if j == 0: + Q[i, j] = Q[i, j + 1] + Q[i, j][2] = -Q[i, j + 1][2] + if j == N - 1: + Q[i, j] = Q[i, j - 1] + Q[i, j][2] = -Q[i, j - 1][2] + + # implement more bc's later + + +@ti.func +def mc_lim(r): + # MC flux limiter + return max(0.0, min(2.0 * r, min(0.5 * (r + 1.0), 2.0))) + + +@ti.func +def w_to_q(w): + # convert primitive variables to conserved variables + q = ti.Vector([0.0, 0.0, 0.0, 0.0]) + q[0] = w[0] # rho + q[1] = w[0] * w[1] # rho*u + q[2] = w[0] * w[2] # rho*v + q[3] = w[0] * (w[3] / ((gamma - 1) * w[0]) + 0.5 * (w[1] ** 2 + w[2] ** 2)) + # rho*e + return q + + +@ti.func +def q_to_w(q): + # convert conserved variables to primitive variables + w = ti.Vector([0.0, 0.0, 0.0, 0.0]) + w[0] = q[0] # rho + w[1] = q[1] / q[0] # u + w[2] = q[2] / q[0] # v + w[3] = (gamma - 1) * (q[3] - 0.5 * (q[1] ** 2 + q[2] ** 2) / q[0]) + # p + return w + + +@ti.func +def HLLC_flux(qL, qR, n): + # normal vector + nx = n[0] + ny = n[1] + + # Left state + rL = qL[0] # rho + uL = qL[1] / qL[0] # u + vL = qL[2] / qL[0] # v + pL = (gamma - 1.0) * (qL[3] - 0.5 * (qL[1] ** 2 + qL[2] ** 2) / qL[0]) + # p + vnL = uL * nx + vL * ny + vtL = -uL * ny + vL * nx + aL = ti.sqrt(gamma * pL / rL) + HL = (qL[3] + pL) / rL + + # Right state + rR = qR[0] # rho + uR = qR[1] / qR[0] # u + vR = qR[2] / qR[0] # v + pR = (gamma - 1.0) * (qR[3] - 0.5 * (qR[1] ** 2 + qR[2] ** 2) / qR[0]) + # p + vnR = uR * nx + vR * ny + vtR = -uR * ny + vR * nx + aR = ti.sqrt(gamma * pR / rR) + HR = (qR[3] + pR) / rR + + # Left and Right fluxes + fL = ti.Vector([rL * vnL, rL * vnL * uL + pL * nx, rL * vnL * vL + pL * ny, rL * vnL * HL]) + fR = ti.Vector([rR * vnR, rR * vnR * uR + pR * nx, rR * vnR * vR + pR * ny, rR * vnR * HR]) + + # Roe Averages + rt = ti.sqrt(rR / rL) + u = (uL + rt * uR) / (1.0 + rt) + v = (vL + rt * vR) / (1.0 + rt) + H = (HL + rt * HR) / (1.0 + rt) + a = ti.sqrt((gamma - 1.0) * (H - (u**2 + v**2) / 2.0)) + vn = u * nx + v * ny + + # wavespeeds + sL = min(vnL - aL, vn - a) + sR = max(vnR + aR, vn + a) + sM = (pL - pR + rR * vnR * (sR - vnR) - rL * vnL * (sL - vnL)) / (rR * (sR - vnR) - rL * (sL - vnL)) + + # HLLC flux. + HLLC = ti.Vector([0.0, 0.0, 0.0, 0.0]) + if 0 <= sL: + HLLC = fL + elif 0 <= sM: + qsL = ( + rL + * (sL - vnL) + / (sL - sM) + * ti.Vector( + [ + 1.0, + sM * nx - vtL * ny, + sM * ny + vtL * nx, + qL[3] / rL + (sM - vnL) * (sM + pL / (rL * (sL - vnL))), + ] + ) + ) + HLLC = fL + sL * (qsL - qL) + elif 0 <= sR: + qsR = ( + rR + * (sR - vnR) + / (sR - sM) + * ti.Vector( + [ + 1.0, + sM * nx - vtR * ny, + sM * ny + vtR * nx, + qR[3] / rR + (sM - vnR) * (sM + pR / (rR * (sR - vnR))), + ] + ) + ) + HLLC = fR + sR * (qsR - qR) + elif 0 >= sR: + HLLC = fR + + return HLLC + + +@ti.kernel +def compute_F_muscl(): + for i, j in Q: + if is_interior_x_face(i, j): + # muscl reconstrucion of left and right states with HLLC flux + wL = ti.Vector([0.0, 0.0, 0.0, 0.0]) + wR = ti.Vector([0.0, 0.0, 0.0, 0.0]) + for f in ti.static(range(4)): + ratio_l = (W[i, j][f] - W[i - 1, j][f]) / (W[i - 1, j][f] - W[i - 2, j][f]) + ratio_r = (W[i, j][f] - W[i - 1, j][f]) / (W[i + 1, j][f] - W[i, j][f]) + wL[f] = W[i - 1, j][f] + 0.5 * mc_lim(ratio_l) * (W[i - 1, j][f] - W[i - 2, j][f]) + wR[f] = W[i, j][f] - 0.5 * mc_lim(ratio_r) * (W[i + 1, j][f] - W[i, j][f]) + F_x[i, j] = HLLC_flux(w_to_q(wL), w_to_q(wR), ti.Vector([1.0, 0.0])) + + elif is_boundary_x_face(i, j): + F_x[i, j] = HLLC_flux(Q[i - 1, j], Q[i, j], ti.Vector([1.0, 0.0])) + + if is_interior_y_face(i, j): + # muscl reconstrucion of left and right states with HLLC flux + wL = ti.Vector([0.0, 0.0, 0.0, 0.0]) + wR = ti.Vector([0.0, 0.0, 0.0, 0.0]) + for f in ti.static(range(4)): + ratio_l = (W[i, j][f] - W[i, j - 1][f]) / (W[i, j - 1][f] - W[i, j - 2][f]) + ratio_r = (W[i, j][f] - W[i, j - 1][f]) / (W[i, j + 1][f] - W[i, j][f]) + wL[f] = W[i, j - 1][f] + 0.5 * mc_lim(ratio_l) * (W[i, j - 1][f] - W[i, j - 2][f]) + wR[f] = W[i, j][f] - 0.5 * mc_lim(ratio_r) * (W[i, j + 1][f] - W[i, j][f]) + F_y[i, j] = HLLC_flux(w_to_q(wL), w_to_q(wR), ti.Vector([0.0, 1.0])) + + elif is_boundary_y_face(i, j): + F_y[i, j] = HLLC_flux(Q[i, j - 1], Q[i, j], ti.Vector([0.0, 1.0])) + + +@ti.func +def sign(a): + sgn = 0.0 + if a > 0.0: + sgn = 1.0 + elif a < 0.0: + sgn = -1.0 + return sgn + + +@ti.func +def cosh(a): + return (ti.exp(a) + ti.exp(-a)) / 2.0 + + +@ti.func +def thinc(wl, wc, wr, beta): + w0 = wc + w1 = wc + if (wr - wc) * (wc - wl) > 0.0: + # use thinc reconstruction + eps = 1.0e-15 + wmin = ti.min(wr, wl) + wmax = ti.max(wr, wl) + wdelta = wmax - wmin + theta = sign(wr - wl) + C = (wc - wmin + eps) / (wdelta + eps) + B = ti.exp(theta * beta * (2 * C - 1)) + A = (B / cosh(beta) - 1) / ti.tanh(beta) + + # reconstructed value on right side of left face + w0 = wmin + wdelta / 2.0 * (1.0 + theta * A) + + # reconstructed value on left side of right face + w1 = wmin + wdelta / 2.0 * (1.0 + theta * (ti.tanh(beta) + A) / (1.0 + A * ti.tanh(beta))) + + return w0, w1 + + +@ti.kernel +def compute_F_thinc(): + # reconstruct primitive variables on interior faces of each cell using + # multiple candidate thinc reconstructions + for i, j in Q: + if is_interior_cell(i, j): + for f in ti.static(range(4)): + # smooth x-dir reconstruction + w0, w1 = thinc(W[i - 1, j][f], W[i, j][f], W[i + 1, j][f], beta_smooth) + W_xr[i, j, 0][f] = w0 + W_xl[i + 1, j, 0][f] = w1 + + # sharp x-dir reconstruction + w0, w1 = thinc(W[i - 1, j][f], W[i, j][f], W[i + 1, j][f], beta_sharp) + W_xr[i, j, 1][f] = w0 + W_xl[i + 1, j, 1][f] = w1 + + # smooth y-dir reconstruction + w0, w1 = thinc(W[i, j - 1][f], W[i, j][f], W[i, j + 1][f], beta_smooth) + W_yr[i, j, 0][f] = w0 + W_yl[i, j + 1, 0][f] = w1 + + # sharp y-dir reconstruction + w0, w1 = thinc(W[i, j - 1][f], W[i, j][f], W[i, j + 1][f], beta_sharp) + W_yr[i, j, 1][f] = w0 + W_yl[i, j + 1, 1][f] = w1 + + for i, j in Q: + # choose the final reconstruction for each cell using the BVD algorithm + if is_interior_cell(i, j): + for f in ti.static(range(4)): + # x-dir + TBV_smooth = abs(W_xl[i, j, 0][f] - W_xr[i, j, 0][f]) + abs(W_xl[i + 1, j, 0][f] - W_xr[i + 1, j, 0][f]) + TBV_sharp = abs(W_xl[i, j, 1][f] - W_xr[i, j, 1][f]) + abs(W_xl[i + 1, j, 1][f] - W_xr[i + 1, j, 1][f]) + + if TBV_smooth < TBV_sharp: + W_xr[i, j, 2][f] = W_xr[i, j, 0][f] + W_xl[i + 1, j, 2][f] = W_xl[i + 1, j, 0][f] + else: + W_xr[i, j, 2][f] = W_xr[i, j, 1][f] + W_xl[i + 1, j, 2][f] = W_xl[i + 1, j, 1][f] + + # y-dir + TBV_smooth = abs(W_yl[i, j, 0][f] - W_yr[i, j, 0][f]) + abs(W_yl[i, j + 1, 0][f] - W_yr[i, j + 1, 0][f]) + TBV_sharp = abs(W_yl[i, j, 1][f] - W_yr[i, j, 1][f]) + abs(W_yl[i, j + 1, 1][f] - W_yr[i, j + 1, 1][f]) + + if TBV_smooth < TBV_sharp: + W_yr[i, j, 2][f] = W_yr[i, j, 0][f] + W_yl[i, j + 1, 2][f] = W_yl[i, j + 1, 0][f] + else: + W_yr[i, j, 2][f] = W_yr[i, j, 1][f] + W_yl[i, j + 1, 2][f] = W_yl[i, j + 1, 1][f] + + for i, j in Q: + # compute numerical fluxes of with Riemann solver + if is_interior_x_face(i, j): + # muscl reconstrucion of left and right states with HLLC flux + F_x[i, j] = HLLC_flux(w_to_q(W_xl[i, j, 2]), w_to_q(W_xr[i, j, 2]), ti.Vector([1.0, 0.0])) + elif is_boundary_x_face(i, j): + F_x[i, j] = HLLC_flux(Q[i - 1, j], Q[i, j], ti.Vector([1.0, 0.0])) + + if is_interior_y_face(i, j): + F_y[i, j] = HLLC_flux(w_to_q(W_yl[i, j, 2]), w_to_q(W_yr[i, j, 2]), ti.Vector([0.0, 1.0])) + elif is_boundary_y_face(i, j): + F_y[i, j] = HLLC_flux(Q[i, j - 1], Q[i, j], ti.Vector([0.0, 1.0])) + + +@ti.kernel +def calc_dt(): + dt[None] = 1.0e5 + for i, j in Q: + w = q_to_w(Q[i, j]) + a = ti.sqrt(gamma * w[3] / w[0]) + vel = ti.sqrt(w[1] ** 2 + w[2] ** 2) + ws = a + vel + ti.atomic_min(dt[None], CFL * h / ws / 2.0) + + +@ti.kernel +def update_Q(rk_step: ti.template()): + for i, j in Q: + if is_interior_cell(i, j): + if ti.static(rk_step == 0): + Q[i, j] = Q[i, j] + dt[None] * (F_x[i, j] - F_x[i + 1, j] + F_y[i, j] - F_y[i, j + 1]) / h + if ti.static(rk_step == 1): + Q[i, j] = (Q[i, j] + Q_old[i, j]) / 2.0 + dt[None] * ( + F_x[i, j] - F_x[i + 1, j] + F_y[i, j] - F_y[i, j + 1] + ) / h + + +@ti.kernel +def paint(): + for i, j in img: + ii = min(max(1, i * N // res), N - 2) + jj = min(max(1, j * N // res), N - 2) + if img_field == 0: # density + img[i, j] = Q[ii, jj][0] + elif img_field == 1: # numerical schlieren + img[i, j] = ti.sqrt( + ((Q[ii + 1, jj][0] - Q[ii - 1, jj][0]) / h) ** 2 + ((Q[ii, jj + 1][0] - Q[ii, jj - 1][0]) / h) ** 2 + ) + elif img_field == 2: # vorticity + img[i, j] = (Q[ii + 1, jj][2] - Q[ii - 1, jj][2]) / h - (Q[ii, jj + 1][1] - Q[ii, jj - 1][1]) / h + elif img_field == 3: # velocity magnitude + img[i, j] = ti.sqrt(Q[ii, jj][1] ** 2 + Q[ii, jj][2] ** 2) + + max_ = -1.0e10 + min_ = 1.0e10 + for i, j in img: + ti.atomic_max(max_, img[i, j]) + ti.atomic_min(min_, img[i, j]) + + for i, j in img: + if use_fixed_caxis: + min_ = fixed_caxis[0] + max_ = fixed_caxis[1] + img[i, j] = (img[i, j] - min_) / (max_ - min_) + + +def main(): + gui = ti.GUI("Euler Equations", (res, res)) + cmap = cm.get_cmap(cmap_name) + set_ic() + set_bc() + + n = 0 + while gui.running: + calc_dt() + copy_to_old() + for rk_step in range(2): + compute_W() + if method == 0: + compute_F_muscl() + else: + compute_F_thinc() + update_Q(rk_step) + set_bc() + + if n % 10 == 0: + paint() + gui.set_image(cmap(img.to_numpy())) + gui.show() + n += 1 + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/eulerfluid2d.py b/modules/createEVENT/TaichiEvent/eulerfluid2d.py new file mode 100644 index 000000000..216755e3a --- /dev/null +++ b/modules/createEVENT/TaichiEvent/eulerfluid2d.py @@ -0,0 +1,310 @@ +# 2D Euler Fluid Simulation using Taichi, originally created by @Lee-abcde +from taichi.examples.patterns import taichi_logo + +import taichi as ti +import taichi.math as tm + +ti.init(arch=ti.gpu) + + +##################### +# Bilinear Interpolation function +##################### +@ti.func +def sample(vf, u, v, shape): + i, j = int(u), int(v) + # Nearest + i = ti.max(0, ti.min(shape[0] - 1, i)) + j = ti.max(0, ti.min(shape[1] - 1, j)) + return vf[i, j] + + +@ti.func +def lerp(vl, vr, frac): + # frac: [0.0, 1.0] + return (1 - frac) * vl + frac * vr + + +@ti.func +def bilerp(vf, u, v, shape): + # use -0.5 to decide where bilerp performs in cells + s, t = u - 0.5, v - 0.5 + iu, iv = int(s), int(t) + a = sample(vf, iu + 0.5, iv + 0.5, shape) + b = sample(vf, iu + 1.5, iv + 0.5, shape) + c = sample(vf, iu + 0.5, iv + 1.5, shape) + d = sample(vf, iu + 1.5, iv + 1.5, shape) + # fract + fu, fv = s - iu, t - iv + return lerp(lerp(a, b, fu), lerp(c, d, fu), fv) + + +##################### +# Simulation parameters +##################### +eulerSimParam = { + "shape": [512, 512], + "dt": 1 / 60.0, + "iteration_step": 20, + "mouse_radius": 0.01, # [0.0,1.0] float + "mouse_speed": 125.0, + "mouse_respondDistance": 0.5, # for every frame, only half the trace of the mouse will influence water + "curl_param": 15, +} + + +# Double Buffer +class TexPair: + def __init__(self, cur, nxt): + self.cur = cur + self.nxt = nxt + + def swap(self): + self.cur, self.nxt = self.nxt, self.cur + + +velocityField = ti.Vector.field(2, float, shape=(eulerSimParam["shape"])) +_new_velocityField = ti.Vector.field(2, float, shape=(eulerSimParam["shape"])) +colorField = ti.Vector.field(3, float, shape=(eulerSimParam["shape"])) +_new_colorField = ti.Vector.field(3, float, shape=(eulerSimParam["shape"])) + +curlField = ti.field(float, shape=(eulerSimParam["shape"])) + +divField = ti.field(float, shape=(eulerSimParam["shape"])) +pressField = ti.field(float, shape=(eulerSimParam["shape"])) +_new_pressField = ti.field(float, shape=(eulerSimParam["shape"])) + +velocities_pair = TexPair(velocityField, _new_velocityField) +pressure_pair = TexPair(pressField, _new_pressField) +color_pair = TexPair(colorField, _new_colorField) + + +@ti.kernel +def init_field(): + # init pressure and velocity fieldfield + pressField.fill(0) + velocityField.fill(0) + for i, j in ti.ndrange(eulerSimParam["shape"][0] * 4, eulerSimParam["shape"][1] * 4): + # 4x4 super sampling: + ret = taichi_logo(ti.Vector([i, j]) / (eulerSimParam["shape"][0] * 4)) + colorField[i // 4, j // 4] += ret / 16 + + +@ti.kernel +def advection(vf: ti.template(), qf: ti.template(), new_qf: ti.template()): + for i, j in vf: + coord_cur = ti.Vector([i, j]) + ti.Vector([0.5, 0.5]) + vel_cur = vf[i, j] + coord_prev = coord_cur - vel_cur * eulerSimParam["dt"] + q_prev = bilerp(qf, coord_prev[0], coord_prev[1], (eulerSimParam["shape"])) + new_qf[i, j] = q_prev + + +@ti.kernel +def curl(vf: ti.template(), cf: ti.template()): + for i, j in vf: + cf[i, j] = 0.5 * ((vf[i + 1, j][1] - vf[i - 1, j][1]) - (vf[i, j + 1][0] - vf[i, j - 1][0])) + + +@ti.kernel +def vorticity_projection(cf: ti.template(), vf: ti.template(), vf_new: ti.template()): + for i, j in cf: + gradcurl = ti.Vector( + [ + 0.5 * (cf[i + 1, j] - cf[i - 1, j]), + 0.5 * (cf[i, j + 1] - cf[i, j - 1]), + 0, + ] + ) + GradCurlLength = tm.length(gradcurl) + if GradCurlLength > 1e-5: + force = eulerSimParam["curl_param"] * tm.cross(gradcurl / GradCurlLength, ti.Vector([0, 0, 1])) + vf_new[i, j] = vf[i, j] + eulerSimParam["dt"] * force[:2] + + +@ti.kernel +def divergence(vf: ti.template(), divf: ti.template()): + for i, j in vf: + divf[i, j] = 0.5 * (vf[i + 1, j][0] - vf[i - 1, j][0] + vf[i, j + 1][1] - vf[i, j - 1][1]) + + +@ti.kernel +def pressure_iteration(divf: ti.template(), pf: ti.template(), new_pf: ti.template()): + for i, j in pf: + new_pf[i, j] = (pf[i + 1, j] + pf[i - 1, j] + pf[i, j - 1] + pf[i, j + 1] - divf[i, j]) / 4 + + +def pressure_solve(presspair: TexPair, divf: ti.template()): + for i in range(eulerSimParam["iteration_step"]): + pressure_iteration(divf, presspair.cur, presspair.nxt) + presspair.swap() + apply_p_bc(presspair.cur) + + +@ti.kernel +def pressure_projection(pf: ti.template(), vf: ti.template(), vf_new: ti.template()): + for i, j in vf: + vf_new[i, j] = vf[i, j] - ti.Vector( + [ + ( + p_with_boundary(pf, i + 1, j, eulerSimParam["shape"]) + - p_with_boundary(pf, i - 1, j, eulerSimParam["shape"]) + ) + / 2.0, + ( + p_with_boundary(pf, i, j + 1, eulerSimParam["shape"]) + - p_with_boundary(pf, i, j - 1, eulerSimParam["shape"]) + ) + / 2.0, + ] + ) + + +##################### +# Boundary Condition +##################### +@ti.func +def vel_with_boundary(vf: ti.template(), i: int, j: int, shape) -> ti.f32: + if (i <= 0) or (i >= shape[0] - 1) or (j >= shape[1] - 1) or (j <= 0): + vf[i, j] = ti.Vector([0.0, 0.0]) + return vf[i, j] + + +@ti.func +def p_with_boundary(pf: ti.template(), i: int, j: int, shape) -> ti.f32: + if ( + (i == j == 0) + or (i == shape[0] - 1 and j == shape[1] - 1) + or (i == 0 and j == shape[1] - 1) + or (i == shape[0] - 1 and j == 0) + ): + pf[i, j] = 0.0 + elif i == 0: + pf[0, j] = pf[1, j] + elif j == 0: + pf[i, 0] = pf[i, 1] + elif i == shape[0] - 1: + pf[shape[0] - 1, j] = pf[shape[0] - 2, j] + elif j == shape[1] - 1: + pf[i, shape[1] - 1] = pf[i, shape[1] - 2] + return pf[i, j] + + +@ti.func +def c_with_boundary(cf: ti.template(), i: int, j: int, shape) -> ti.f32: + if (i <= 0) or (i >= shape[0] - 1) or (j >= shape[1] - 1) or (j <= 0): + cf[i, j] = 0.0 + return cf[i, j] + + +@ti.kernel +def apply_vel_bc(vf: ti.template()): + for i, j in vf: + vel_with_boundary(vf, i, j, eulerSimParam["shape"]) + + +@ti.kernel +def apply_p_bc(pf: ti.template()): + for i, j in pf: + p_with_boundary(pf, i, j, eulerSimParam["shape"]) + + +@ti.kernel +def apply_c_bc(cf: ti.template()): + for i, j in cf: + c_with_boundary(cf, i, j, eulerSimParam["shape"]) + + +def mouse_interaction(prev_posx: int, prev_posy: int): + mouse_x, mouse_y = window.get_cursor_pos() + mousePos_x = int(mouse_x * eulerSimParam["shape"][0]) + mousePos_y = int(mouse_y * eulerSimParam["shape"][1]) + if prev_posx == 0 and prev_posy == 0: + prev_posx = mousePos_x + prev_posy = mousePos_y + mouseRadius = eulerSimParam["mouse_radius"] * min(eulerSimParam["shape"][0], eulerSimParam["shape"][1]) + + mouse_addspeed( + mousePos_x, + mousePos_y, + prev_posx, + prev_posy, + mouseRadius, + velocities_pair.cur, + velocities_pair.nxt, + ) + velocities_pair.swap() + prev_posx = mousePos_x + prev_posy = mousePos_y + return prev_posx, prev_posy + + +@ti.kernel +def mouse_addspeed( + cur_posx: int, + cur_posy: int, + prev_posx: int, + prev_posy: int, + mouseRadius: float, + vf: ti.template(), + new_vf: ti.template(), +): + for i, j in vf: + vec1 = ti.Vector([cur_posx - prev_posx, cur_posy - prev_posy]) + vec2 = ti.Vector([i - prev_posx, j - prev_posy]) + dotans = tm.dot(vec1, vec2) + distance = abs(tm.cross(vec1, vec2)) / (tm.length(vec1) + 0.001) + if ( + dotans >= 0 + and dotans <= eulerSimParam["mouse_respondDistance"] * tm.length(vec1) + and distance <= mouseRadius + ): + new_vf[i, j] = vf[i, j] + vec1 * eulerSimParam["mouse_speed"] + else: + new_vf[i, j] = vf[i, j] + + +def advaction_step(): + advection(velocities_pair.cur, color_pair.cur, color_pair.nxt) + advection(velocities_pair.cur, velocities_pair.cur, velocities_pair.nxt) + color_pair.swap() + velocities_pair.swap() + apply_vel_bc(velocities_pair.cur) + + +def voricity_step(): + curl(velocities_pair.cur, curlField) + apply_c_bc(curlField) + vorticity_projection(curlField, velocities_pair.cur, velocities_pair.nxt) + velocities_pair.swap() + apply_vel_bc(velocities_pair.cur) + + +def pressure_step(): + divergence(velocities_pair.cur, divField) + pressure_solve(pressure_pair, divField) + pressure_projection(pressure_pair.cur, velocities_pair.cur, velocities_pair.nxt) + velocities_pair.swap() + apply_vel_bc(velocities_pair.cur) + + +##################### +# The Euler Simulation starts! +##################### +init_field() +apply_vel_bc(velocities_pair.cur) +window = ti.GUI("Euler 2D Simulation", res=(eulerSimParam["shape"][0], eulerSimParam["shape"][1])) + +mouse_prevposx, mouse_prevposy = 0, 0 +while window.running: + # advection + advaction_step() + # curl + voricity_step() + # mouse interact with fluid + mouse_prevposx, mouse_prevposy = mouse_interaction(mouse_prevposx, mouse_prevposy) + # pressure iteration and projection + pressure_step() + + window.set_image(colorField) + window.show() diff --git a/modules/createEVENT/TaichiEvent/fem128.py b/modules/createEVENT/TaichiEvent/fem128.py new file mode 100644 index 000000000..3f8d3d8d7 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/fem128.py @@ -0,0 +1,144 @@ +import taichi as ti + +ti.init(arch=ti.gpu) + +N = 12 +dt = 5e-5 +dx = 1 / N +rho = 4e1 +NF = 2 * N**2 # number of faces +NV = (N + 1) ** 2 # number of vertices +E, nu = 4e4, 0.2 # Young's modulus and Poisson's ratio +mu, lam = E / 2 / (1 + nu), E * nu / (1 + nu) / (1 - 2 * nu) # Lame parameters +ball_pos, ball_radius = ti.Vector([0.5, 0.0]), 0.31 +damping = 14.5 + +pos = ti.Vector.field(2, float, NV, needs_grad=True) +vel = ti.Vector.field(2, float, NV) +f2v = ti.Vector.field(3, int, NF) # ids of three vertices of each face +B = ti.Matrix.field(2, 2, float, NF) +F = ti.Matrix.field(2, 2, float, NF, needs_grad=True) +V = ti.field(float, NF) +phi = ti.field(float, NF) # potential energy of each face (Neo-Hookean) +U = ti.field(float, (), needs_grad=True) # total potential energy + +gravity = ti.Vector.field(2, float, ()) +attractor_pos = ti.Vector.field(2, float, ()) +attractor_strength = ti.field(float, ()) + + +@ti.kernel +def update_U(): + for i in range(NF): + ia, ib, ic = f2v[i] + a, b, c = pos[ia], pos[ib], pos[ic] + V[i] = abs((a - c).cross(b - c)) + D_i = ti.Matrix.cols([a - c, b - c]) + F[i] = D_i @ B[i] + for i in range(NF): + F_i = F[i] + log_J_i = ti.log(F_i.determinant()) + phi_i = mu / 2 * ((F_i.transpose() @ F_i).trace() - 2) + phi_i -= mu * log_J_i + phi_i += lam / 2 * log_J_i**2 + phi[i] = phi_i + U[None] += V[i] * phi_i + + +@ti.kernel +def advance(): + for i in range(NV): + acc = -pos.grad[i] / (rho * dx**2) + g = gravity[None] * 0.8 + attractor_strength[None] * (attractor_pos[None] - pos[i]).normalized(1e-5) + vel[i] += dt * (acc + g * 40) + vel[i] *= ti.exp(-dt * damping) + for i in range(NV): + # ball boundary condition: + disp = pos[i] - ball_pos + disp2 = disp.norm_sqr() + if disp2 <= ball_radius**2: + NoV = vel[i].dot(disp) + if NoV < 0: + vel[i] -= NoV * disp / disp2 + cond = (pos[i] < 0) & (vel[i] < 0) | (pos[i] > 1) & (vel[i] > 0) + # rect boundary condition: + for j in ti.static(range(pos.n)): + if cond[j]: + vel[i][j] = 0 + pos[i] += dt * vel[i] + + +@ti.kernel +def init_pos(): + for i, j in ti.ndrange(N + 1, N + 1): + k = i * (N + 1) + j + pos[k] = ti.Vector([i, j]) / N * 0.25 + ti.Vector([0.45, 0.45]) + vel[k] = ti.Vector([0, 0]) + for i in range(NF): + ia, ib, ic = f2v[i] + a, b, c = pos[ia], pos[ib], pos[ic] + B_i_inv = ti.Matrix.cols([a - c, b - c]) + B[i] = B_i_inv.inverse() + + +@ti.kernel +def init_mesh(): + for i, j in ti.ndrange(N, N): + k = (i * N + j) * 2 + a = i * (N + 1) + j + b = a + 1 + c = a + N + 2 + d = a + N + 1 + f2v[k + 0] = [a, b, c] + f2v[k + 1] = [c, d, a] + + +def paint_phi(gui): + pos_ = pos.to_numpy() + phi_ = phi.to_numpy() + f2v_ = f2v.to_numpy() + a, b, c = pos_[f2v_[:, 0]], pos_[f2v_[:, 1]], pos_[f2v_[:, 2]] + k = phi_ * (10 / E) + gb = (1 - k) * 0.5 + gui.triangles(a, b, c, color=ti.rgb_to_hex([k + gb, gb, gb])) + + +def main(): + init_mesh() + init_pos() + gravity[None] = [0, -1] + + gui = ti.GUI("FEM128") + print( + "[Hint] Use WSAD/arrow keys to control gravity. Use left/right mouse buttons to attract/repel. Press R to reset." + ) + while gui.running: + for e in gui.get_events(gui.PRESS): + if e.key == gui.ESCAPE: + gui.running = False + elif e.key == "r": + init_pos() + elif e.key in ("a", gui.LEFT): + gravity[None] = [-1, 0] + elif e.key in ("d", gui.RIGHT): + gravity[None] = [+1, 0] + elif e.key in ("s", gui.DOWN): + gravity[None] = [0, -1] + elif e.key in ("w", gui.UP): + gravity[None] = [0, +1] + mouse_pos = gui.get_cursor_pos() + attractor_pos[None] = mouse_pos + attractor_strength[None] = gui.is_pressed(gui.LMB) - gui.is_pressed(gui.RMB) + for i in range(50): + with ti.ad.Tape(loss=U): + update_U() + advance() + paint_phi(gui) + gui.circle(mouse_pos, radius=15, color=0x336699) + gui.circle(ball_pos, radius=ball_radius * 512, color=0x666666) + gui.circles(pos.to_numpy(), radius=2, color=0xFFAA33) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/fem99.py b/modules/createEVENT/TaichiEvent/fem99.py new file mode 100644 index 000000000..b7d9620f5 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/fem99.py @@ -0,0 +1,112 @@ +import taichi as ti + +ti.init(arch=ti.gpu) + +N = 64 +dt = 2.5e-5 +dx = 1 / N +rho = 4e1 +NF = 2 * N**2 # number of faces +NV = (N + 1) ** 2 # number of vertices +E, nu = 8e4, 0.2 # Young's modulus and Poisson's ratio +mu, lam = E / 2 / (1 + nu), E * nu / (1 + nu) / (1 - 2 * nu) # Lame parameters +ball_pos, ball_radius = ti.Vector([0.5, 0.0]), 0.32 +gravity = ti.Vector([0, -40]) +damping = 12.5 + +pos = ti.Vector.field(2, float, NV, needs_grad=True) +vel = ti.Vector.field(2, float, NV) +f2v = ti.Vector.field(3, int, NF) # ids of three vertices of each face +B = ti.Matrix.field(2, 2, float, NF) +F = ti.Matrix.field(2, 2, float, NF, needs_grad=True) +V = ti.field(float, NF) +phi = ti.field(float, NF) # potential energy of each face (Neo-Hookean) +U = ti.field(float, (), needs_grad=True) # total potential energy + + +@ti.kernel +def update_U(): + for i in range(NF): + ia, ib, ic = f2v[i] + a, b, c = pos[ia], pos[ib], pos[ic] + V[i] = abs((a - c).cross(b - c)) + D_i = ti.Matrix.cols([a - c, b - c]) + F[i] = D_i @ B[i] + for i in range(NF): + F_i = F[i] + log_J_i = ti.log(F_i.determinant()) + phi_i = mu / 2 * ((F_i.transpose() @ F_i).trace() - 2) + phi_i -= mu * log_J_i + phi_i += lam / 2 * log_J_i**2 + phi[i] = phi_i + U[None] += V[i] * phi_i + + +@ti.kernel +def advance(): + for i in range(NV): + acc = -pos.grad[i] / (rho * dx**2) + vel[i] += dt * (acc + gravity) + vel[i] *= ti.exp(-dt * damping) + for i in range(NV): + # ball boundary condition: + disp = pos[i] - ball_pos + disp2 = disp.norm_sqr() + if disp2 <= ball_radius**2: + NoV = vel[i].dot(disp) + if NoV < 0: + vel[i] -= NoV * disp / disp2 + # rect boundary condition: + cond = (pos[i] < 0) & (vel[i] < 0) | (pos[i] > 1) & (vel[i] > 0) + for j in ti.static(range(pos.n)): + if cond[j]: + vel[i][j] = 0 + pos[i] += dt * vel[i] + + +@ti.kernel +def init_pos(): + for i, j in ti.ndrange(N + 1, N + 1): + k = i * (N + 1) + j + pos[k] = ti.Vector([i, j]) / N * 0.25 + ti.Vector([0.45, 0.45]) + vel[k] = ti.Vector([0, 0]) + for i in range(NF): + ia, ib, ic = f2v[i] + a, b, c = pos[ia], pos[ib], pos[ic] + B_i_inv = ti.Matrix.cols([a - c, b - c]) + B[i] = B_i_inv.inverse() + + +@ti.kernel +def init_mesh(): + for i, j in ti.ndrange(N, N): + k = (i * N + j) * 2 + a = i * (N + 1) + j + b = a + 1 + c = a + N + 2 + d = a + N + 1 + f2v[k + 0] = [a, b, c] + f2v[k + 1] = [c, d, a] + + +def main(): + init_mesh() + init_pos() + gui = ti.GUI("FEM99") + while gui.running: + for e in gui.get_events(): + if e.key == gui.ESCAPE: + gui.running = False + elif e.key == "r": + init_pos() + for i in range(30): + with ti.ad.Tape(loss=U): + update_U() + advance() + gui.circles(pos.to_numpy(), radius=2, color=0xFFAA33) + gui.circle(ball_pos, radius=ball_radius * 512, color=0x666666) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/fractal.py b/modules/createEVENT/TaichiEvent/fractal.py new file mode 100644 index 000000000..dc63c7b22 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/fractal.py @@ -0,0 +1,37 @@ +import taichi as ti + +ti.init(arch=ti.gpu) + +n = 320 +pixels = ti.field(dtype=float, shape=(n * 2, n)) + + +@ti.func +def complex_sqr(z): + return ti.Vector([z[0] ** 2 - z[1] ** 2, z[1] * z[0] * 2]) + + +@ti.kernel +def paint(t: float): + for i, j in pixels: # Parallelized over all pixels + c = ti.Vector([-0.8, ti.cos(t) * 0.2]) + z = ti.Vector([i / n - 1, j / n - 0.5]) * 2 + iterations = 0 + while z.norm() < 20 and iterations < 50: + z = complex_sqr(z) + c + iterations += 1 + pixels[i, j] = 1 - iterations * 0.02 + + +def main(): + gui = ti.GUI("Julia Set", res=(n * 2, n)) + t = 0.0 + while not gui.get_event(ti.GUI.ESCAPE, ti.GUI.EXIT): + paint(t) + t += 0.03 + gui.set_image(pixels) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/game_of_life.py b/modules/createEVENT/TaichiEvent/game_of_life.py new file mode 100644 index 000000000..eeec3f677 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/game_of_life.py @@ -0,0 +1,105 @@ +# Game of Life written in 100 lines of Taichi +# In memory of John Horton Conway (1937 - 2020) + +import numpy as np + +import taichi as ti + +ti.init() + +n = 64 +cell_size = 8 +img_size = n * cell_size +alive = ti.field(int, shape=(n, n)) # alive = 1, dead = 0 +count = ti.field(int, shape=(n, n)) # count of neighbours + + +@ti.func +def get_alive(i, j): + return alive[i, j] if 0 <= i < n and 0 <= j < n else 0 + + +@ti.func +def get_count(i, j): + return ( + get_alive(i - 1, j) + + get_alive(i + 1, j) + + get_alive(i, j - 1) + + get_alive(i, j + 1) + + get_alive(i - 1, j - 1) + + get_alive(i + 1, j - 1) + + get_alive(i - 1, j + 1) + + get_alive(i + 1, j + 1) + ) + + +# See https://www.conwaylife.com/wiki/Cellular_automaton#Rules for more rules +B, S = [3], [2, 3] +# B, S = [2], [0] + + +@ti.func +def calc_rule(a, c): + if a == 0: + for t in ti.static(B): + if c == t: + a = 1 + elif a == 1: + a = 0 + for t in ti.static(S): + if c == t: + a = 1 + return a + + +@ti.kernel +def run(): + for i, j in alive: + count[i, j] = get_count(i, j) + + for i, j in alive: + alive[i, j] = calc_rule(alive[i, j], count[i, j]) + + +@ti.kernel +def init(): + for i, j in alive: + if ti.random() > 0.8: + alive[i, j] = 1 + else: + alive[i, j] = 0 + + +def main(): + gui = ti.GUI("Game of Life", (img_size, img_size)) + gui.fps_limit = 15 + + print("[Hint] Press `r` to reset") + print("[Hint] Press SPACE to pause") + print("[Hint] Click LMB, RMB and drag to add alive / dead cells") + + init() + paused = False + while gui.running: + for e in gui.get_events(gui.PRESS, gui.MOTION): + if e.key == gui.ESCAPE: + gui.running = False + elif e.key == gui.SPACE: + paused = not paused + elif e.key == "r": + alive.fill(0) + + if gui.is_pressed(gui.LMB, gui.RMB): + mx, my = gui.get_cursor_pos() + alive[int(mx * n), int(my * n)] = gui.is_pressed(gui.LMB) + paused = True + + if not paused: + run() + + gui.set_image(ti.tools.imresize(alive, img_size).astype(np.uint8) * 255) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/implicit_fem.py b/modules/createEVENT/TaichiEvent/implicit_fem.py new file mode 100644 index 000000000..5ddfc766e --- /dev/null +++ b/modules/createEVENT/TaichiEvent/implicit_fem.py @@ -0,0 +1,436 @@ +import argparse + +import numpy as np +from taichi._lib import core as _ti_core + +import taichi as ti + +parser = argparse.ArgumentParser() +parser.add_argument("--exp", choices=["implicit", "explicit"], default="implicit") +parser.add_argument("--dim", type=int, default=3) +parser.add_argument("--gui", choices=["auto", "ggui", "cpu"], default="auto") +parser.add_argument( + "-s", + "--use_sparse", + action="store_true", + help="Use sparse matrix and sparse solver", +) +parser.add_argument("place_holder", nargs="*") +args = parser.parse_args() + +ti.init(arch=ti.cuda) + +if args.gui == "auto": + if _ti_core.GGUI_AVAILABLE and ti.lang.impl.current_cfg().arch == ti.cuda: + args.gui = "ggui" + else: + args.gui = "cpu" + +E, nu = 5e4, 0.0 +mu, la = E / (2 * (1 + nu)), E * nu / ((1 + nu) * (1 - 2 * nu)) # lambda = 0 +density = 1000.0 +dt = 2e-4 + +if args.exp == "implicit": + dt = 1e-2 + +use_sparse = args.use_sparse + +n_cube = np.array([5] * 3) +n_verts = np.product(n_cube) +n_cells = 5 * np.product(n_cube - 1) +dx = 1 / (n_cube.max() - 1) + +F_vertices = ti.Vector.field(4, dtype=ti.i32, shape=n_cells) + +F_x = ti.Vector.field(args.dim, dtype=ti.f32, shape=n_verts) +F_ox = ti.Vector.field(args.dim, dtype=ti.f32, shape=n_verts) +F_v = ti.Vector.field(args.dim, dtype=ti.f32, shape=n_verts) +F_f = ti.Vector.field(args.dim, dtype=ti.f32, shape=n_verts) +F_mul_ans = ti.Vector.field(args.dim, dtype=ti.f32, shape=n_verts) +F_m = ti.field(dtype=ti.f32, shape=n_verts) + +n_cells = (n_cube - 1).prod() * 5 +F_B = ti.Matrix.field(args.dim, args.dim, dtype=ti.f32, shape=n_cells) +F_W = ti.field(dtype=ti.f32, shape=n_cells) + + +@ti.func +def i2p(I): + return (I.x * n_cube[1] + I.y) * n_cube[2] + I.z + + +@ti.func +def set_element(e, I, verts): + for i in ti.static(range(args.dim + 1)): + F_vertices[e][i] = i2p(I + (([verts[i] >> k for k in range(3)] ^ I) & 1)) + + +@ti.kernel +def get_vertices(): + """ + This kernel partitions the cube into tetrahedrons. + Each unit cube is divided into 5 tetrahedrons. + """ + for I in ti.grouped(ti.ndrange(*(n_cube - 1))): + e = ((I.x * (n_cube[1] - 1) + I.y) * (n_cube[2] - 1) + I.z) * 5 + for i, j in ti.static(enumerate([0, 3, 5, 6])): + set_element(e + i, I, (j, j ^ 1, j ^ 2, j ^ 4)) + set_element(e + 4, I, (1, 2, 4, 7)) + for I in ti.grouped(ti.ndrange(*(n_cube))): + F_ox[i2p(I)] = I * dx + + +@ti.func +def Ds(verts): + return ti.Matrix.cols([F_x[verts[i]] - F_x[verts[3]] for i in range(3)]) + + +@ti.func +def ssvd(F): + U, sig, V = ti.svd(F) + if U.determinant() < 0: + for i in ti.static(range(3)): + U[i, 2] *= -1 + sig[2, 2] = -sig[2, 2] + if V.determinant() < 0: + for i in ti.static(range(3)): + V[i, 2] *= -1 + sig[2, 2] = -sig[2, 2] + return U, sig, V + + +@ti.func +def get_force_func(c, verts): + F = Ds(verts) @ F_B[c] + P = ti.Matrix.zero(ti.f32, 3, 3) + U, sig, V = ssvd(F) + P = 2 * mu * (F - U @ V.transpose()) + H = -F_W[c] * P @ F_B[c].transpose() + for i in ti.static(range(3)): + force = ti.Vector([H[j, i] for j in range(3)]) + F_f[verts[i]] += force + F_f[verts[3]] -= force + + +@ti.kernel +def get_force(): + for c in F_vertices: + get_force_func(c, F_vertices[c]) + for u in F_f: + F_f[u].y -= 9.8 * F_m[u] + + +@ti.kernel +def matmul_cell(ret: ti.template(), vel: ti.template()): + for i in ret: + ret[i] = vel[i] * F_m[i] + for c in F_vertices: + verts = F_vertices[c] + W_c = F_W[c] + B_c = F_B[c] + for u in range(4): + for d in range(3): + dD = ti.Matrix.zero(ti.f32, 3, 3) + if u == 3: + for j in range(3): + dD[d, j] = -1 + else: + dD[d, u] = 1 + dF = dD @ B_c + dP = 2.0 * mu * dF + dH = -W_c * dP @ B_c.transpose() + for i in range(3): + for j in range(3): + tmp = vel[verts[i]][j] - vel[verts[3]][j] + ret[verts[u]][d] += -(dt**2) * dH[j, i] * tmp + + +@ti.kernel +def add(input: ti.template(), a: ti.template(), k: ti.f32, b: ti.template()): + for i in input: + input[i] = a[i] + k * b[i] + + +@ti.kernel +def dot(a: ti.template(), b: ti.template()) -> ti.f32: + value = 0.0 + for i in a: + value += a[i].dot(b[i]) + return value + + +F_b = ti.Vector.field(3, dtype=ti.f32, shape=n_verts) +F_r0 = ti.Vector.field(3, dtype=ti.f32, shape=n_verts) +F_p0 = ti.Vector.field(3, dtype=ti.f32, shape=n_verts) + +# ndarray version of F_b +F_b_ndarr = ti.ndarray(dtype=ti.f32, shape=3 * n_verts) +# stiffness matrix +A_builder = ti.linalg.SparseMatrixBuilder(3 * n_verts, 3 * n_verts, max_num_triplets=50000) +solver = ti.linalg.SparseSolver(ti.f32, "LLT") + + +@ti.kernel +def get_b(): + for i in F_b: + F_b[i] = F_m[i] * F_v[i] + dt * F_f[i] + + +def cg(): + def mul(x): + matmul_cell(F_mul_ans, x) + return F_mul_ans + + get_force() + get_b() + mul(F_v) + add(F_r0, F_b, -1, mul(F_v)) + + d = F_p0 + d.copy_from(F_r0) + r_2 = dot(F_r0, F_r0) + n_iter = 50 + epsilon = 1e-6 + r_2_init = r_2 + r_2_new = r_2 + for _ in range(n_iter): + q = mul(d) + alpha = r_2_new / dot(d, q) + add(F_v, F_v, alpha, d) + add(F_r0, F_r0, -alpha, q) + r_2 = r_2_new + r_2_new = dot(F_r0, F_r0) + if r_2_new <= r_2_init * epsilon**2: + break + beta = r_2_new / r_2 + add(d, F_r0, beta, d) + F_f.fill(0) + add(F_x, F_x, dt, F_v) + + +@ti.kernel +def compute_A(A: ti.types.sparse_matrix_builder()): + # A = M - dt * dt * K + for i in range(n_verts): + for j in range(3): + A[3 * i + j, 3 * i + j] += F_m[i] + for c in F_vertices: + verts = F_vertices[c] + W_c = F_W[c] + B_c = F_B[c] + for u in range(4): + for d in range(3): + dD = ti.Matrix.zero(ti.f32, 3, 3) + if u == 3: + for j in range(3): + dD[d, j] = -1 + else: + dD[d, u] = 1 + dF = dD @ B_c + dP = 2.0 * mu * dF + dH = -W_c * dP @ B_c.transpose() + for i in range(3): + for j in range(3): + A[3 * verts[u] + d, 3 * verts[i] + j] += -(dt**2) * dH[j, i] + for i in range(3): + A[3 * verts[u] + d, 3 * verts[3] + i] += -(dt**2) * (-dH[i, 0] - dH[i, 1] - dH[i, 2]) + + +@ti.kernel +def flatten(dest: ti.types.ndarray(), src: ti.template()): + for i in range(n_verts): + for j in range(3): + dest[3 * i + j] = src[i][j] + + +@ti.kernel +def aggregate(dest: ti.template(), src: ti.types.ndarray()): + for i in range(n_verts): + for j in range(3): + dest[i][j] = src[3 * i + j] + + +def direct(): + get_force() + get_b() + flatten(F_b_ndarr, F_b) + v = solver.solve(F_b_ndarr) + aggregate(F_v, v) + F_f.fill(0) + add(F_x, F_x, dt, F_v) + + +@ti.kernel +def advect(): + for p in F_x: + F_v[p] += dt * (F_f[p] / F_m[p]) + F_x[p] += dt * F_v[p] + F_f[p] = ti.Vector([0, 0, 0]) + + +@ti.kernel +def init(): + for u in F_x: + F_x[u] = F_ox[u] + F_v[u] = [0.0] * 3 + F_f[u] = [0.0] * 3 + F_m[u] = 0.0 + for c in F_vertices: + F = Ds(F_vertices[c]) + F_B[c] = F.inverse() + F_W[c] = ti.abs(F.determinant()) / 6 + for i in range(4): + F_m[F_vertices[c][i]] += F_W[c] / 4 * density + for u in F_x: + F_x[u].y += 1.0 + + +def init_A(): + compute_A(A_builder) + A = A_builder.build() + solver.analyze_pattern(A) + solver.factorize(A) + + +@ti.kernel +def floor_bound(): + for u in F_x: + if F_x[u].y < 0: + F_x[u].y = 0 + if F_v[u].y < 0: + F_v[u].y = 0 + + +@ti.func +def check(u): + value = 0 + rest = u + for i in ti.static(range(3)): + k = rest % n_cube[2 - i] + rest = rest // n_cube[2 - i] + if k == 0: + value |= 1 << (i * 2) + if k == n_cube[2 - i] - 1: + value |= 1 << (i * 2 + 1) + return value + + +def gen_indices(): + su = 0 + for i in range(3): + su += (n_cube[i] - 1) * (n_cube[(i + 1) % 3] - 1) + return ti.field(ti.i32, shape=2 * su * 2 * 3) + + +indices = gen_indices() + + +@ti.kernel +def get_indices(): + # calculate all the meshes on surface + cnt = 0 + for c in F_vertices: + if c % 5 != 4: + for i in ti.static([0, 2, 3]): + verts = [F_vertices[c][(i + j) % 4] for j in range(3)] + sum_ = check(verts[0]) & check(verts[1]) & check(verts[2]) + if sum_: + m = ti.atomic_add(cnt, 1) + det = ti.Matrix.rows([F_x[verts[i]] - [0.5, 1.5, 0.5] for i in range(3)]).determinant() + if det < 0: + tmp = verts[1] + verts[1] = verts[2] + verts[2] = tmp + indices[m * 3] = verts[0] + indices[m * 3 + 1] = verts[1] + indices[m * 3 + 2] = verts[2] + + +def substep(): + if args.exp == "explicit": + for i in range(10): + get_force() + advect() + else: + if use_sparse: + direct() + else: + cg() + floor_bound() + + +def main(): + get_vertices() + init() + get_indices() + + if use_sparse: + init_A() + + if args.gui == "ggui": + res = (800, 600) + window = ti.ui.Window("Implicit FEM", res, vsync=True) + + frame_id = 0 + canvas = window.get_canvas() + scene = window.get_scene() + camera = ti.ui.Camera() + camera.position(2.0, 2.0, 3.95) + camera.lookat(0.5, 0.5, 0.5) + camera.fov(55) + + def render(): + camera.track_user_inputs(window, movement_speed=0.03, hold_key=ti.ui.RMB) + scene.set_camera(camera) + + scene.ambient_light((0.1,) * 3) + + scene.point_light(pos=(0.5, 10.0, 0.5), color=(0.5, 0.5, 0.5)) + scene.point_light(pos=(10.0, 10.0, 10.0), color=(0.5, 0.5, 0.5)) + + scene.mesh(F_x, indices, color=(0.73, 0.33, 0.23)) + + canvas.scene(scene) + + while window.running: + frame_id += 1 + frame_id = frame_id % 256 + substep() + if window.is_pressed("r"): + init() + if window.is_pressed(ti.GUI.ESCAPE): + break + + render() + + window.show() + + else: + + def T(a): + phi, theta = np.radians(28), np.radians(32) + + a = a - 0.2 + x, y, z = a[:, 0], a[:, 1], a[:, 2] + c, s = np.cos(phi), np.sin(phi) + C, S = np.cos(theta), np.sin(theta) + x, z = x * c + z * s, z * c - x * s + u, v = x, y * C + z * S + return np.array([u, v]).swapaxes(0, 1) + 0.5 + + gui = ti.GUI("Implicit FEM") + while gui.running: + substep() + if gui.get_event(ti.GUI.PRESS): + if gui.event.key in [ti.GUI.ESCAPE, ti.GUI.EXIT]: + break + if gui.is_pressed("r"): + init() + gui.clear(0x000000) + gui.circles(T(F_x.to_numpy() / 3), radius=1.5, color=0xBA543A) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/implicit_mass_spring.py b/modules/createEVENT/TaichiEvent/implicit_mass_spring.py new file mode 100644 index 000000000..9f7213228 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/implicit_mass_spring.py @@ -0,0 +1,293 @@ +# https://www.cs.cmu.edu/~baraff/papers/sig98.pdf +import argparse + +import numpy as np + +import taichi as ti + + +@ti.data_oriented +class Cloth: + def __init__(self, N): + self.N = N + self.NF = 2 * N**2 # number of faces + self.NV = (N + 1) ** 2 # number of vertices + self.NE = 2 * N * (N + 1) + 2 * N * N # numbser of edges + self.pos = ti.Vector.field(2, ti.f32, self.NV) + self.initPos = ti.Vector.field(2, ti.f32, self.NV) + self.vel = ti.Vector.field(2, ti.f32, self.NV) + self.force = ti.Vector.field(2, ti.f32, self.NV) + self.mass = ti.field(ti.f32, self.NV) + self.vel_1D = ti.ndarray(ti.f32, 2 * self.NV) + self.force_1D = ti.ndarray(ti.f32, 2 * self.NV) + self.b = ti.ndarray(ti.f32, 2 * self.NV) + + self.spring = ti.Vector.field(2, ti.i32, self.NE) + self.indices = ti.field(ti.i32, 2 * self.NE) + self.Jx = ti.Matrix.field(2, 2, ti.f32, self.NE) # Jacobian with respect to position + self.Jv = ti.Matrix.field(2, 2, ti.f32, self.NE) # Jacobian with respect to velocity + self.rest_len = ti.field(ti.f32, self.NE) + self.ks = 1000.0 # spring stiffness + self.kd = 0.5 # damping constant + self.kf = 1.0e5 # fix point stiffness + + self.gravity = ti.Vector([0.0, -2.0]) + self.init_pos() + self.init_edges() + self.MassBuilder = ti.linalg.SparseMatrixBuilder(2 * self.NV, 2 * self.NV, max_num_triplets=10000) + self.DBuilder = ti.linalg.SparseMatrixBuilder(2 * self.NV, 2 * self.NV, max_num_triplets=10000) + self.KBuilder = ti.linalg.SparseMatrixBuilder(2 * self.NV, 2 * self.NV, max_num_triplets=10000) + self.init_mass_sp(self.MassBuilder) + self.M = self.MassBuilder.build() + self.fix_vertex = [self.N, self.NV - 1] + self.Jf = ti.Matrix.field(2, 2, ti.f32, len(self.fix_vertex)) # fix constraint hessian + + @ti.kernel + def init_pos(self): + for i, j in ti.ndrange(self.N + 1, self.N + 1): + k = i * (self.N + 1) + j + self.pos[k] = ti.Vector([i, j]) / self.N * 0.5 + ti.Vector([0.25, 0.25]) + self.initPos[k] = self.pos[k] + self.vel[k] = ti.Vector([0, 0]) + self.mass[k] = 1.0 + + @ti.kernel + def init_edges(self): + pos, spring, N, rest_len = ti.static(self.pos, self.spring, self.N, self.rest_len) + for i, j in ti.ndrange(N + 1, N): + idx, idx1 = i * N + j, i * (N + 1) + j + spring[idx] = ti.Vector([idx1, idx1 + 1]) + rest_len[idx] = (pos[idx1] - pos[idx1 + 1]).norm() + start = N * (N + 1) + for i, j in ti.ndrange(N, N + 1): + idx, idx1, idx2 = ( + start + i + j * N, + i * (N + 1) + j, + i * (N + 1) + j + N + 1, + ) + spring[idx] = ti.Vector([idx1, idx2]) + rest_len[idx] = (pos[idx1] - pos[idx2]).norm() + start = 2 * N * (N + 1) + for i, j in ti.ndrange(N, N): + idx, idx1, idx2 = ( + start + i * N + j, + i * (N + 1) + j, + (i + 1) * (N + 1) + j + 1, + ) + spring[idx] = ti.Vector([idx1, idx2]) + rest_len[idx] = (pos[idx1] - pos[idx2]).norm() + start = 2 * N * (N + 1) + N * N + for i, j in ti.ndrange(N, N): + idx, idx1, idx2 = ( + start + i * N + j, + i * (N + 1) + j + 1, + (i + 1) * (N + 1) + j, + ) + spring[idx] = ti.Vector([idx1, idx2]) + rest_len[idx] = (pos[idx1] - pos[idx2]).norm() + + @ti.kernel + def init_mass_sp(self, M: ti.types.sparse_matrix_builder()): + for i in range(self.NV): + mass = self.mass[i] + M[2 * i + 0, 2 * i + 0] += mass + M[2 * i + 1, 2 * i + 1] += mass + + @ti.func + def clear_force(self): + for i in self.force: + self.force[i] = ti.Vector([0.0, 0.0]) + + @ti.kernel + def compute_force(self): + self.clear_force() + for i in self.force: + self.force[i] += self.gravity * self.mass[i] + + for i in self.spring: + idx1, idx2 = self.spring[i][0], self.spring[i][1] + pos1, pos2 = self.pos[idx1], self.pos[idx2] + dis = pos2 - pos1 + force = self.ks * (dis.norm() - self.rest_len[i]) * dis.normalized() + self.force[idx1] += force + self.force[idx2] -= force + # fix constraint gradient + self.force[self.N] += self.kf * (self.initPos[self.N] - self.pos[self.N]) + self.force[self.NV - 1] += self.kf * (self.initPos[self.NV - 1] - self.pos[self.NV - 1]) + + @ti.kernel + def compute_Jacobians(self): + for i in self.spring: + idx1, idx2 = self.spring[i][0], self.spring[i][1] + pos1, pos2 = self.pos[idx1], self.pos[idx2] + dx = pos1 - pos2 + I = ti.Matrix([[1.0, 0.0], [0.0, 1.0]]) + dxtdx = ti.Matrix([[dx[0] * dx[0], dx[0] * dx[1]], [dx[1] * dx[0], dx[1] * dx[1]]]) + l = dx.norm() + if l != 0.0: + l = 1.0 / l + self.Jx[i] = (I - self.rest_len[i] * l * (I - dxtdx * l**2)) * self.ks + self.Jv[i] = self.kd * I + + # fix point constraint hessian + self.Jf[0] = ti.Matrix([[-self.kf, 0], [0, -self.kf]]) + self.Jf[1] = ti.Matrix([[-self.kf, 0], [0, -self.kf]]) + + @ti.kernel + def assemble_K(self, K: ti.types.sparse_matrix_builder()): + for i in self.spring: + idx1, idx2 = self.spring[i][0], self.spring[i][1] + for m, n in ti.static(ti.ndrange(2, 2)): + K[2 * idx1 + m, 2 * idx1 + n] -= self.Jx[i][m, n] + K[2 * idx1 + m, 2 * idx2 + n] += self.Jx[i][m, n] + K[2 * idx2 + m, 2 * idx1 + n] += self.Jx[i][m, n] + K[2 * idx2 + m, 2 * idx2 + n] -= self.Jx[i][m, n] + for m, n in ti.static(ti.ndrange(2, 2)): + K[2 * self.N + m, 2 * self.N + n] += self.Jf[0][m, n] + K[2 * (self.NV - 1) + m, 2 * (self.NV - 1) + n] += self.Jf[1][m, n] + + @ti.kernel + def assemble_D(self, D: ti.types.sparse_matrix_builder()): + for i in self.spring: + idx1, idx2 = self.spring[i][0], self.spring[i][1] + for m, n in ti.static(ti.ndrange(2, 2)): + D[2 * idx1 + m, 2 * idx1 + n] -= self.Jv[i][m, n] + D[2 * idx1 + m, 2 * idx2 + n] += self.Jv[i][m, n] + D[2 * idx2 + m, 2 * idx1 + n] += self.Jv[i][m, n] + D[2 * idx2 + m, 2 * idx2 + n] -= self.Jv[i][m, n] + + @ti.kernel + def updatePosVel(self, h: ti.f32, dv: ti.types.ndarray()): + for i in self.pos: + self.vel[i] += ti.Vector([dv[2 * i], dv[2 * i + 1]]) + self.pos[i] += h * self.vel[i] + + @ti.kernel + def copy_to(self, des: ti.types.ndarray(), source: ti.template()): + for i in range(self.NV): + des[2 * i] = source[i][0] + des[2 * i + 1] = source[i][1] + + @ti.kernel + def compute_b( + self, + b: ti.types.ndarray(), + f: ti.types.ndarray(), + Kv: ti.types.ndarray(), + h: ti.f32, + ): + for i in range(2 * self.NV): + b[i] = (f[i] + Kv[i] * h) * h + + def update(self, h): + self.compute_force() + + self.compute_Jacobians() + # Assemble global system + + self.assemble_D(self.DBuilder) + D = self.DBuilder.build() + + self.assemble_K(self.KBuilder) + K = self.KBuilder.build() + + A = self.M - h * D - h**2 * K + + self.copy_to(self.vel_1D, self.vel) + self.copy_to(self.force_1D, self.force) + + # b = (force + h * K @ vel) * h + Kv = K @ self.vel_1D + self.compute_b(self.b, self.force_1D, Kv, h) + + # Sparse solver + solver = ti.linalg.SparseSolver(solver_type="LDLT") + solver.analyze_pattern(A) + solver.factorize(A) + # Solve the linear system + dv = solver.solve(self.b) + self.updatePosVel(h, dv) + + def display(self, gui, radius=5, color=0xFFFFFF): + lines = self.spring.to_numpy() + pos = self.pos.to_numpy() + edgeBegin = np.zeros(shape=(lines.shape[0], 2)) + edgeEnd = np.zeros(shape=(lines.shape[0], 2)) + for i in range(lines.shape[0]): + idx1, idx2 = lines[i][0], lines[i][1] + edgeBegin[i] = pos[idx1] + edgeEnd[i] = pos[idx2] + gui.lines(edgeBegin, edgeEnd, radius=2, color=0x0000FF) + gui.circles(self.pos.to_numpy(), radius, color) + + @ti.kernel + def spring2indices(self): + for i in self.spring: + self.indices[2 * i + 0] = self.spring[i][0] + self.indices[2 * i + 1] = self.spring[i][1] + + def displayGGUI(self, canvas, radius=0.01, color=(1.0, 1.0, 1.0)): + self.spring2indices() + canvas.lines(self.pos, width=0.005, indices=self.indices, color=(0.0, 0.0, 1.0)) + canvas.circles(self.pos, radius, color) + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("-g", "--use-ggui", action="store_true", help="Display with GGUI") + parser.add_argument( + "-a", + "--arch", + required=False, + default="cpu", + dest="arch", + type=str, + help="The arch (backend) to run this example on", + ) + args, unknowns = parser.parse_known_args() + arch = args.arch + if arch in ["x64", "cpu", "arm64"]: + ti.init(arch=ti.cpu) + elif arch in ["cuda", "gpu"]: + ti.init(arch=ti.cuda) + else: + raise ValueError("Only CPU and CUDA backends are supported for now.") + + h = 0.01 + pause = False + cloth = Cloth(N=5) + + use_ggui = args.use_ggui + if not use_ggui: + gui = ti.GUI("Implicit Mass Spring System", res=(500, 500)) + while gui.running: + for e in gui.get_events(): + if e.key == gui.ESCAPE: + gui.running = False + elif e.key == gui.SPACE: + pause = not pause + + if not pause: + cloth.update(h) + + cloth.display(gui) + gui.show() + else: + window = ti.ui.Window("Implicit Mass Spring System", res=(500, 500)) + while window.running: + if window.get_event(ti.ui.PRESS): + if window.event.key == ti.ui.ESCAPE: + break + if window.is_pressed(ti.ui.SPACE): + pause = not pause + + if not pause: + cloth.update(h) + + canvas = window.get_canvas() + cloth.displayGGUI(canvas) + window.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/initial_value_problem.py b/modules/createEVENT/TaichiEvent/initial_value_problem.py new file mode 100644 index 000000000..2361b6abb --- /dev/null +++ b/modules/createEVENT/TaichiEvent/initial_value_problem.py @@ -0,0 +1,50 @@ +import time + +import numpy as np + +import taichi as ti + + +def init(): + a = [] + for i in np.linspace(0, 1, n, False): + for j in np.linspace(0, 1, n, False): + a.append([i, j]) + return np.array(a).astype(np.float32) + + +ti.init(arch=ti.gpu) +n = 50 +dirs = ti.field(dtype=float, shape=(n * n, 2)) +locations_np = init() + +locations = ti.field(dtype=float, shape=(n * n, 2)) +locations.from_numpy(locations_np) + + +@ti.kernel +def paint(t: float): + (o, p) = locations_np.shape + for i in range(0, o): # Parallelized over all pixels + x = locations[i, 0] + y = locations[i, 1] + dirs[i, 0] = ti.sin((t * x - y)) + dirs[i, 1] = ti.cos(t * y - x) + l = (dirs[i, 0] ** 2 + dirs[i, 1] ** 2) ** 0.5 + dirs[i, 0] /= l * 40 + dirs[i, 1] /= l * 40 + + +def main(): + gui = ti.GUI("Vector Field", res=(500, 500)) + + beginning = time.time_ns() + for k in range(1000000): + paint((time.time_ns() - beginning) * 0.00000001) + dirs_np = dirs.to_numpy() + gui.arrows(locations_np, dirs_np, radius=1) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/karman_vortex_street.py b/modules/createEVENT/TaichiEvent/karman_vortex_street.py new file mode 100644 index 000000000..8341ccbfd --- /dev/null +++ b/modules/createEVENT/TaichiEvent/karman_vortex_street.py @@ -0,0 +1,177 @@ +# Fluid solver based on lattice boltzmann method using taichi language +# Author : Wang (hietwll@gmail.com) +# Original code at https://github.com/hietwll/LBM_Taichi + +import matplotlib +import numpy as np +from matplotlib import cm + +import taichi as ti +import taichi.math as tm + +ti.init(arch=ti.gpu) + +vec9 = ti.types.vector(9, float) +imat9x2 = ti.types.matrix(9, 2, int) + + +@ti.data_oriented +class lbm_solver: + def __init__( + self, + nx, # domain size + ny, + niu, # viscosity of fluid + bc_type, # [left,top,right,bottom] boundary conditions: 0 -> Dirichlet ; 1 -> Neumann + bc_value, # if bc_type = 0, we need to specify the velocity in bc_value + cy=0, # whether to place a cylindrical obstacle + cy_para=[0.0, 0.0, 0.0], # location and radius of the cylinder + ): + self.nx = nx # by convention, dx = dy = dt = 1.0 (lattice units) + self.ny = ny + self.niu = niu + self.tau = 3.0 * niu + 0.5 + self.inv_tau = 1.0 / self.tau + self.rho = ti.field(float, shape=(nx, ny)) + self.vel = ti.Vector.field(2, float, shape=(nx, ny)) + self.mask = ti.field(float, shape=(nx, ny)) + self.f_old = ti.Vector.field(9, float, shape=(nx, ny)) + self.f_new = ti.Vector.field(9, float, shape=(nx, ny)) + self.w = vec9(4, 1, 1, 1, 1, 1 / 4, 1 / 4, 1 / 4, 1 / 4) / 9.0 + self.e = imat9x2([0, 0], [1, 0], [0, 1], [-1, 0], [0, -1], [1, 1], [-1, 1], [-1, -1], [1, -1]) + self.bc_type = ti.field(int, 4) + self.bc_type.from_numpy(np.array(bc_type, dtype=np.int32)) + self.bc_value = ti.Vector.field(2, float, shape=4) + self.bc_value.from_numpy(np.array(bc_value, dtype=np.float32)) + self.cy = cy + self.cy_para = tm.vec3(cy_para) + + @ti.func # compute equilibrium distribution function + def f_eq(self, i, j): + eu = self.e @ self.vel[i, j] + uv = tm.dot(self.vel[i, j], self.vel[i, j]) + return self.w * self.rho[i, j] * (1 + 3 * eu + 4.5 * eu * eu - 1.5 * uv) + + @ti.kernel + def init(self): + self.vel.fill(0) + self.rho.fill(1) + self.mask.fill(0) + for i, j in self.rho: + self.f_old[i, j] = self.f_new[i, j] = self.f_eq(i, j) + if self.cy == 1: + if (i - self.cy_para[0]) ** 2 + (j - self.cy_para[1]) ** 2 <= self.cy_para[2] ** 2: + self.mask[i, j] = 1.0 + + @ti.kernel + def collide_and_stream(self): # lbm core equation + for i, j in ti.ndrange((1, self.nx - 1), (1, self.ny - 1)): + for k in ti.static(range(9)): + ip = i - self.e[k, 0] + jp = j - self.e[k, 1] + feq = self.f_eq(ip, jp) + self.f_new[i, j][k] = (1 - self.inv_tau) * self.f_old[ip, jp][k] + feq[k] * self.inv_tau + + @ti.kernel + def update_macro_var(self): # compute rho u v + for i, j in ti.ndrange((1, self.nx - 1), (1, self.ny - 1)): + self.rho[i, j] = 0 + self.vel[i, j] = 0, 0 + for k in ti.static(range(9)): + self.f_old[i, j][k] = self.f_new[i, j][k] + self.rho[i, j] += self.f_new[i, j][k] + self.vel[i, j] += tm.vec2(self.e[k, 0], self.e[k, 1]) * self.f_new[i, j][k] + + self.vel[i, j] /= self.rho[i, j] + + @ti.kernel + def apply_bc(self): # impose boundary conditions + # left and right + for j in range(1, self.ny - 1): + # left: dr = 0; ibc = 0; jbc = j; inb = 1; jnb = j + self.apply_bc_core(1, 0, 0, j, 1, j) + + # right: dr = 2; ibc = nx-1; jbc = j; inb = nx-2; jnb = j + self.apply_bc_core(1, 2, self.nx - 1, j, self.nx - 2, j) + + # top and bottom + for i in range(self.nx): + # top: dr = 1; ibc = i; jbc = ny-1; inb = i; jnb = ny-2 + self.apply_bc_core(1, 1, i, self.ny - 1, i, self.ny - 2) + + # bottom: dr = 3; ibc = i; jbc = 0; inb = i; jnb = 1 + self.apply_bc_core(1, 3, i, 0, i, 1) + + # cylindrical obstacle + # Note: for cuda backend, putting 'if statement' inside loops can be much faster! + for i, j in ti.ndrange(self.nx, self.ny): + if self.cy == 1 and self.mask[i, j] == 1: + self.vel[i, j] = 0, 0 # velocity is zero at solid boundary + inb = 0 + jnb = 0 + if i >= self.cy_para[0]: + inb = i + 1 + else: + inb = i - 1 + if j >= self.cy_para[1]: + jnb = j + 1 + else: + jnb = j - 1 + self.apply_bc_core(0, 0, i, j, inb, jnb) + + @ti.func + def apply_bc_core(self, outer, dr, ibc, jbc, inb, jnb): + if outer == 1: # handle outer boundary + if self.bc_type[dr] == 0: + self.vel[ibc, jbc] = self.bc_value[dr] + + elif self.bc_type[dr] == 1: + self.vel[ibc, jbc] = self.vel[inb, jnb] + + self.rho[ibc, jbc] = self.rho[inb, jnb] + self.f_old[ibc, jbc] = self.f_eq(ibc, jbc) - self.f_eq(inb, jnb) + self.f_old[inb, jnb] + + def solve(self): + gui = ti.GUI("Karman Vortex Street", (self.nx, 2 * self.ny)) + self.init() + while not gui.get_event(ti.GUI.ESCAPE, ti.GUI.EXIT): + for _ in range(10): + self.collide_and_stream() + self.update_macro_var() + self.apply_bc() + + ## code fragment displaying vorticity is contributed by woclass + vel = self.vel.to_numpy() + ugrad = np.gradient(vel[:, :, 0]) + vgrad = np.gradient(vel[:, :, 1]) + vorticity = ugrad[1] - vgrad[0] + vel_mag = (vel[:, :, 0] ** 2.0 + vel[:, :, 1] ** 2.0) ** 0.5 + ## color map + colors = [ + (1, 1, 0), + (0.953, 0.490, 0.016), + (0, 0, 0), + (0.176, 0.976, 0.529), + (0, 1, 1), + ] + my_cmap = matplotlib.colors.LinearSegmentedColormap.from_list("my_cmap", colors) + vor_img = cm.ScalarMappable(norm=matplotlib.colors.Normalize(vmin=-0.02, vmax=0.02), cmap=my_cmap).to_rgba( + vorticity + ) + vel_img = cm.plasma(vel_mag / 0.15) + img = np.concatenate((vor_img, vel_img), axis=1) + gui.set_image(img) + gui.show() + + +if __name__ == "__main__": + lbm = lbm_solver( + 801, + 201, + 0.01, + [0, 0, 1, 0], + [[0.1, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], + 1, + [160.0, 100.0, 20.0], + ) + lbm.solve() diff --git a/modules/createEVENT/TaichiEvent/laplace_equation.py b/modules/createEVENT/TaichiEvent/laplace_equation.py new file mode 100644 index 000000000..16fd2f9a9 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/laplace_equation.py @@ -0,0 +1,287 @@ +import taichi as ti +import taichi.math as tm + +ti.init() +vec2 = tm.vec2 + +# Author: bismarckkk +# 2D-Simulation of fundamental solutions of Laplace equation +# Refer to "Aerodynamics" Chapter 3 (ISBN:9787030582287) + +# Pressing left mouse button with 's/v/d' key will create new 'source/vortex/dipole' on cursor position, +# and pressing right mouse button will create negative one. +# Pressing left and right mouse buttons will delete element around cursor position. +# Pressing left/right mouse button without other keys will increase/decrease the intensity of element around cursor. + +screen = ( + 30, + 20, +) # density of point generation positions on the boundary, also decide window size +arrowField = [int(it / 2) for it in screen] # number of arrow in window +meshSpace = 20 # screen * meshSpace = windowSize +maxElements = 20 # the max number of each element(source/sink, vortex, dipole) +refillFrame = 20 # frame interval for each refill points +refillVelThresh = ( + 0.2 # points will be refilled when the absolute value of the velocity on boundary is greater than this value +) +V = vec2(1.0, 0) # the velocity of uniform stream +dt = 0.002 +fadeMax = 10 # frames of fade in/out +fade = fadeMax # control arrows' fade in and fade out +maxPoints = screen[0] * screen[1] * 2 +gui_ = ti.GUI("demo", tuple(it * meshSpace for it in screen)) +guiHeight = meshSpace * screen[1] + +points = ti.Vector.field(2, float, maxPoints) +sources = ti.Struct.field({"pos": vec2, "q": ti.f32}, shape=maxElements) +vortexes = ti.Struct.field({"pos": vec2, "q": ti.f32}, shape=maxElements) +dipoles = ti.Struct.field({"pos": vec2, "m": ti.f32}, shape=maxElements) +arrows = ti.Struct.field({"base": vec2, "dir": vec2, "vel": ti.f32}, shape=arrowField) +points.fill(-100) + + +def initPoints(): + for x, y in ti.ndrange(*arrowField): + arrows[x, y].base = vec2( + (x + 1) * (1 - 1 / arrowField[0]) / arrowField[0], + (y + 1) * (1 - 1 / arrowField[1]) / arrowField[1], + ) + dipoles[0].pos = vec2(0.5, 0.5) + dipoles[0].m = 0.01 + vortexes[0].pos = vec2(0.5, 0.5) + vortexes[0].q = -0.5 + + +@ti.func +def getVel(pos): + vel = V # add uniform stream velocity + for i in range(maxElements): + # add source/sink velocity + uv = pos - sources[i].pos + uv[0] *= screen[1] / screen[0] + vel += uv * sources[i].q / (2 * tm.pi * (uv[0] ** 2 + uv[1] ** 2)) + # add vortex velocity + uv = pos - vortexes[i].pos + uv = vec2(-uv[1], uv[0]) + uv[0] *= screen[1] / screen[0] + vel += uv * vortexes[i].q / (2 * tm.pi * (uv[0] ** 2 + uv[1] ** 2)) + # add dipole velocity + uv = pos - dipoles[i].pos + uv[0] *= screen[1] / screen[0] + vel_t = vec2(uv[1] ** 2 - uv[0] ** 2, -2 * uv[0] * uv[1]) + vel += vel_t * dipoles[i].m / (uv[0] ** 2 + uv[1] ** 2) ** 2 + return vel + + +# @param.boundary: 0: left, 1: bottom, 2: right, 3: top +# @param.index: last access point index +@ti.func +def refillPointsOnOneBoundary(boundary, index): + pointsNumber = screen[0] + if boundary == 0 or boundary == 2: + pointsNumber = screen[1] + for i in range(pointsNumber): + found = False + for _ in range(maxPoints): + if points[index][0] == -100: + found = True + if boundary == 0: + points[index] = vec2(-dt * refillVelThresh, (i + 0.5) / screen[1]) + elif boundary == 1: + points[index] = vec2((i + 0.5) / screen[0], -dt * refillVelThresh) + elif boundary == 2: + points[index] = vec2(1 + dt * refillVelThresh, (i + 0.5) / screen[1]) + elif boundary == 3: + points[index] = vec2((i + 0.5) / screen[0], 1 + dt * refillVelThresh) + break + index += 1 + if index >= maxPoints: + index = 0 + if not found: + break + return index + + +@ti.kernel +def refillPoints(): + # traverse positions and refill points on the boundary. + # if the normal velocity is less than thresh, refill point will be deleted when update. + index = 0 + found = True + ti.loop_config(serialize=True) + for i in range(4): + _index = refillPointsOnOneBoundary(i, index) + if _index == index: + break + else: + index = _index + + +@ti.kernel +def updatePoints(): + for i in points: + if points[i][0] != -100: + vel = getVel(points[i]) + points[i] += vel * dt + if not (0 <= points[i][0] <= 1 and 0 <= points[i][1] <= 1): + points[i] = vec2(-100, -100) + for j in range(maxElements): + if sources[j].q < 0 and tm.length(points[i] - sources[j].pos) < 0.025: + points[i] = vec2(-100, -100) + if dipoles[j].m != 0 and tm.length(points[i] - dipoles[j].pos) < 0.05: + points[i] = vec2(-100, -100) + + +@ti.kernel +def updateArrows(): + for x, y in ti.ndrange(*arrowField): + vel = getVel(arrows[x, y].base) + arrows[x, y].vel = vel.norm() + arrows[x, y].dir = vel / vel.norm() / meshSpace / 1.5 + + +def drawArrows(_gui): + global fade + if fade < fadeMax: + fade += 1 + if fade == 0: + updateArrows() # after fade out, update arrow filed + else: + arr = arrows.to_numpy() + vel = arr["vel"].reshape(1, -1)[0] + vel = (vel / vel.max() * 0xDD + 0x11) * (abs(fade / fadeMax)) + mean = vel.mean() + if mean > 0x7F: + vel /= mean / 0x7F # make uniform stream more beautiful + vel = vel.astype(int) + vel *= 2**16 + 2**8 + 1 + _gui.arrows( + arr["base"].reshape(arr["base"].shape[0] * arr["base"].shape[1], 2), + arr["dir"].reshape(arr["dir"].shape[0] * arr["dir"].shape[1], 2), + radius=1.5, + color=vel, + ) + + +def drawMark(_gui, _frame): + triangleTrans = [ + vec2(0, 1) / guiHeight, + vec2(tm.cos(7.0 / 6.0 * tm.pi), tm.sin(7.0 / 6.0 * tm.pi)) / guiHeight, + vec2(tm.cos(-1.0 / 6.0 * tm.pi), tm.sin(-1.0 / 6.0 * tm.pi)) / guiHeight, + ] + rectTrans = [ + vec2(1 * screen[1] / screen[0], 1) / guiHeight, + vec2(-1 * screen[1] / screen[0], -1) / guiHeight, + ] + for i in range(maxElements): + if dipoles[i].m > 0: + _gui.circle(dipoles[i].pos, 0xFFFDC0, dipoles[i].m * 2000) + elif dipoles[i].m < 0: + _gui.circle(dipoles[i].pos, 0xD25565, dipoles[i].m * -2000) + if sources[i].q > 0: + _gui.rect( + sources[i].pos + rectTrans[0] * 2 * sources[i].q, + sources[i].pos + rectTrans[1] * 2 * sources[i].q, + 2 * sources[i].q + 1, + 0xFFFDC0, + ) + elif sources[i].q < 0: + _gui.rect( + sources[i].pos + rectTrans[0] * 2 * sources[i].q, + sources[i].pos + rectTrans[1] * 2 * sources[i].q, + -2 * sources[i].q + 1, + 0xD25565, + ) + if vortexes[i].q != 0: + theta = vortexes[i].q * dt * 40 * _frame + ct, st = ti.cos(theta), ti.sin(theta) + rotateMatrix = ti.Matrix([[ct, -st], [st, ct]]) + trans = [rotateMatrix @ it for it in triangleTrans] + for it in trans: + it[0] *= screen[1] / screen[0] + _gui.triangle( + vortexes[i].pos + trans[0] * 16, + vortexes[i].pos + trans[1] * 16, + vortexes[i].pos + trans[2] * 16, + 0xD25565, + ) + + +def processGuiEvent(_gui): + global fade + while _gui.get_event((ti.GUI.PRESS, ti.GUI.LMB), (ti.GUI.PRESS, ti.GUI.RMB)): + if _gui.is_pressed(ti.GUI.LMB) and _gui.is_pressed(ti.GUI.RMB): # process delete event + for i in range(maxElements): + if sources[i].q != 0 and (sources[i].pos - vec2(*_gui.get_cursor_pos())).norm() < 5 / guiHeight: + sources[i].q = 0 + if vortexes[i].q != 0 and (vortexes[i].pos - vec2(*_gui.get_cursor_pos())).norm() < 5 / guiHeight: + vortexes[i].q = 0 + if dipoles[i].m != 0 and (dipoles[i].pos - vec2(*_gui.get_cursor_pos())).norm() < 5 / guiHeight: + dipoles[i].m = 0 + else: + if _gui.is_pressed("s"): # process add source/sink event + for i in range(maxElements): + if sources[i].q == 0: + if _gui.is_pressed(ti.GUI.RMB): + sources[i].q -= 1 + else: + sources[i].q += 1 + sources[i].pos = vec2(*_gui.get_cursor_pos()) + break + elif _gui.is_pressed("v"): # process add vortex event + for i in range(maxElements): + if vortexes[i].q == 0: + if _gui.is_pressed(ti.GUI.RMB): + vortexes[i].q -= 0.5 + else: + vortexes[i].q += 0.5 + vortexes[i].pos = vec2(*_gui.get_cursor_pos()) + break + elif _gui.is_pressed("d"): # process add dipole event + for i in range(maxElements): + if dipoles[i].m == 0: + if _gui.is_pressed(ti.GUI.RMB): + dipoles[i].m -= 0.01 + else: + dipoles[i].m += 0.01 + dipoles[i].pos = vec2(*_gui.get_cursor_pos()) + break + else: + for i in range(maxElements): # process increase/decrease event + if sources[i].q != 0 and (sources[i].pos - vec2(*_gui.get_cursor_pos())).norm() < 5 / guiHeight: + if _gui.is_pressed(ti.GUI.RMB): + sources[i].q -= 0.5 * int((sources[i].q >= 0.0) - (sources[i].q <= 0.0)) + else: + sources[i].q += 0.5 * int((sources[i].q >= 0.0) - (sources[i].q <= 0.0)) + if vortexes[i].q != 0 and (vortexes[i].pos - vec2(*_gui.get_cursor_pos())).norm() < 5 / guiHeight: + if _gui.is_pressed(ti.GUI.RMB): + vortexes[i].q -= 0.1 * int((vortexes[i].q >= 0.0) - (vortexes[i].q <= 0.0)) + else: + vortexes[i].q += 0.1 * int((vortexes[i].q >= 0.0) - (vortexes[i].q <= 0.0)) + if dipoles[i].m != 0 and (dipoles[i].pos - vec2(*_gui.get_cursor_pos())).norm() < 5 / guiHeight: + if _gui.is_pressed(ti.GUI.RMB): + dipoles[i].m -= 0.001 * int((dipoles[i].m >= 0.0) - (dipoles[i].m <= 0.0)) + else: + dipoles[i].m += 0.001 * int((dipoles[i].m >= 0.0) - (dipoles[i].m <= 0.0)) + fade = -abs(fade) # fade out arrow field + + +if __name__ == "__main__": + initPoints() + updateArrows() + refillPoints() + refillCount = 0 + frame_count = 0 + while gui_.running: + processGuiEvent(gui_) + drawArrows(gui_) + updatePoints() + ps = points.to_numpy() + gui_.circles(ps, 3, color=0x2E94B9) + if refillCount > refillFrame: + refillCount = 0 + refillPoints() + drawMark(gui_, frame_count) + gui_.show() + refillCount += 1 + frame_count += 1 diff --git a/modules/createEVENT/TaichiEvent/mandelbrot_zoom.py b/modules/createEVENT/TaichiEvent/mandelbrot_zoom.py new file mode 100644 index 000000000..8943ea43f --- /dev/null +++ b/modules/createEVENT/TaichiEvent/mandelbrot_zoom.py @@ -0,0 +1,55 @@ +import taichi as ti +from taichi.math import cmul, dot, log2, vec2, vec3 + +ti.init(arch=ti.gpu) + +MAXITERS = 100 +width, height = 800, 640 +pixels = ti.Vector.field(3, ti.f32, shape=(width, height)) + + +@ti.func +def setcolor(z, i): + v = log2(i + 1 - log2(log2(z.norm()))) / 5 + col = vec3(0.0) + if v < 1.0: + col = vec3(v**4, v**2.5, v) + else: + v = ti.max(0.0, 2 - v) + col = vec3(v, v**1.5, v**3) + return col + + +@ti.kernel +def render(time: ti.f32): + zoo = 0.64 + 0.36 * ti.cos(0.02 * time) + zoo = ti.pow(zoo, 8.0) + ca = ti.cos(0.15 * (1.0 - zoo) * time) + sa = ti.sin(0.15 * (1.0 - zoo) * time) + for i, j in pixels: + c = 2.0 * vec2(i, j) / height - vec2(1) + # c *= 1.16 + xy = vec2(c.x * ca - c.y * sa, c.x * sa + c.y * ca) + c = vec2(-0.745, 0.186) + xy * zoo + z = vec2(0.0) + count = 0.0 + while count < MAXITERS and dot(z, z) < 50: + z = cmul(z, z) + c + count += 1.0 + + if count == MAXITERS: + pixels[i, j] = [0, 0, 0] + else: + pixels[i, j] = setcolor(z, count) + + +def main(): + gui = ti.GUI("Mandelbrot set zoom", res=(width, height)) + for i in range(100000): + render(i * 0.03) + gui.set_image(pixels) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/mass_spring_game.py b/modules/createEVENT/TaichiEvent/mass_spring_game.py new file mode 100644 index 000000000..8da3c6b47 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/mass_spring_game.py @@ -0,0 +1,184 @@ +# Tutorials (Chinese): +# - https://www.bilibili.com/video/BV1UK4y177iH +# - https://www.bilibili.com/video/BV1DK411A771 + +import taichi as ti + +ti.init(arch=ti.cpu) + +spring_Y = ti.field(dtype=ti.f32, shape=()) # Young's modulus +paused = ti.field(dtype=ti.i32, shape=()) +drag_damping = ti.field(dtype=ti.f32, shape=()) +dashpot_damping = ti.field(dtype=ti.f32, shape=()) + +max_num_particles = 1024 +particle_mass = 1.0 +dt = 1e-3 +substeps = 10 + +num_particles = ti.field(dtype=ti.i32, shape=()) +x = ti.Vector.field(2, dtype=ti.f32, shape=max_num_particles) +v = ti.Vector.field(2, dtype=ti.f32, shape=max_num_particles) +f = ti.Vector.field(2, dtype=ti.f32, shape=max_num_particles) +fixed = ti.field(dtype=ti.i32, shape=max_num_particles) + +# rest_length[i, j] == 0 means i and j are NOT connected +rest_length = ti.field(dtype=ti.f32, shape=(max_num_particles, max_num_particles)) + + +@ti.kernel +def substep(): + n = num_particles[None] + + # Compute force + for i in range(n): + # Gravity + f[i] = ti.Vector([0, -9.8]) * particle_mass + for j in range(n): + if rest_length[i, j] != 0: + x_ij = x[i] - x[j] + d = x_ij.normalized() + + # Spring force + f[i] += -spring_Y[None] * (x_ij.norm() / rest_length[i, j] - 1) * d + + # Dashpot damping + v_rel = (v[i] - v[j]).dot(d) + f[i] += -dashpot_damping[None] * v_rel * d + + # We use a semi-implicit Euler (aka symplectic Euler) time integrator + for i in range(n): + if not fixed[i]: + v[i] += dt * f[i] / particle_mass + v[i] *= ti.exp(-dt * drag_damping[None]) # Drag damping + + x[i] += v[i] * dt + else: + v[i] = ti.Vector([0, 0]) + + # Collide with four walls + for d in ti.static(range(2)): + # d = 0: treating X (horizontal) component + # d = 1: treating Y (vertical) component + + if x[i][d] < 0: # Bottom and left + x[i][d] = 0 # move particle inside + v[i][d] = 0 # stop it from moving further + + if x[i][d] > 1: # Top and right + x[i][d] = 1 # move particle inside + v[i][d] = 0 # stop it from moving further + + +@ti.kernel +def new_particle(pos_x: ti.f32, pos_y: ti.f32, fixed_: ti.i32): + # Taichi doesn't support using vectors as kernel arguments yet, so we pass scalars + new_particle_id = num_particles[None] + x[new_particle_id] = [pos_x, pos_y] + v[new_particle_id] = [0, 0] + fixed[new_particle_id] = fixed_ + num_particles[None] += 1 + + # Connect with existing particles + for i in range(new_particle_id): + dist = (x[new_particle_id] - x[i]).norm() + connection_radius = 0.15 + if dist < connection_radius: + # Connect the new particle with particle i + rest_length[i, new_particle_id] = 0.1 + rest_length[new_particle_id, i] = 0.1 + + +@ti.kernel +def attract(pos_x: ti.f32, pos_y: ti.f32): + for i in range(num_particles[None]): + p = ti.Vector([pos_x, pos_y]) + v[i] += -dt * substeps * (x[i] - p) * 100 + + +def main(): + gui = ti.GUI("Explicit Mass Spring System", res=(512, 512), background_color=0xDDDDDD) + + spring_Y[None] = 1000 + drag_damping[None] = 1 + dashpot_damping[None] = 100 + + new_particle(0.3, 0.3, False) + new_particle(0.3, 0.4, False) + new_particle(0.4, 0.4, False) + + while True: + for e in gui.get_events(ti.GUI.PRESS): + if e.key in [ti.GUI.ESCAPE, ti.GUI.EXIT]: + exit() + elif e.key == gui.SPACE: + paused[None] = not paused[None] + elif e.key == ti.GUI.LMB: + new_particle(e.pos[0], e.pos[1], int(gui.is_pressed(ti.GUI.SHIFT))) + elif e.key == "c": + num_particles[None] = 0 + rest_length.fill(0) + elif e.key == "y": + if gui.is_pressed("Shift"): + spring_Y[None] /= 1.1 + else: + spring_Y[None] *= 1.1 + elif e.key == "d": + if gui.is_pressed("Shift"): + drag_damping[None] /= 1.1 + else: + drag_damping[None] *= 1.1 + elif e.key == "x": + if gui.is_pressed("Shift"): + dashpot_damping[None] /= 1.1 + else: + dashpot_damping[None] *= 1.1 + + if gui.is_pressed(ti.GUI.RMB): + cursor_pos = gui.get_cursor_pos() + attract(cursor_pos[0], cursor_pos[1]) + + if not paused[None]: + for step in range(substeps): + substep() + + X = x.to_numpy() + n = num_particles[None] + + # Draw the springs + for i in range(n): + for j in range(i + 1, n): + if rest_length[i, j] != 0: + gui.line(begin=X[i], end=X[j], radius=2, color=0x444444) + + # Draw the particles + for i in range(n): + c = 0xFF0000 if fixed[i] else 0x111111 + gui.circle(pos=X[i], color=c, radius=5) + + gui.text( + content="Left click: add mass point (with shift to fix); Right click: attract", + pos=(0, 0.99), + color=0x0, + ) + gui.text(content="C: clear all; Space: pause", pos=(0, 0.95), color=0x0) + gui.text( + content=f"Y: Spring Young's modulus {spring_Y[None]:.1f}", + pos=(0, 0.9), + color=0x0, + ) + gui.text( + content=f"D: Drag damping {drag_damping[None]:.2f}", + pos=(0, 0.85), + color=0x0, + ) + gui.text( + content=f"X: Dashpot damping {dashpot_damping[None]:.2f}", + pos=(0, 0.8), + color=0x0, + ) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/mpm128.py b/modules/createEVENT/TaichiEvent/mpm128.py new file mode 100644 index 000000000..83b013d03 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/mpm128.py @@ -0,0 +1,159 @@ +import taichi as ti + +ti.init(arch=ti.gpu) # Try to run on GPU + +quality = 1 # Use a larger value for higher-res simulations +n_particles, n_grid = 9000 * quality**2, 128 * quality +dx, inv_dx = 1 / n_grid, float(n_grid) +dt = 1e-4 / quality +p_vol, p_rho = (dx * 0.5) ** 2, 1 +p_mass = p_vol * p_rho +E, nu = 5e3, 0.2 # Young's modulus and Poisson's ratio +mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / ((1 + nu) * (1 - 2 * nu)) # Lame parameters + +x = ti.Vector.field(2, dtype=float, shape=n_particles) # position +v = ti.Vector.field(2, dtype=float, shape=n_particles) # velocity +C = ti.Matrix.field(2, 2, dtype=float, shape=n_particles) # affine velocity field +F = ti.Matrix.field(2, 2, dtype=float, shape=n_particles) # deformation gradient +material = ti.field(dtype=int, shape=n_particles) # material id +Jp = ti.field(dtype=float, shape=n_particles) # plastic deformation +grid_v = ti.Vector.field(2, dtype=float, shape=(n_grid, n_grid)) # grid node momentum/velocity +grid_m = ti.field(dtype=float, shape=(n_grid, n_grid)) # grid node mass +gravity = ti.Vector.field(2, dtype=float, shape=()) +attractor_strength = ti.field(dtype=float, shape=()) +attractor_pos = ti.Vector.field(2, dtype=float, shape=()) + + +@ti.kernel +def substep(): + for i, j in grid_m: + grid_v[i, j] = [0, 0] + grid_m[i, j] = 0 + for p in x: # Particle state update and scatter to grid (P2G) + base = (x[p] * inv_dx - 0.5).cast(int) + fx = x[p] * inv_dx - base.cast(float) + # Quadratic kernels [http://mpm.graphics Eqn. 123, with x=fx, fx-1,fx-2] + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2] + # deformation gradient update + F[p] = (ti.Matrix.identity(float, 2) + dt * C[p]) @ F[p] + # Hardening coefficient: snow gets harder when compressed + h = ti.max(0.1, ti.min(5, ti.exp(10 * (1.0 - Jp[p])))) + if material[p] == 1: # jelly, make it softer + h = 0.3 + mu, la = mu_0 * h, lambda_0 * h + if material[p] == 0: # liquid + mu = 0.0 + U, sig, V = ti.svd(F[p]) + J = 1.0 + for d in ti.static(range(2)): + new_sig = sig[d, d] + if material[p] == 2: # Snow + new_sig = min(max(sig[d, d], 1 - 2.5e-2), 1 + 4.5e-3) # Plasticity + Jp[p] *= sig[d, d] / new_sig + sig[d, d] = new_sig + J *= new_sig + if material[p] == 0: + # Reset deformation gradient to avoid numerical instability + F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J) + elif material[p] == 2: + # Reconstruct elastic deformation gradient after plasticity + F[p] = U @ sig @ V.transpose() + stress = 2 * mu * (F[p] - U @ V.transpose()) @ F[p].transpose() + ti.Matrix.identity(float, 2) * la * J * ( + J - 1 + ) + stress = (-dt * p_vol * 4 * inv_dx * inv_dx) * stress + affine = stress + p_mass * C[p] + for i, j in ti.static(ti.ndrange(3, 3)): + # Loop over 3x3 grid node neighborhood + offset = ti.Vector([i, j]) + dpos = (offset.cast(float) - fx) * dx + weight = w[i][0] * w[j][1] + grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos) + grid_m[base + offset] += weight * p_mass + for i, j in grid_m: + if grid_m[i, j] > 0: # No need for epsilon here + # Momentum to velocity + grid_v[i, j] = (1 / grid_m[i, j]) * grid_v[i, j] + grid_v[i, j] += dt * gravity[None] * 30 # gravity + dist = attractor_pos[None] - dx * ti.Vector([i, j]) + grid_v[i, j] += dist / (0.01 + dist.norm()) * attractor_strength[None] * dt * 100 + if i < 3 and grid_v[i, j][0] < 0: + grid_v[i, j][0] = 0 # Boundary conditions + if i > n_grid - 3 and grid_v[i, j][0] > 0: + grid_v[i, j][0] = 0 + if j < 3 and grid_v[i, j][1] < 0: + grid_v[i, j][1] = 0 + if j > n_grid - 3 and grid_v[i, j][1] > 0: + grid_v[i, j][1] = 0 + for p in x: # grid to particle (G2P) + base = (x[p] * inv_dx - 0.5).cast(int) + fx = x[p] * inv_dx - base.cast(float) + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1.0) ** 2, 0.5 * (fx - 0.5) ** 2] + new_v = ti.Vector.zero(float, 2) + new_C = ti.Matrix.zero(float, 2, 2) + for i, j in ti.static(ti.ndrange(3, 3)): + # loop over 3x3 grid node neighborhood + dpos = ti.Vector([i, j]).cast(float) - fx + g_v = grid_v[base + ti.Vector([i, j])] + weight = w[i][0] * w[j][1] + new_v += weight * g_v + new_C += 4 * inv_dx * weight * g_v.outer_product(dpos) + v[p], C[p] = new_v, new_C + x[p] += dt * v[p] # advection + + +@ti.kernel +def reset(): + group_size = n_particles // 3 + for i in range(n_particles): + x[i] = [ + ti.random() * 0.2 + 0.3 + 0.10 * (i // group_size), + ti.random() * 0.2 + 0.05 + 0.32 * (i // group_size), + ] + material[i] = i // group_size # 0: fluid 1: jelly 2: snow + v[i] = [0, 0] + F[i] = ti.Matrix([[1, 0], [0, 1]]) + Jp[i] = 1 + C[i] = ti.Matrix.zero(float, 2, 2) + + +print("[Hint] Use WSAD/arrow keys to control gravity. Use left/right mouse buttons to attract/repel. Press R to reset.") +gui = ti.GUI("Taichi MLS-MPM-128", res=512, background_color=0x112F41) +reset() +gravity[None] = [0, -1] + +for frame in range(20000): + if gui.get_event(ti.GUI.PRESS): + if gui.event.key == "r": + reset() + elif gui.event.key in [ti.GUI.ESCAPE, ti.GUI.EXIT]: + break + if gui.event is not None: + gravity[None] = [0, 0] # if had any event + if gui.is_pressed(ti.GUI.LEFT, "a"): + gravity[None][0] = -1 + if gui.is_pressed(ti.GUI.RIGHT, "d"): + gravity[None][0] = 1 + if gui.is_pressed(ti.GUI.UP, "w"): + gravity[None][1] = 1 + if gui.is_pressed(ti.GUI.DOWN, "s"): + gravity[None][1] = -1 + mouse = gui.get_cursor_pos() + gui.circle((mouse[0], mouse[1]), color=0x336699, radius=15) + attractor_pos[None] = [mouse[0], mouse[1]] + attractor_strength[None] = 0 + if gui.is_pressed(ti.GUI.LMB): + attractor_strength[None] = 1 + if gui.is_pressed(ti.GUI.RMB): + attractor_strength[None] = -1 + for s in range(int(2e-3 // dt)): + substep() + gui.circles( + x.to_numpy(), + radius=1.5, + palette=[0x068587, 0xED553B, 0xEEEEF0], + palette_indices=material, + ) + + # Change to gui.show(f'{frame:06d}.png') to write images to disk + gui.show() diff --git a/modules/createEVENT/TaichiEvent/mpm3d.py b/modules/createEVENT/TaichiEvent/mpm3d.py new file mode 100644 index 000000000..d1a1e2bbc --- /dev/null +++ b/modules/createEVENT/TaichiEvent/mpm3d.py @@ -0,0 +1,122 @@ +export_file = "" # use '/tmp/mpm3d.ply' for exporting result to disk + +import numpy as np + +import taichi as ti + +ti.init(arch=ti.gpu) + +# dim, n_grid, steps, dt = 2, 128, 20, 2e-4 +# dim, n_grid, steps, dt = 2, 256, 32, 1e-4 +dim, n_grid, steps, dt = 3, 32, 25, 4e-4 +# dim, n_grid, steps, dt = 3, 64, 25, 2e-4 +# dim, n_grid, steps, dt = 3, 128, 25, 8e-5 + +n_particles = n_grid**dim // 2 ** (dim - 1) +dx = 1 / n_grid + +p_rho = 1 +p_vol = (dx * 0.5) ** 2 +p_mass = p_vol * p_rho +gravity = 9.8 +bound = 3 +E = 400 + +F_x = ti.Vector.field(dim, float, n_particles) +F_v = ti.Vector.field(dim, float, n_particles) +F_C = ti.Matrix.field(dim, dim, float, n_particles) +F_J = ti.field(float, n_particles) + +F_grid_v = ti.Vector.field(dim, float, (n_grid,) * dim) +F_grid_m = ti.field(float, (n_grid,) * dim) + +neighbour = (3,) * dim + + +@ti.kernel +def substep(): + for I in ti.grouped(F_grid_m): + F_grid_v[I] = ti.zero(F_grid_v[I]) + F_grid_m[I] = 0 + ti.loop_config(block_dim=n_grid) + for p in F_x: + Xp = F_x[p] / dx + base = int(Xp - 0.5) + fx = Xp - base + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2] + stress = -dt * 4 * E * p_vol * (F_J[p] - 1) / dx**2 + affine = ti.Matrix.identity(float, dim) * stress + p_mass * F_C[p] + for offset in ti.static(ti.grouped(ti.ndrange(*neighbour))): + dpos = (offset - fx) * dx + weight = 1.0 + for i in ti.static(range(dim)): + weight *= w[offset[i]][i] + F_grid_v[base + offset] += weight * (p_mass * F_v[p] + affine @ dpos) + F_grid_m[base + offset] += weight * p_mass + for I in ti.grouped(F_grid_m): + if F_grid_m[I] > 0: + F_grid_v[I] /= F_grid_m[I] + F_grid_v[I][1] -= dt * gravity + cond = (I < bound) & (F_grid_v[I] < 0) | (I > n_grid - bound) & (F_grid_v[I] > 0) + F_grid_v[I] = ti.select(cond, 0, F_grid_v[I]) + ti.loop_config(block_dim=n_grid) + for p in F_x: + Xp = F_x[p] / dx + base = int(Xp - 0.5) + fx = Xp - base + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2] + new_v = ti.zero(F_v[p]) + new_C = ti.zero(F_C[p]) + for offset in ti.static(ti.grouped(ti.ndrange(*neighbour))): + dpos = (offset - fx) * dx + weight = 1.0 + for i in ti.static(range(dim)): + weight *= w[offset[i]][i] + g_v = F_grid_v[base + offset] + new_v += weight * g_v + new_C += 4 * weight * g_v.outer_product(dpos) / dx**2 + F_v[p] = new_v + F_x[p] += dt * F_v[p] + F_J[p] *= 1 + dt * new_C.trace() + F_C[p] = new_C + + +@ti.kernel +def init(): + for i in range(n_particles): + F_x[i] = ti.Vector([ti.random() for i in range(dim)]) * 0.4 + 0.15 + F_J[i] = 1 + + +def T(a): + if dim == 2: + return a + + phi, theta = np.radians(28), np.radians(32) + + a = a - 0.5 + x, y, z = a[:, 0], a[:, 1], a[:, 2] + cp, sp = np.cos(phi), np.sin(phi) + ct, st = np.cos(theta), np.sin(theta) + x, z = x * cp + z * sp, z * cp - x * sp + u, v = x, y * ct + z * st + return np.array([u, v]).swapaxes(0, 1) + 0.5 + + +def main(): + init() + gui = ti.GUI("MPM3D", background_color=0x112F41) + while gui.running and not gui.get_event(gui.ESCAPE): + for s in range(steps): + substep() + pos = F_x.to_numpy() + if export_file: + writer = ti.tools.PLYWriter(num_vertices=n_particles) + writer.add_vertex_pos(pos[:, 0], pos[:, 1], pos[:, 2]) + writer.export_frame(gui.frame, export_file) + gui.circles(T(pos), radius=1.5, color=0x66CCFF) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/mpm88.py b/modules/createEVENT/TaichiEvent/mpm88.py new file mode 100644 index 000000000..875d5386f --- /dev/null +++ b/modules/createEVENT/TaichiEvent/mpm88.py @@ -0,0 +1,92 @@ +# MPM-MLS in 88 lines of Taichi code, originally created by @yuanming-hu +import taichi as ti + +ti.init(arch=ti.gpu) + +n_particles = 8192 +n_grid = 128 +dx = 1 / n_grid +dt = 2e-4 + +p_rho = 1 +p_vol = (dx * 0.5) ** 2 +p_mass = p_vol * p_rho +gravity = 9.8 +bound = 3 +E = 400 + +x = ti.Vector.field(2, float, n_particles) +v = ti.Vector.field(2, float, n_particles) +C = ti.Matrix.field(2, 2, float, n_particles) +J = ti.field(float, n_particles) + +grid_v = ti.Vector.field(2, float, (n_grid, n_grid)) +grid_m = ti.field(float, (n_grid, n_grid)) + + +@ti.kernel +def substep(): + for i, j in grid_m: + grid_v[i, j] = [0, 0] + grid_m[i, j] = 0 + for p in x: + Xp = x[p] / dx + base = int(Xp - 0.5) + fx = Xp - base + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2] + stress = -dt * 4 * E * p_vol * (J[p] - 1) / dx**2 + affine = ti.Matrix([[stress, 0], [0, stress]]) + p_mass * C[p] + for i, j in ti.static(ti.ndrange(3, 3)): + offset = ti.Vector([i, j]) + dpos = (offset - fx) * dx + weight = w[i].x * w[j].y + grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos) + grid_m[base + offset] += weight * p_mass + for i, j in grid_m: + if grid_m[i, j] > 0: + grid_v[i, j] /= grid_m[i, j] + grid_v[i, j].y -= dt * gravity + if i < bound and grid_v[i, j].x < 0: + grid_v[i, j].x = 0 + if i > n_grid - bound and grid_v[i, j].x > 0: + grid_v[i, j].x = 0 + if j < bound and grid_v[i, j].y < 0: + grid_v[i, j].y = 0 + if j > n_grid - bound and grid_v[i, j].y > 0: + grid_v[i, j].y = 0 + for p in x: + Xp = x[p] / dx + base = int(Xp - 0.5) + fx = Xp - base + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2] + new_v = ti.Vector.zero(float, 2) + new_C = ti.Matrix.zero(float, 2, 2) + for i, j in ti.static(ti.ndrange(3, 3)): + offset = ti.Vector([i, j]) + dpos = (offset - fx) * dx + weight = w[i].x * w[j].y + g_v = grid_v[base + offset] + new_v += weight * g_v + new_C += 4 * weight * g_v.outer_product(dpos) / dx**2 + v[p] = new_v + x[p] += dt * v[p] + J[p] *= 1 + dt * new_C.trace() + C[p] = new_C + + +@ti.kernel +def init(): + for i in range(n_particles): + x[i] = [ti.random() * 0.4 + 0.2, ti.random() * 0.4 + 0.2] + v[i] = [0, -1] + J[i] = 1 + + +init() +gui = ti.GUI("MPM88") +while gui.running and not gui.get_event(gui.ESCAPE): + for s in range(50): + substep() + gui.clear(0x112F41) + gui.circles(x.to_numpy(), radius=1.5, color=0x068587) + gui.show() diff --git a/modules/createEVENT/TaichiEvent/mpm99.py b/modules/createEVENT/TaichiEvent/mpm99.py new file mode 100644 index 000000000..99b157cd4 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/mpm99.py @@ -0,0 +1,133 @@ +import taichi as ti + +ti.init(arch=ti.gpu) # Try to run on GPU +quality = 1 # Use a larger value for higher-res simulations +n_particles, n_grid = 9000 * quality**2, 128 * quality +dx, inv_dx = 1 / n_grid, float(n_grid) +dt = 1e-4 / quality +p_vol, p_rho = (dx * 0.5) ** 2, 1 +p_mass = p_vol * p_rho +E, nu = 0.1e4, 0.2 # Young's modulus and Poisson's ratio +mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / ((1 + nu) * (1 - 2 * nu)) # Lame parameters +x = ti.Vector.field(2, dtype=float, shape=n_particles) # position +v = ti.Vector.field(2, dtype=float, shape=n_particles) # velocity +C = ti.Matrix.field(2, 2, dtype=float, shape=n_particles) # affine velocity field +F = ti.Matrix.field(2, 2, dtype=float, shape=n_particles) # deformation gradient +material = ti.field(dtype=int, shape=n_particles) # material id +Jp = ti.field(dtype=float, shape=n_particles) # plastic deformation +grid_v = ti.Vector.field(2, dtype=float, shape=(n_grid, n_grid)) # grid node momentum/velocity +grid_m = ti.field(dtype=float, shape=(n_grid, n_grid)) # grid node mass + + +@ti.kernel +def substep(): + for i, j in grid_m: + grid_v[i, j] = [0, 0] + grid_m[i, j] = 0 + for p in x: # Particle state update and scatter to grid (P2G) + base = (x[p] * inv_dx - 0.5).cast(int) + fx = x[p] * inv_dx - base.cast(float) + # Quadratic kernels [http://mpm.graphics Eqn. 123, with x=fx, fx-1,fx-2] + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2] + # F[p]: deformation gradient update + F[p] = (ti.Matrix.identity(float, 2) + dt * C[p]) @ F[p] + # h: Hardening coefficient: snow gets harder when compressed + h = ti.exp(10 * (1.0 - Jp[p])) + if material[p] == 1: # jelly, make it softer + h = 0.3 + mu, la = mu_0 * h, lambda_0 * h + if material[p] == 0: # liquid + mu = 0.0 + U, sig, V = ti.svd(F[p]) + # Avoid zero eigenvalues because of numerical errors + for d in ti.static(range(2)): + sig[d, d] = ti.max(sig[d, d], 1e-6) + J = 1.0 + for d in ti.static(range(2)): + new_sig = sig[d, d] + if material[p] == 2: # Snow + new_sig = ti.min(ti.max(sig[d, d], 1 - 2.5e-2), 1 + 4.5e-3) # Plasticity + Jp[p] *= sig[d, d] / new_sig + sig[d, d] = new_sig + J *= new_sig + if material[p] == 0: + # Reset deformation gradient to avoid numerical instability + F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J) + elif material[p] == 2: + # Reconstruct elastic deformation gradient after plasticity + F[p] = U @ sig @ V.transpose() + stress = 2 * mu * (F[p] - U @ V.transpose()) @ F[p].transpose() + ti.Matrix.identity(float, 2) * la * J * ( + J - 1 + ) + stress = (-dt * p_vol * 4 * inv_dx * inv_dx) * stress + affine = stress + p_mass * C[p] + # Loop over 3x3 grid node neighborhood + for i, j in ti.static(ti.ndrange(3, 3)): + offset = ti.Vector([i, j]) + dpos = (offset.cast(float) - fx) * dx + weight = w[i][0] * w[j][1] + grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos) + grid_m[base + offset] += weight * p_mass + for i, j in grid_m: + if grid_m[i, j] > 0: # No need for epsilon here + grid_v[i, j] = (1 / grid_m[i, j]) * grid_v[i, j] # Momentum to velocity + grid_v[i, j][1] -= dt * 50 # gravity + if i < 3 and grid_v[i, j][0] < 0: + grid_v[i, j][0] = 0 # Boundary conditions + if i > n_grid - 3 and grid_v[i, j][0] > 0: + grid_v[i, j][0] = 0 + if j < 3 and grid_v[i, j][1] < 0: + grid_v[i, j][1] = 0 + if j > n_grid - 3 and grid_v[i, j][1] > 0: + grid_v[i, j][1] = 0 + for p in x: # grid to particle (G2P) + base = (x[p] * inv_dx - 0.5).cast(int) + fx = x[p] * inv_dx - base.cast(float) + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1.0) ** 2, 0.5 * (fx - 0.5) ** 2] + new_v = ti.Vector.zero(float, 2) + new_C = ti.Matrix.zero(float, 2, 2) + for i, j in ti.static(ti.ndrange(3, 3)): + # loop over 3x3 grid node neighborhood + dpos = ti.Vector([i, j]).cast(float) - fx + g_v = grid_v[base + ti.Vector([i, j])] + weight = w[i][0] * w[j][1] + new_v += weight * g_v + new_C += 4 * inv_dx * weight * g_v.outer_product(dpos) + v[p], C[p] = new_v, new_C + x[p] += dt * v[p] # advection + + +group_size = n_particles // 3 + + +@ti.kernel +def initialize(): + for i in range(n_particles): + x[i] = [ + ti.random() * 0.2 + 0.3 + 0.10 * (i // group_size), + ti.random() * 0.2 + 0.05 + 0.32 * (i // group_size), + ] + material[i] = i // group_size # 0: fluid 1: jelly 2: snow + v[i] = ti.Matrix([0, 0]) + F[i] = ti.Matrix([[1, 0], [0, 1]]) + Jp[i] = 1 + + +def main(): + initialize() + gui = ti.GUI("Taichi MLS-MPM-99", res=512, background_color=0x112F41) + while not gui.get_event(ti.GUI.ESCAPE, ti.GUI.EXIT): + for s in range(int(2e-3 // dt)): + substep() + gui.circles( + x.to_numpy(), + radius=1.5, + palette=[0x068587, 0xED553B, 0xEEEEF0], + palette_indices=material, + ) + # Change to gui.show(f'{frame:06d}.png') to write images to disk + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/mpm_fluid_debri_2d_reu_2024.06.10.py b/modules/createEVENT/TaichiEvent/mpm_fluid_debri_2d_reu_2024.06.10.py new file mode 100644 index 000000000..83b013d03 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/mpm_fluid_debri_2d_reu_2024.06.10.py @@ -0,0 +1,159 @@ +import taichi as ti + +ti.init(arch=ti.gpu) # Try to run on GPU + +quality = 1 # Use a larger value for higher-res simulations +n_particles, n_grid = 9000 * quality**2, 128 * quality +dx, inv_dx = 1 / n_grid, float(n_grid) +dt = 1e-4 / quality +p_vol, p_rho = (dx * 0.5) ** 2, 1 +p_mass = p_vol * p_rho +E, nu = 5e3, 0.2 # Young's modulus and Poisson's ratio +mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / ((1 + nu) * (1 - 2 * nu)) # Lame parameters + +x = ti.Vector.field(2, dtype=float, shape=n_particles) # position +v = ti.Vector.field(2, dtype=float, shape=n_particles) # velocity +C = ti.Matrix.field(2, 2, dtype=float, shape=n_particles) # affine velocity field +F = ti.Matrix.field(2, 2, dtype=float, shape=n_particles) # deformation gradient +material = ti.field(dtype=int, shape=n_particles) # material id +Jp = ti.field(dtype=float, shape=n_particles) # plastic deformation +grid_v = ti.Vector.field(2, dtype=float, shape=(n_grid, n_grid)) # grid node momentum/velocity +grid_m = ti.field(dtype=float, shape=(n_grid, n_grid)) # grid node mass +gravity = ti.Vector.field(2, dtype=float, shape=()) +attractor_strength = ti.field(dtype=float, shape=()) +attractor_pos = ti.Vector.field(2, dtype=float, shape=()) + + +@ti.kernel +def substep(): + for i, j in grid_m: + grid_v[i, j] = [0, 0] + grid_m[i, j] = 0 + for p in x: # Particle state update and scatter to grid (P2G) + base = (x[p] * inv_dx - 0.5).cast(int) + fx = x[p] * inv_dx - base.cast(float) + # Quadratic kernels [http://mpm.graphics Eqn. 123, with x=fx, fx-1,fx-2] + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2] + # deformation gradient update + F[p] = (ti.Matrix.identity(float, 2) + dt * C[p]) @ F[p] + # Hardening coefficient: snow gets harder when compressed + h = ti.max(0.1, ti.min(5, ti.exp(10 * (1.0 - Jp[p])))) + if material[p] == 1: # jelly, make it softer + h = 0.3 + mu, la = mu_0 * h, lambda_0 * h + if material[p] == 0: # liquid + mu = 0.0 + U, sig, V = ti.svd(F[p]) + J = 1.0 + for d in ti.static(range(2)): + new_sig = sig[d, d] + if material[p] == 2: # Snow + new_sig = min(max(sig[d, d], 1 - 2.5e-2), 1 + 4.5e-3) # Plasticity + Jp[p] *= sig[d, d] / new_sig + sig[d, d] = new_sig + J *= new_sig + if material[p] == 0: + # Reset deformation gradient to avoid numerical instability + F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J) + elif material[p] == 2: + # Reconstruct elastic deformation gradient after plasticity + F[p] = U @ sig @ V.transpose() + stress = 2 * mu * (F[p] - U @ V.transpose()) @ F[p].transpose() + ti.Matrix.identity(float, 2) * la * J * ( + J - 1 + ) + stress = (-dt * p_vol * 4 * inv_dx * inv_dx) * stress + affine = stress + p_mass * C[p] + for i, j in ti.static(ti.ndrange(3, 3)): + # Loop over 3x3 grid node neighborhood + offset = ti.Vector([i, j]) + dpos = (offset.cast(float) - fx) * dx + weight = w[i][0] * w[j][1] + grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos) + grid_m[base + offset] += weight * p_mass + for i, j in grid_m: + if grid_m[i, j] > 0: # No need for epsilon here + # Momentum to velocity + grid_v[i, j] = (1 / grid_m[i, j]) * grid_v[i, j] + grid_v[i, j] += dt * gravity[None] * 30 # gravity + dist = attractor_pos[None] - dx * ti.Vector([i, j]) + grid_v[i, j] += dist / (0.01 + dist.norm()) * attractor_strength[None] * dt * 100 + if i < 3 and grid_v[i, j][0] < 0: + grid_v[i, j][0] = 0 # Boundary conditions + if i > n_grid - 3 and grid_v[i, j][0] > 0: + grid_v[i, j][0] = 0 + if j < 3 and grid_v[i, j][1] < 0: + grid_v[i, j][1] = 0 + if j > n_grid - 3 and grid_v[i, j][1] > 0: + grid_v[i, j][1] = 0 + for p in x: # grid to particle (G2P) + base = (x[p] * inv_dx - 0.5).cast(int) + fx = x[p] * inv_dx - base.cast(float) + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1.0) ** 2, 0.5 * (fx - 0.5) ** 2] + new_v = ti.Vector.zero(float, 2) + new_C = ti.Matrix.zero(float, 2, 2) + for i, j in ti.static(ti.ndrange(3, 3)): + # loop over 3x3 grid node neighborhood + dpos = ti.Vector([i, j]).cast(float) - fx + g_v = grid_v[base + ti.Vector([i, j])] + weight = w[i][0] * w[j][1] + new_v += weight * g_v + new_C += 4 * inv_dx * weight * g_v.outer_product(dpos) + v[p], C[p] = new_v, new_C + x[p] += dt * v[p] # advection + + +@ti.kernel +def reset(): + group_size = n_particles // 3 + for i in range(n_particles): + x[i] = [ + ti.random() * 0.2 + 0.3 + 0.10 * (i // group_size), + ti.random() * 0.2 + 0.05 + 0.32 * (i // group_size), + ] + material[i] = i // group_size # 0: fluid 1: jelly 2: snow + v[i] = [0, 0] + F[i] = ti.Matrix([[1, 0], [0, 1]]) + Jp[i] = 1 + C[i] = ti.Matrix.zero(float, 2, 2) + + +print("[Hint] Use WSAD/arrow keys to control gravity. Use left/right mouse buttons to attract/repel. Press R to reset.") +gui = ti.GUI("Taichi MLS-MPM-128", res=512, background_color=0x112F41) +reset() +gravity[None] = [0, -1] + +for frame in range(20000): + if gui.get_event(ti.GUI.PRESS): + if gui.event.key == "r": + reset() + elif gui.event.key in [ti.GUI.ESCAPE, ti.GUI.EXIT]: + break + if gui.event is not None: + gravity[None] = [0, 0] # if had any event + if gui.is_pressed(ti.GUI.LEFT, "a"): + gravity[None][0] = -1 + if gui.is_pressed(ti.GUI.RIGHT, "d"): + gravity[None][0] = 1 + if gui.is_pressed(ti.GUI.UP, "w"): + gravity[None][1] = 1 + if gui.is_pressed(ti.GUI.DOWN, "s"): + gravity[None][1] = -1 + mouse = gui.get_cursor_pos() + gui.circle((mouse[0], mouse[1]), color=0x336699, radius=15) + attractor_pos[None] = [mouse[0], mouse[1]] + attractor_strength[None] = 0 + if gui.is_pressed(ti.GUI.LMB): + attractor_strength[None] = 1 + if gui.is_pressed(ti.GUI.RMB): + attractor_strength[None] = -1 + for s in range(int(2e-3 // dt)): + substep() + gui.circles( + x.to_numpy(), + radius=1.5, + palette=[0x068587, 0xED553B, 0xEEEEF0], + palette_indices=material, + ) + + # Change to gui.show(f'{frame:06d}.png') to write images to disk + gui.show() diff --git a/modules/createEVENT/TaichiEvent/mpm_fluid_debris_2d_reu_nogui_2024.06.10.py b/modules/createEVENT/TaichiEvent/mpm_fluid_debris_2d_reu_nogui_2024.06.10.py new file mode 100644 index 000000000..99b157cd4 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/mpm_fluid_debris_2d_reu_nogui_2024.06.10.py @@ -0,0 +1,133 @@ +import taichi as ti + +ti.init(arch=ti.gpu) # Try to run on GPU +quality = 1 # Use a larger value for higher-res simulations +n_particles, n_grid = 9000 * quality**2, 128 * quality +dx, inv_dx = 1 / n_grid, float(n_grid) +dt = 1e-4 / quality +p_vol, p_rho = (dx * 0.5) ** 2, 1 +p_mass = p_vol * p_rho +E, nu = 0.1e4, 0.2 # Young's modulus and Poisson's ratio +mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / ((1 + nu) * (1 - 2 * nu)) # Lame parameters +x = ti.Vector.field(2, dtype=float, shape=n_particles) # position +v = ti.Vector.field(2, dtype=float, shape=n_particles) # velocity +C = ti.Matrix.field(2, 2, dtype=float, shape=n_particles) # affine velocity field +F = ti.Matrix.field(2, 2, dtype=float, shape=n_particles) # deformation gradient +material = ti.field(dtype=int, shape=n_particles) # material id +Jp = ti.field(dtype=float, shape=n_particles) # plastic deformation +grid_v = ti.Vector.field(2, dtype=float, shape=(n_grid, n_grid)) # grid node momentum/velocity +grid_m = ti.field(dtype=float, shape=(n_grid, n_grid)) # grid node mass + + +@ti.kernel +def substep(): + for i, j in grid_m: + grid_v[i, j] = [0, 0] + grid_m[i, j] = 0 + for p in x: # Particle state update and scatter to grid (P2G) + base = (x[p] * inv_dx - 0.5).cast(int) + fx = x[p] * inv_dx - base.cast(float) + # Quadratic kernels [http://mpm.graphics Eqn. 123, with x=fx, fx-1,fx-2] + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2] + # F[p]: deformation gradient update + F[p] = (ti.Matrix.identity(float, 2) + dt * C[p]) @ F[p] + # h: Hardening coefficient: snow gets harder when compressed + h = ti.exp(10 * (1.0 - Jp[p])) + if material[p] == 1: # jelly, make it softer + h = 0.3 + mu, la = mu_0 * h, lambda_0 * h + if material[p] == 0: # liquid + mu = 0.0 + U, sig, V = ti.svd(F[p]) + # Avoid zero eigenvalues because of numerical errors + for d in ti.static(range(2)): + sig[d, d] = ti.max(sig[d, d], 1e-6) + J = 1.0 + for d in ti.static(range(2)): + new_sig = sig[d, d] + if material[p] == 2: # Snow + new_sig = ti.min(ti.max(sig[d, d], 1 - 2.5e-2), 1 + 4.5e-3) # Plasticity + Jp[p] *= sig[d, d] / new_sig + sig[d, d] = new_sig + J *= new_sig + if material[p] == 0: + # Reset deformation gradient to avoid numerical instability + F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J) + elif material[p] == 2: + # Reconstruct elastic deformation gradient after plasticity + F[p] = U @ sig @ V.transpose() + stress = 2 * mu * (F[p] - U @ V.transpose()) @ F[p].transpose() + ti.Matrix.identity(float, 2) * la * J * ( + J - 1 + ) + stress = (-dt * p_vol * 4 * inv_dx * inv_dx) * stress + affine = stress + p_mass * C[p] + # Loop over 3x3 grid node neighborhood + for i, j in ti.static(ti.ndrange(3, 3)): + offset = ti.Vector([i, j]) + dpos = (offset.cast(float) - fx) * dx + weight = w[i][0] * w[j][1] + grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos) + grid_m[base + offset] += weight * p_mass + for i, j in grid_m: + if grid_m[i, j] > 0: # No need for epsilon here + grid_v[i, j] = (1 / grid_m[i, j]) * grid_v[i, j] # Momentum to velocity + grid_v[i, j][1] -= dt * 50 # gravity + if i < 3 and grid_v[i, j][0] < 0: + grid_v[i, j][0] = 0 # Boundary conditions + if i > n_grid - 3 and grid_v[i, j][0] > 0: + grid_v[i, j][0] = 0 + if j < 3 and grid_v[i, j][1] < 0: + grid_v[i, j][1] = 0 + if j > n_grid - 3 and grid_v[i, j][1] > 0: + grid_v[i, j][1] = 0 + for p in x: # grid to particle (G2P) + base = (x[p] * inv_dx - 0.5).cast(int) + fx = x[p] * inv_dx - base.cast(float) + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1.0) ** 2, 0.5 * (fx - 0.5) ** 2] + new_v = ti.Vector.zero(float, 2) + new_C = ti.Matrix.zero(float, 2, 2) + for i, j in ti.static(ti.ndrange(3, 3)): + # loop over 3x3 grid node neighborhood + dpos = ti.Vector([i, j]).cast(float) - fx + g_v = grid_v[base + ti.Vector([i, j])] + weight = w[i][0] * w[j][1] + new_v += weight * g_v + new_C += 4 * inv_dx * weight * g_v.outer_product(dpos) + v[p], C[p] = new_v, new_C + x[p] += dt * v[p] # advection + + +group_size = n_particles // 3 + + +@ti.kernel +def initialize(): + for i in range(n_particles): + x[i] = [ + ti.random() * 0.2 + 0.3 + 0.10 * (i // group_size), + ti.random() * 0.2 + 0.05 + 0.32 * (i // group_size), + ] + material[i] = i // group_size # 0: fluid 1: jelly 2: snow + v[i] = ti.Matrix([0, 0]) + F[i] = ti.Matrix([[1, 0], [0, 1]]) + Jp[i] = 1 + + +def main(): + initialize() + gui = ti.GUI("Taichi MLS-MPM-99", res=512, background_color=0x112F41) + while not gui.get_event(ti.GUI.ESCAPE, ti.GUI.EXIT): + for s in range(int(2e-3 // dt)): + substep() + gui.circles( + x.to_numpy(), + radius=1.5, + palette=[0x068587, 0xED553B, 0xEEEEF0], + palette_indices=material, + ) + # Change to gui.show(f'{frame:06d}.png') to write images to disk + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/mpm_lagrangian_forces.py b/modules/createEVENT/TaichiEvent/mpm_lagrangian_forces.py index 975f5999b..cbc002f6a 100755 --- a/modules/createEVENT/TaichiEvent/mpm_lagrangian_forces.py +++ b/modules/createEVENT/TaichiEvent/mpm_lagrangian_forces.py @@ -1,204 +1,186 @@ -#!/usr/bin/env python3 # noqa: D100 - -import numpy as np -import taichi as ti - -ti.init(arch=ti.gpu) - -dim = 2 -quality = 1 # Use a larger integral number for higher quality -n_particle_x = 100 * quality -n_particle_y = 8 * quality -n_particles = n_particle_x * n_particle_y -n_elements = (n_particle_x - 1) * (n_particle_y - 1) * 2 -n_grid = 64 * quality -dx = 1 / n_grid -inv_dx = 1 / dx -dt = 1e-4 / quality -E = 25000 -p_mass = 1 -p_vol = 1 -mu = 1 -la = 1 - -x = ti.Vector.field(dim, dtype=float, shape=n_particles, needs_grad=True) -v = ti.Vector.field(dim, dtype=float, shape=n_particles) -C = ti.Matrix.field(dim, dim, dtype=float, shape=n_particles) -grid_v = ti.Vector.field(dim, dtype=float, shape=(n_grid, n_grid)) -grid_m = ti.field(dtype=float, shape=(n_grid, n_grid)) -restT = ti.Matrix.field(dim, dim, dtype=float, shape=n_elements) # noqa: N816 -total_energy = ti.field(dtype=float, shape=(), needs_grad=True) -vertices = ti.field(dtype=ti.i32, shape=(n_elements, 3)) - - -@ti.func -def mesh(i, j): # noqa: D103 - return i * n_particle_y + j - - -@ti.func -def compute_T(i): # noqa: N802, D103 - a = vertices[i, 0] - b = vertices[i, 1] - c = vertices[i, 2] - ab = x[b] - x[a] - ac = x[c] - x[a] - return ti.Matrix([[ab[0], ac[0]], [ab[1], ac[1]]]) - - -@ti.kernel -def initialize(): # noqa: D103 - for i in range(n_particle_x): - for j in range(n_particle_y): - t = mesh(i, j) - x[t] = [0.1 + i * dx * 0.5, 0.7 + j * dx * 0.5] - v[t] = [0, -1] - - # build mesh - for i in range(n_particle_x - 1): - for j in range(n_particle_y - 1): - # element id - eid = (i * (n_particle_y - 1) + j) * 2 - vertices[eid, 0] = mesh(i, j) - vertices[eid, 1] = mesh(i + 1, j) - vertices[eid, 2] = mesh(i, j + 1) - - eid = (i * (n_particle_y - 1) + j) * 2 + 1 - vertices[eid, 0] = mesh(i, j + 1) - vertices[eid, 1] = mesh(i + 1, j + 1) - vertices[eid, 2] = mesh(i + 1, j) - - for i in range(n_elements): - restT[i] = compute_T(i) # Compute rest T - - -@ti.kernel -def compute_total_energy(): # noqa: D103 - for i in range(n_elements): - currentT = compute_T(i) # noqa: N806 - F = currentT @ restT[i].inverse() # noqa: N806 - # NeoHookean - I1 = (F @ F.transpose()).trace() # noqa: N806 - J = F.determinant() # noqa: N806 - element_energy = ( - 0.5 * mu * (I1 - 2) - mu * ti.log(J) + 0.5 * la * ti.log(J) ** 2 - ) - total_energy[None] += E * element_energy * dx * dx - - -@ti.kernel -def p2g(): # noqa: D103 - for p in x: - base = ti.cast(x[p] * inv_dx - 0.5, ti.i32) - fx = x[p] * inv_dx - ti.cast(base, float) - w = [ - 0.5 * (1.5 - fx) ** 2, - 0.75 - (fx - 1) ** 2, - 0.5 * (fx - 0.5) ** 2, - ] - affine = p_mass * C[p] - for i in ti.static(range(3)): - for j in ti.static(range(3)): - I = ti.Vector([i, j]) # noqa: N806, E741 - dpos = (float(I) - fx) * dx - weight = w[i].x * w[j].y - grid_v[base + I] += weight * ( - p_mass * v[p] - dt * x.grad[p] + affine @ dpos - ) - grid_m[base + I] += weight * p_mass - - -bound = 3 - - -@ti.kernel -def grid_op(): # noqa: D103 - for i, j in grid_m: - if grid_m[i, j] > 0: - inv_m = 1 / grid_m[i, j] - grid_v[i, j] = inv_m * grid_v[i, j] - grid_v[i, j].y -= dt * 9.8 - - # center collision circle - dist = ti.Vector([i * dx - 0.5, j * dx - 0.5]) - if dist.norm_sqr() < 0.005: # noqa: PLR2004 - dist = dist.normalized() - grid_v[i, j] -= dist * ti.min(0, grid_v[i, j].dot(dist)) - - # box - if i < bound and grid_v[i, j].x < 0: - grid_v[i, j].x = 0 - if i > n_grid - bound and grid_v[i, j].x > 0: - grid_v[i, j].x = 0 - if j < bound and grid_v[i, j].y < 0: - grid_v[i, j].y = 0 - if j > n_grid - bound and grid_v[i, j].y > 0: - grid_v[i, j].y = 0 - - -@ti.kernel -def g2p(): # noqa: D103 - for p in x: - base = ti.cast(x[p] * inv_dx - 0.5, ti.i32) - fx = x[p] * inv_dx - float(base) - w = [ - 0.5 * (1.5 - fx) ** 2, - 0.75 - (fx - 1.0) ** 2, - 0.5 * (fx - 0.5) ** 2, - ] - new_v = ti.Vector([0.0, 0.0]) - new_C = ti.Matrix([[0.0, 0.0], [0.0, 0.0]]) # noqa: N806 - - for i in ti.static(range(3)): - for j in ti.static(range(3)): - I = ti.Vector([i, j]) # noqa: N806, E741 - dpos = float(I) - fx - g_v = grid_v[base + I] - weight = w[i].x * w[j].y - new_v += weight * g_v - new_C += 4 * weight * g_v.outer_product(dpos) * inv_dx # noqa: N806 - - v[p] = new_v - x[p] += dt * v[p] - C[p] = new_C - - -gui = ti.GUI('MPM', (640, 640), background_color=0x112F41) - - -def main(): # noqa: D103 - initialize() - - vertices_ = vertices.to_numpy() - - while gui.running and not gui.get_event(gui.ESCAPE): - for s in range(int(1e-2 // dt)): # noqa: B007 - grid_m.fill(0) - grid_v.fill(0) - # Note that we are now differentiating the total energy w.r.t. the particle position. - # Recall that F = - \partial (total_energy) / \partial x - with ti.ad.Tape(total_energy): - # Do the forward computation of total energy and backward propagation for x.grad, which is later used in p2g - compute_total_energy() - # It's OK not to use the computed total_energy at all, since we only need x.grad - p2g() - grid_op() - g2p() - - gui.circle((0.5, 0.5), radius=45, color=0x068587) - particle_pos = x.to_numpy() - a = vertices_.reshape(n_elements * 3) - b = np.roll(vertices_, shift=1, axis=1).reshape(n_elements * 3) - gui.lines(particle_pos[a], particle_pos[b], radius=1, color=0x4FB99F) - gui.circles(particle_pos, radius=1.5, color=0xF2B134) - gui.line( - (0.00, 0.03 / quality), - (1.0, 0.03 / quality), - color=0xFFFFFF, - radius=3, - ) - gui.show() - - -if __name__ == '__main__': - main() +import numpy as np + +import taichi as ti + +ti.init(arch=ti.gpu) + +dim = 2 +quality = 1 # Use a larger integral number for higher quality +n_particle_x = 100 * quality +n_particle_y = 8 * quality +n_particles = n_particle_x * n_particle_y +n_elements = (n_particle_x - 1) * (n_particle_y - 1) * 2 +n_grid = 64 * quality +dx = 1 / n_grid +inv_dx = 1 / dx +dt = 1e-4 / quality +E = 25000 +p_mass = 1 +p_vol = 1 +mu = 1 +la = 1 + +x = ti.Vector.field(dim, dtype=float, shape=n_particles, needs_grad=True) +v = ti.Vector.field(dim, dtype=float, shape=n_particles) +C = ti.Matrix.field(dim, dim, dtype=float, shape=n_particles) +grid_v = ti.Vector.field(dim, dtype=float, shape=(n_grid, n_grid)) +grid_m = ti.field(dtype=float, shape=(n_grid, n_grid)) +restT = ti.Matrix.field(dim, dim, dtype=float, shape=n_elements) +total_energy = ti.field(dtype=float, shape=(), needs_grad=True) +vertices = ti.field(dtype=ti.i32, shape=(n_elements, 3)) + + +@ti.func +def mesh(i, j): + return i * n_particle_y + j + + +@ti.func +def compute_T(i): + a = vertices[i, 0] + b = vertices[i, 1] + c = vertices[i, 2] + ab = x[b] - x[a] + ac = x[c] - x[a] + return ti.Matrix([[ab[0], ac[0]], [ab[1], ac[1]]]) + + +@ti.kernel +def initialize(): + for i in range(n_particle_x): + for j in range(n_particle_y): + t = mesh(i, j) + x[t] = [0.1 + i * dx * 0.5, 0.7 + j * dx * 0.5] + v[t] = [0, -1] + + # build mesh + for i in range(n_particle_x - 1): + for j in range(n_particle_y - 1): + # element id + eid = (i * (n_particle_y - 1) + j) * 2 + vertices[eid, 0] = mesh(i, j) + vertices[eid, 1] = mesh(i + 1, j) + vertices[eid, 2] = mesh(i, j + 1) + + eid = (i * (n_particle_y - 1) + j) * 2 + 1 + vertices[eid, 0] = mesh(i, j + 1) + vertices[eid, 1] = mesh(i + 1, j + 1) + vertices[eid, 2] = mesh(i + 1, j) + + for i in range(n_elements): + restT[i] = compute_T(i) # Compute rest T + + +@ti.kernel +def compute_total_energy(): + for i in range(n_elements): + currentT = compute_T(i) + F = currentT @ restT[i].inverse() + # NeoHookean + I1 = (F @ F.transpose()).trace() + J = F.determinant() + element_energy = 0.5 * mu * (I1 - 2) - mu * ti.log(J) + 0.5 * la * ti.log(J) ** 2 + total_energy[None] += E * element_energy * dx * dx + + +@ti.kernel +def p2g(): + for p in x: + base = ti.cast(x[p] * inv_dx - 0.5, ti.i32) + fx = x[p] * inv_dx - ti.cast(base, float) + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2] + affine = p_mass * C[p] + for i in ti.static(range(3)): + for j in ti.static(range(3)): + I = ti.Vector([i, j]) + dpos = (float(I) - fx) * dx + weight = w[i].x * w[j].y + grid_v[base + I] += weight * (p_mass * v[p] - dt * x.grad[p] + affine @ dpos) + grid_m[base + I] += weight * p_mass + + +bound = 3 + + +@ti.kernel +def grid_op(): + for i, j in grid_m: + if grid_m[i, j] > 0: + inv_m = 1 / grid_m[i, j] + grid_v[i, j] = inv_m * grid_v[i, j] + grid_v[i, j].y -= dt * 9.8 + + # center collision circle + dist = ti.Vector([i * dx - 0.5, j * dx - 0.5]) + if dist.norm_sqr() < 0.005: + dist = dist.normalized() + grid_v[i, j] -= dist * ti.min(0, grid_v[i, j].dot(dist)) + + # box + if i < bound and grid_v[i, j].x < 0: + grid_v[i, j].x = 0 + if i > n_grid - bound and grid_v[i, j].x > 0: + grid_v[i, j].x = 0 + if j < bound and grid_v[i, j].y < 0: + grid_v[i, j].y = 0 + if j > n_grid - bound and grid_v[i, j].y > 0: + grid_v[i, j].y = 0 + + +@ti.kernel +def g2p(): + for p in x: + base = ti.cast(x[p] * inv_dx - 0.5, ti.i32) + fx = x[p] * inv_dx - float(base) + w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1.0) ** 2, 0.5 * (fx - 0.5) ** 2] + new_v = ti.Vector([0.0, 0.0]) + new_C = ti.Matrix([[0.0, 0.0], [0.0, 0.0]]) + + for i in ti.static(range(3)): + for j in ti.static(range(3)): + I = ti.Vector([i, j]) + dpos = float(I) - fx + g_v = grid_v[base + I] + weight = w[i].x * w[j].y + new_v += weight * g_v + new_C += 4 * weight * g_v.outer_product(dpos) * inv_dx + + v[p] = new_v + x[p] += dt * v[p] + C[p] = new_C + + +gui = ti.GUI("MPM", (640, 640), background_color=0x112F41) + + +def main(): + initialize() + + vertices_ = vertices.to_numpy() + + while gui.running and not gui.get_event(gui.ESCAPE): + for s in range(int(1e-2 // dt)): + grid_m.fill(0) + grid_v.fill(0) + # Note that we are now differentiating the total energy w.r.t. the particle position. + # Recall that F = - \partial (total_energy) / \partial x + with ti.ad.Tape(total_energy): + # Do the forward computation of total energy and backward propagation for x.grad, which is later used in p2g + compute_total_energy() + # It's OK not to use the computed total_energy at all, since we only need x.grad + p2g() + grid_op() + g2p() + + gui.circle((0.5, 0.5), radius=45, color=0x068587) + particle_pos = x.to_numpy() + a = vertices_.reshape(n_elements * 3) + b = np.roll(vertices_, shift=1, axis=1).reshape(n_elements * 3) + gui.lines(particle_pos[a], particle_pos[b], radius=1, color=0x4FB99F) + gui.circles(particle_pos, radius=1.5, color=0xF2B134) + gui.line((0.00, 0.03 / quality), (1.0, 0.03 / quality), color=0xFFFFFF, radius=3) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/nbody.py b/modules/createEVENT/TaichiEvent/nbody.py new file mode 100644 index 000000000..68c487314 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/nbody.py @@ -0,0 +1,105 @@ +# Authored by Tiantian Liu, Taichi Graphics. +import math + +import taichi as ti + +ti.init(arch=ti.cuda) + +# global control +paused = ti.field(ti.i32, ()) + +# gravitational constant 6.67408e-11, using 1 for simplicity +G = 1 + +# number of planets +N = 3000 +# unit mass +m = 1 +# galaxy size +galaxy_size = 0.4 +# planet radius (for rendering) +planet_radius = 2 +# init vel +init_vel = 120 + +# time-step size +h = 1e-4 +# substepping +substepping = 10 + +# center of the screen +center = ti.Vector.field(2, ti.f32, ()) + +# pos, vel and force of the planets +# Nx2 vectors +pos = ti.Vector.field(2, ti.f32, N) +vel = ti.Vector.field(2, ti.f32, N) +force = ti.Vector.field(2, ti.f32, N) + + +@ti.kernel +def initialize(): + center[None] = [0.5, 0.5] + for i in range(N): + theta = ti.random() * 2 * math.pi + r = (ti.sqrt(ti.random()) * 0.6 + 0.4) * galaxy_size + offset = r * ti.Vector([ti.cos(theta), ti.sin(theta)]) + pos[i] = center[None] + offset + vel[i] = [-offset.y, offset.x] + vel[i] *= init_vel + + +@ti.kernel +def compute_force(): + # clear force + for i in range(N): + force[i] = [0.0, 0.0] + + # compute gravitational force + for i in range(N): + p = pos[i] + for j in range(N): + if i != j: # double the computation for a better memory footprint and load balance + diff = p - pos[j] + r = diff.norm(1e-5) + + # gravitational force -(GMm / r^2) * (diff/r) for i + f = -G * m * m * (1.0 / r) ** 3 * diff + + # assign to each particle + force[i] += f + + +@ti.kernel +def update(): + dt = h / substepping + for i in range(N): + # symplectic euler + vel[i] += dt * force[i] / m + pos[i] += dt * vel[i] + + +def main(): + gui = ti.GUI("N-body problem", (800, 800)) + + initialize() + while gui.running: + for e in gui.get_events(ti.GUI.PRESS): + if e.key in [ti.GUI.ESCAPE, ti.GUI.EXIT]: + exit() + elif e.key == "r": + initialize() + elif e.key == ti.GUI.SPACE: + paused[None] = not paused[None] + + if not paused[None]: + for i in range(substepping): + compute_force() + update() + + gui.circles(pos.to_numpy(), color=0xFFFFFF, radius=planet_radius) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/odop_solar.py b/modules/createEVENT/TaichiEvent/odop_solar.py new file mode 100644 index 000000000..7c62247e8 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/odop_solar.py @@ -0,0 +1,71 @@ +import math + +import taichi as ti + +ti.init() + + +@ti.data_oriented +class SolarSystem: + def __init__(self, n, dt): # Initializer of the solar system simulator + self.n = n + self.dt = dt + self.x = ti.Vector.field(2, dtype=ti.f32, shape=n) + self.v = ti.Vector.field(2, dtype=ti.f32, shape=n) + self.center = ti.Vector.field(2, dtype=ti.f32, shape=()) + + @staticmethod + @ti.func + def random_vector(radius): # Create a random vector in circle + theta = ti.random() * 2 * math.pi + r = ti.random() * radius + return r * ti.Vector([ti.cos(theta), ti.sin(theta)]) + + @ti.kernel + def initialize_particles(self): + # (Re)initialize particle position/velocities + for i in range(self.n): + offset = self.random_vector(0.5) + self.x[i] = self.center[None] + offset # Offset from center + self.v[i] = [-offset.y, offset.x] # Perpendicular to offset + self.v[i] += self.random_vector(0.02) # Random velocity noise + self.v[i] *= 1 / offset.norm() ** 1.5 # Kepler's third law + + @ti.func + def gravity(self, pos): # Compute gravity at pos + offset = -(pos - self.center[None]) + return offset / offset.norm() ** 3 + + @ti.kernel + def integrate(self): # Semi-implicit Euler time integration + for i in range(self.n): + self.v[i] += self.dt * self.gravity(self.x[i]) + self.x[i] += self.dt * self.v[i] + + @staticmethod + def render(gui): # Render the scene on GUI + gui.circle([0.5, 0.5], radius=10, color=0xFFAA88) + gui.circles(solar.x.to_numpy(), radius=3, color=0xFFFFFF) + + +def main(): + global solar + + solar = SolarSystem(8, 0.0001) + solar.center[None] = [0.5, 0.5] + solar.initialize_particles() + + gui = ti.GUI("Solar System", background_color=0x0071A) + while gui.running: + if gui.get_event() and gui.is_pressed(gui.SPACE): + solar.initialize_particles() # reinitialize when space bar pressed. + + for _ in range(10): # Time integration + solar.integrate() + + solar.render(gui) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/pbf2d.py b/modules/createEVENT/TaichiEvent/pbf2d.py new file mode 100644 index 000000000..08b237807 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/pbf2d.py @@ -0,0 +1,306 @@ +# Macklin, M. and Müller, M., 2013. Position based fluids. ACM Transactions on Graphics (TOG), 32(4), p.104. +# Taichi implementation by Ye Kuang (k-ye) + +import math + +import numpy as np + +import taichi as ti + +ti.init(arch=ti.gpu) + +screen_res = (800, 400) +screen_to_world_ratio = 10.0 +boundary = ( + screen_res[0] / screen_to_world_ratio, + screen_res[1] / screen_to_world_ratio, +) +cell_size = 2.51 +cell_recpr = 1.0 / cell_size + + +def round_up(f, s): + return (math.floor(f * cell_recpr / s) + 1) * s + + +grid_size = (round_up(boundary[0], 1), round_up(boundary[1], 1)) + +dim = 2 +bg_color = 0x112F41 +particle_color = 0x068587 +boundary_color = 0xEBACA2 +num_particles_x = 60 +num_particles = num_particles_x * 20 +max_num_particles_per_cell = 100 +max_num_neighbors = 100 +time_delta = 1.0 / 20.0 +epsilon = 1e-5 +particle_radius = 3.0 +particle_radius_in_world = particle_radius / screen_to_world_ratio + +# PBF params +h_ = 1.1 +mass = 1.0 +rho0 = 1.0 +lambda_epsilon = 100.0 +pbf_num_iters = 5 +corr_deltaQ_coeff = 0.3 +corrK = 0.001 +# Need ti.pow() +# corrN = 4.0 +neighbor_radius = h_ * 1.05 + +poly6_factor = 315.0 / 64.0 / math.pi +spiky_grad_factor = -45.0 / math.pi + +old_positions = ti.Vector.field(dim, float) +positions = ti.Vector.field(dim, float) +velocities = ti.Vector.field(dim, float) +grid_num_particles = ti.field(int) +grid2particles = ti.field(int) +particle_num_neighbors = ti.field(int) +particle_neighbors = ti.field(int) +lambdas = ti.field(float) +position_deltas = ti.Vector.field(dim, float) +# 0: x-pos, 1: timestep in sin() +board_states = ti.Vector.field(2, float) + +ti.root.dense(ti.i, num_particles).place(old_positions, positions, velocities) +grid_snode = ti.root.dense(ti.ij, grid_size) +grid_snode.place(grid_num_particles) +grid_snode.dense(ti.k, max_num_particles_per_cell).place(grid2particles) +nb_node = ti.root.dense(ti.i, num_particles) +nb_node.place(particle_num_neighbors) +nb_node.dense(ti.j, max_num_neighbors).place(particle_neighbors) +ti.root.dense(ti.i, num_particles).place(lambdas, position_deltas) +ti.root.place(board_states) + + +@ti.func +def poly6_value(s, h): + result = 0.0 + if 0 < s and s < h: + x = (h * h - s * s) / (h * h * h) + result = poly6_factor * x * x * x + return result + + +@ti.func +def spiky_gradient(r, h): + result = ti.Vector([0.0, 0.0]) + r_len = r.norm() + if 0 < r_len and r_len < h: + x = (h - r_len) / (h * h * h) + g_factor = spiky_grad_factor * x * x + result = r * g_factor / r_len + return result + + +@ti.func +def compute_scorr(pos_ji): + # Eq (13) + x = poly6_value(pos_ji.norm(), h_) / poly6_value(corr_deltaQ_coeff * h_, h_) + # pow(x, 4) + x = x * x + x = x * x + return (-corrK) * x + + +@ti.func +def get_cell(pos): + return int(pos * cell_recpr) + + +@ti.func +def is_in_grid(c): + # @c: Vector(i32) + return 0 <= c[0] and c[0] < grid_size[0] and 0 <= c[1] and c[1] < grid_size[1] + + +@ti.func +def confine_position_to_boundary(p): + bmin = particle_radius_in_world + bmax = ti.Vector([board_states[None][0], boundary[1]]) - particle_radius_in_world + for i in ti.static(range(dim)): + # Use randomness to prevent particles from sticking into each other after clamping + if p[i] <= bmin: + p[i] = bmin + epsilon * ti.random() + elif bmax[i] <= p[i]: + p[i] = bmax[i] - epsilon * ti.random() + return p + + +@ti.kernel +def move_board(): + # probably more accurate to exert force on particles according to hooke's law. + b = board_states[None] + b[1] += 1.0 + period = 90 + vel_strength = 8.0 + if b[1] >= 2 * period: + b[1] = 0 + b[0] += -ti.sin(b[1] * np.pi / period) * vel_strength * time_delta + board_states[None] = b + + +@ti.kernel +def prologue(): + # save old positions + for i in positions: + old_positions[i] = positions[i] + # apply gravity within boundary + for i in positions: + g = ti.Vector([0.0, -9.8]) + pos, vel = positions[i], velocities[i] + vel += g * time_delta + pos += vel * time_delta + positions[i] = confine_position_to_boundary(pos) + + # clear neighbor lookup table + for I in ti.grouped(grid_num_particles): + grid_num_particles[I] = 0 + for I in ti.grouped(particle_neighbors): + particle_neighbors[I] = -1 + + # update grid + for p_i in positions: + cell = get_cell(positions[p_i]) + # ti.Vector doesn't seem to support unpacking yet + # but we can directly use int Vectors as indices + offs = ti.atomic_add(grid_num_particles[cell], 1) + grid2particles[cell, offs] = p_i + # find particle neighbors + for p_i in positions: + pos_i = positions[p_i] + cell = get_cell(pos_i) + nb_i = 0 + for offs in ti.static(ti.grouped(ti.ndrange((-1, 2), (-1, 2)))): + cell_to_check = cell + offs + if is_in_grid(cell_to_check): + for j in range(grid_num_particles[cell_to_check]): + p_j = grid2particles[cell_to_check, j] + if nb_i < max_num_neighbors and p_j != p_i and (pos_i - positions[p_j]).norm() < neighbor_radius: + particle_neighbors[p_i, nb_i] = p_j + nb_i += 1 + particle_num_neighbors[p_i] = nb_i + + +@ti.kernel +def substep(): + # compute lambdas + # Eq (8) ~ (11) + for p_i in positions: + pos_i = positions[p_i] + + grad_i = ti.Vector([0.0, 0.0]) + sum_gradient_sqr = 0.0 + density_constraint = 0.0 + + for j in range(particle_num_neighbors[p_i]): + p_j = particle_neighbors[p_i, j] + if p_j < 0: + break + pos_ji = pos_i - positions[p_j] + grad_j = spiky_gradient(pos_ji, h_) + grad_i += grad_j + sum_gradient_sqr += grad_j.dot(grad_j) + # Eq(2) + density_constraint += poly6_value(pos_ji.norm(), h_) + + # Eq(1) + density_constraint = (mass * density_constraint / rho0) - 1.0 + + sum_gradient_sqr += grad_i.dot(grad_i) + lambdas[p_i] = (-density_constraint) / (sum_gradient_sqr + lambda_epsilon) + # compute position deltas + # Eq(12), (14) + for p_i in positions: + pos_i = positions[p_i] + lambda_i = lambdas[p_i] + + pos_delta_i = ti.Vector([0.0, 0.0]) + for j in range(particle_num_neighbors[p_i]): + p_j = particle_neighbors[p_i, j] + if p_j < 0: + break + lambda_j = lambdas[p_j] + pos_ji = pos_i - positions[p_j] + scorr_ij = compute_scorr(pos_ji) + pos_delta_i += (lambda_i + lambda_j + scorr_ij) * spiky_gradient(pos_ji, h_) + + pos_delta_i /= rho0 + position_deltas[p_i] = pos_delta_i + # apply position deltas + for i in positions: + positions[i] += position_deltas[i] + + +@ti.kernel +def epilogue(): + # confine to boundary + for i in positions: + pos = positions[i] + positions[i] = confine_position_to_boundary(pos) + # update velocities + for i in positions: + velocities[i] = (positions[i] - old_positions[i]) / time_delta + # no vorticity/xsph because we cannot do cross product in 2D... + + +def run_pbf(): + prologue() + for _ in range(pbf_num_iters): + substep() + epilogue() + + +def render(gui): + gui.clear(bg_color) + pos_np = positions.to_numpy() + for j in range(dim): + pos_np[:, j] *= screen_to_world_ratio / screen_res[j] + gui.circles(pos_np, radius=particle_radius, color=particle_color) + gui.rect( + (0, 0), + (board_states[None][0] / boundary[0], 1), + radius=1.5, + color=boundary_color, + ) + gui.show() + + +@ti.kernel +def init_particles(): + for i in range(num_particles): + delta = h_ * 0.8 + offs = ti.Vector([(boundary[0] - delta * num_particles_x) * 0.5, boundary[1] * 0.02]) + positions[i] = ti.Vector([i % num_particles_x, i // num_particles_x]) * delta + offs + for c in ti.static(range(dim)): + velocities[i][c] = (ti.random() - 0.5) * 4 + board_states[None] = ti.Vector([boundary[0] - epsilon, -0.0]) + + +def print_stats(): + print("PBF stats:") + num = grid_num_particles.to_numpy() + avg, max_ = np.mean(num), np.max(num) + print(f" #particles per cell: avg={avg:.2f} max={max_}") + num = particle_num_neighbors.to_numpy() + avg, max_ = np.mean(num), np.max(num) + print(f" #neighbors per particle: avg={avg:.2f} max={max_}") + + +def main(): + init_particles() + print(f"boundary={boundary} grid={grid_size} cell_size={cell_size}") + gui = ti.GUI("PBF2D", screen_res) + while gui.running and not gui.get_event(gui.ESCAPE): + move_board() + run_pbf() + if gui.frame % 20 == 1: + print_stats() + render(gui) + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/physarum.py b/modules/createEVENT/TaichiEvent/physarum.py new file mode 100644 index 000000000..701cf56b9 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/physarum.py @@ -0,0 +1,86 @@ +"""Physarum simulation example. + +See https://sagejenson.com/physarum for the details.""" + +import numpy as np + +import taichi as ti + +ti.init(arch=ti.gpu) + +PARTICLE_N = 1024 +GRID_SIZE = 512 +SENSE_ANGLE = 0.20 * np.pi +SENSE_DIST = 4.0 +EVAPORATION = 0.95 +MOVE_ANGLE = 0.1 * np.pi +MOVE_STEP = 2.0 + +grid = ti.field(dtype=ti.f32, shape=[2, GRID_SIZE, GRID_SIZE]) +position = ti.Vector.field(2, dtype=ti.f32, shape=[PARTICLE_N]) +heading = ti.field(dtype=ti.f32, shape=[PARTICLE_N]) + + +@ti.kernel +def init(): + for p in ti.grouped(grid): + grid[p] = 0.0 + for i in position: + position[i] = ti.Vector([ti.random(), ti.random()]) * GRID_SIZE + heading[i] = ti.random() * np.pi * 2.0 + + +@ti.func +def sense(phase, pos, ang): + p = pos + ti.Vector([ti.cos(ang), ti.sin(ang)]) * SENSE_DIST + return grid[phase, p.cast(int) % GRID_SIZE] + + +@ti.kernel +def step(phase: ti.i32): + # move + for i in position: + pos, ang = position[i], heading[i] + l = sense(phase, pos, ang - SENSE_ANGLE) + c = sense(phase, pos, ang) + r = sense(phase, pos, ang + SENSE_ANGLE) + if l < c < r: + ang += MOVE_ANGLE + elif l > c > r: + ang -= MOVE_ANGLE + elif c < l and c < r: + ang += MOVE_ANGLE * (2 * (ti.random() < 0.5) - 1) + pos += ti.Vector([ti.cos(ang), ti.sin(ang)]) * MOVE_STEP + position[i], heading[i] = pos, ang + + # deposit + for i in position: + ipos = position[i].cast(int) % GRID_SIZE + grid[phase, ipos] += 1.0 + + # diffuse + for i, j in ti.ndrange(GRID_SIZE, GRID_SIZE): + a = 0.0 + for di in ti.static(range(-1, 2)): + for dj in ti.static(range(-1, 2)): + a += grid[phase, (i + di) % GRID_SIZE, (j + dj) % GRID_SIZE] + a *= EVAPORATION / 9.0 + grid[1 - phase, i, j] = a + + +def main(): + print("[Hint] Use slider to change simulation speed.") + gui = ti.GUI("Physarum") + init() + i = 0 + step_per_frame = gui.slider("step_per_frame", 1, 100, 1) + while gui.running and not gui.get_event(gui.ESCAPE): + for _ in range(int(step_per_frame.value)): + step(i % 2) + i += 1 + gui.set_image(grid.to_numpy()[0]) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/snow_phaseField.py b/modules/createEVENT/TaichiEvent/snow_phaseField.py new file mode 100644 index 000000000..a25479121 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/snow_phaseField.py @@ -0,0 +1,235 @@ +""" phase field model for snow dendrite evolution, which simulates the snow growing from water vapor. + space discretization: finite difference method, time integration: Runge-Kutta method + repo's link: https://github.com/mo-hanxuan/Snow-PhaseField + more details about physical interpretation refer to [Physica D 63(3-4): 410-423]""" +import numpy as np + +import taichi as ti + + +@ti.data_oriented +class Dendrite: + def __init__( + self, + dx=0.03, # grid space + dt=2.0e-4, # time step + n=512, # field shape + dtype=ti.f64, # data type + n_fold_symmetry=6, # dendrites number + angle0=0.0, # initial angle + ): + self.n = n # the size of the field + ### phase field indicate solid phase (phi=1) and water vapor phase (phi=0) + self.phi = ti.field(dtype=dtype, shape=(n, n)) + self.temperature = ti.field(dtype=dtype, shape=(n, n)) + self.phi_old = ti.field(dtype=dtype, shape=(n, n)) + self.temperature_old = ti.field(dtype=dtype, shape=(n, n)) + self.dEnergy_dGrad_term1 = ti.Vector.field(2, dtype, shape=(n, n)) + self.epsilons = ti.field(dtype=dtype, shape=(n, n)) # anisotropic gradient energy coefficient + self.phiRate = ti.Vector.field(4, dtype=dtype, shape=(n, n)) # rate of phi, with RK4 method + self.temperatureRate = ti.Vector.field(4, dtype=dtype, shape=(n, n)) + + self.dx = dx + self.dt = dt # maximum dt before going unstable is 2.9e-4 + self.mobility = 128.0e2 + + self.grad_energy_coef = 0.006 # magnitude of gradient energy coefficient + self.latent_heat_coef = 1.5 # latent heat coefficient + self.aniso_magnitude = 0.12 # the magnitude of anisotropy + self.n_fold_symmetry = n_fold_symmetry # n-fold rotational symmetry of the crystalline structure + self.angle0 = angle0 / 180.0 * np.pi # initial tilt angle of the crystal + + ### m(T) = (α / π) * arctan[γ(T_equi - T)], m determines the relative potential between water vapor and crystal + self.alpha = 0.9 / np.pi # α + self.gamma = 10.0 # γ + self.temperature_equi = 1.0 # temperature of equilibrium state + + ### parameters for RK4 method + self.dtRatio_rk4 = ti.field(ti.f64, shape=(4)) + self.dtRatio_rk4.from_numpy(np.array([0.0, 0.5, 0.5, 1.0])) + self.weights_rk4 = ti.field(ti.f64, shape=(4)) + self.weights_rk4.from_numpy(np.array([1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0])) + + self.showFrameFrequency = int(4 * 1.0e-4 / self.dt) + + @ti.kernel + def initialize( + self, + ): + radius = 4.0 # 1. + center = ti.Vector([self.n // 2, self.n // 2]) + for I in ti.grouped(self.phi): + if ((ti.Vector([I[0], I[1]]) - center) ** 2).sum() < radius**2: + self.phi[I] = 1.0 + else: + self.phi[I] = 0.0 + self.phi_old[I] = self.phi[I] + + self.temperature[I] = 0.0 # temperature + self.temperature_old[I] = self.temperature[I] + + @ti.func + def neighbor_index(self, i, j): + """use periodic boundary condition to get neighbor index""" + im = i - 1 if i - 1 >= 0 else self.n - 1 + jm = j - 1 if j - 1 >= 0 else self.n - 1 + ip = i + 1 if i + 1 < self.n else 0 + jp = j + 1 if j + 1 < self.n else 0 + return im, jm, ip, jp + + @ti.kernel + def rk4_intermediate_update(self, rk_loop: int): + """update field variables at the intermediate step of RK4""" + phi, phi_old, temperature, temperature_old, dt = ti.static( + self.phi, self.phi_old, self.temperature, self.temperature_old, self.dt + ) + for I in ti.grouped(phi): + phi[I] = phi_old[I] + self.dtRatio_rk4[rk_loop] * dt * self.phiRate[I][rk_loop - 1] + temperature[I] = temperature_old[I] + self.dtRatio_rk4[rk_loop] * dt * self.temperatureRate[I][rk_loop - 1] + + @ti.kernel + def get_rate(self, rk_loop: int): + """get rate of phi and temperature""" + ( + phi, + temperature, + dx, + mobility, + aniso_magnitude, + epsilons, + n_fold_symmetry, + grad_energy_coef, + angle0, + ) = ti.static( + self.phi, + self.temperature, + self.dx, + self.mobility, + self.aniso_magnitude, + self.epsilons, + self.n_fold_symmetry, + self.grad_energy_coef, + self.angle0, + ) + ### first, get epsilons and dEnergy_dGrad_term1 + for i, j in phi: + im, jm, ip, jp = self.neighbor_index(i, j) + grad = ti.Vector( + [ + (phi[ip, j] - phi[im, j]) / (2.0 * dx), + (phi[i, jp] - phi[i, jm]) / (2.0 * dx), + ] + ) + gradNorm = (grad**2).sum() + if gradNorm < 1.0e-8: + self.dEnergy_dGrad_term1[i, j] = ti.Vector([0.0, 0.0]) + angle = ti.atan2(grad[1], grad[0]) + epsilons[i, j] = grad_energy_coef * (1.0 + aniso_magnitude * ti.cos(n_fold_symmetry * (angle - angle0))) + else: + angle = ti.atan2(grad[1], grad[0]) + epsilon = grad_energy_coef * (1.0 + aniso_magnitude * ti.cos(n_fold_symmetry * (angle - angle0))) + epsilons[i, j] = epsilon + dAngle_dGradX = -grad[1] / gradNorm + dAngle_dGradY = grad[0] / gradNorm + tmp = grad_energy_coef * aniso_magnitude * -ti.sin(n_fold_symmetry * (angle - angle0)) * n_fold_symmetry + depsilon_dGrad = tmp * ti.Vector([dAngle_dGradX, dAngle_dGradY]) + self.dEnergy_dGrad_term1[i, j] = epsilon * depsilon_dGrad * gradNorm + + ### then, get the phi rate and temperature rate + for i, j in phi: + im, jm, ip, jp = self.neighbor_index(i, j) + + lapla_phi = ( # laplacian of phi + 2 * (phi[im, j] + phi[i, jm] + phi[ip, j] + phi[i, jp]) + + (phi[im, jm] + phi[im, jp] + phi[ip, jm] + phi[ip, jp]) + - 12 * phi[i, j] + ) / (3.0 * dx * dx) + lapla_tp = ( # laplacian of temperature + 2 * (temperature[im, j] + temperature[i, jm] + temperature[ip, j] + temperature[i, jp]) + + (temperature[im, jm] + temperature[im, jp] + temperature[ip, jm] + temperature[ip, jp]) + - 12 * temperature[i, j] + ) / (3.0 * dx * dx) + + m_chem = self.alpha * ti.atan2(self.gamma * (self.temperature_equi - temperature[i, j]), 1.0) + chemicalForce = phi[i, j] * (1.0 - phi[i, j]) * (phi[i, j] - 0.5 + m_chem) + gradForce_term1 = self.divergence_dEnergy_dGrad_term1(i, j) + grad_epsilon2 = ti.Vector( + [ + (epsilons[ip, j] ** 2 - epsilons[im, j] ** 2) / (2.0 * dx), + (epsilons[i, jp] ** 2 - epsilons[i, jm] ** 2) / (2.0 * dx), + ] + ) + grad_phi = ti.Vector( + [ + (phi[ip, j] - phi[im, j]) / (2.0 * dx), + (phi[i, jp] - phi[i, jm]) / (2.0 * dx), + ] + ) + gradForce_term2 = ( + grad_epsilon2[0] * grad_phi[0] + grad_epsilon2[1] * grad_phi[1] + epsilons[i, j] ** 2 * lapla_phi + ) + + self.phiRate[i, j][rk_loop] = mobility * (chemicalForce + gradForce_term1 + gradForce_term2) + self.temperatureRate[i, j][rk_loop] = lapla_tp + self.latent_heat_coef * self.phiRate[i, j][rk_loop] + + @ti.kernel + def rk4_total_update( + self, + ): + """the final step in RK4 (Runge-Kutta) process""" + ( + dt, + phi, + phi_old, + temperature, + temperature_old, + phiRate, + temperatureRate, + ) = ti.static( + self.dt, + self.phi, + self.phi_old, + self.temperature, + self.temperature_old, + self.phiRate, + self.temperatureRate, + ) + for I in ti.grouped(phi): + for k in ti.static(range(4)): + phi_old[I] = phi_old[I] + self.weights_rk4[k] * dt * phiRate[I][k] + temperature_old[I] = temperature_old[I] + self.weights_rk4[k] * dt * temperatureRate[I][k] + for I in ti.grouped(phi): + phi[I] = phi_old[I] + temperature[I] = temperature_old[I] + + @ti.func + def divergence_dEnergy_dGrad_term1(self, i, j): + im, jm, ip, jp = self.neighbor_index(i, j) + return (self.dEnergy_dGrad_term1[ip, j][0] - self.dEnergy_dGrad_term1[im, j][0]) / (2.0 * self.dx) + ( + self.dEnergy_dGrad_term1[i, jp][1] - self.dEnergy_dGrad_term1[i, jm][1] + ) / (2.0 * self.dx) + + def advance( + self, + ): # advance a time step + self.get_rate(rk_loop=0) + for rk_loop in range(1, 4): + self.rk4_intermediate_update(rk_loop=rk_loop) + self.get_rate(rk_loop=rk_loop) + self.rk4_total_update() + + def getDendritic(self, steps=2048): + self.initialize() + gui = ti.GUI("phase field", res=(self.n, self.n)) + while not gui.get_event(ti.GUI.ESCAPE, ti.GUI.EXIT): + for _ in range(self.showFrameFrequency): + self.advance() + + gui.set_image(self.phi) + gui.show() + return self.phi + + +if __name__ == "__main__": + ti.init(arch=ti.cuda, default_fp=ti.f64) + Dendrite().getDendritic(steps=10000) diff --git a/modules/createEVENT/TaichiEvent/stable_fluid.py b/modules/createEVENT/TaichiEvent/stable_fluid.py new file mode 100644 index 000000000..316d33218 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/stable_fluid.py @@ -0,0 +1,395 @@ +# References: +# http://developer.download.nvidia.com/books/HTML/gpugems/gpugems_ch38.html +# https://github.com/PavelDoGreat/WebGL-Fluid-Simulation +# https://www.bilibili.com/video/BV1ZK411H7Hc?p=4 +# https://github.com/ShaneFX/GAMES201/tree/master/HW01 + +import argparse + +import numpy as np + +import taichi as ti + +# How to run: +# `python stable_fluid.py`: use the jacobi iteration to solve the linear system. +# `python stable_fluid.py -S`: use a sparse matrix to do so. +parser = argparse.ArgumentParser() +parser.add_argument( + "-S", + "--use-sp-mat", + action="store_true", + help="Solve Poisson's equation by using a sparse matrix", +) +parser.add_argument( + "-a", + "--arch", + required=False, + default="cpu", + dest="arch", + type=str, + help="The arch (backend) to run this example on", +) +args, unknowns = parser.parse_known_args() + +res = 512 +dt = 0.03 +p_jacobi_iters = 500 # 40 for a quicker but less accurate result +f_strength = 10000.0 +curl_strength = 0 +time_c = 2 +maxfps = 60 +dye_decay = 1 - 1 / (maxfps * time_c) +force_radius = res / 2.0 +debug = False + +use_sparse_matrix = args.use_sp_mat +arch = args.arch +if arch in ["x64", "cpu", "arm64"]: + ti.init(arch=ti.cpu) +elif arch in ["cuda", "gpu"]: + ti.init(arch=ti.cuda) +else: + raise ValueError("Only CPU and CUDA backends are supported for now.") + +if use_sparse_matrix: + print("Using sparse matrix") +else: + print("Using jacobi iteration") + +_velocities = ti.Vector.field(2, float, shape=(res, res)) +_new_velocities = ti.Vector.field(2, float, shape=(res, res)) +velocity_divs = ti.field(float, shape=(res, res)) +velocity_curls = ti.field(float, shape=(res, res)) +_pressures = ti.field(float, shape=(res, res)) +_new_pressures = ti.field(float, shape=(res, res)) +_dye_buffer = ti.Vector.field(3, float, shape=(res, res)) +_new_dye_buffer = ti.Vector.field(3, float, shape=(res, res)) + + +class TexPair: + def __init__(self, cur, nxt): + self.cur = cur + self.nxt = nxt + + def swap(self): + self.cur, self.nxt = self.nxt, self.cur + + +velocities_pair = TexPair(_velocities, _new_velocities) +pressures_pair = TexPair(_pressures, _new_pressures) +dyes_pair = TexPair(_dye_buffer, _new_dye_buffer) + +if use_sparse_matrix: + # use a sparse matrix to solve Poisson's pressure equation. + @ti.kernel + def fill_laplacian_matrix(A: ti.types.sparse_matrix_builder()): + for i, j in ti.ndrange(res, res): + row = i * res + j + center = 0.0 + if j != 0: + A[row, row - 1] += -1.0 + center += 1.0 + if j != res - 1: + A[row, row + 1] += -1.0 + center += 1.0 + if i != 0: + A[row, row - res] += -1.0 + center += 1.0 + if i != res - 1: + A[row, row + res] += -1.0 + center += 1.0 + A[row, row] += center + + N = res * res + K = ti.linalg.SparseMatrixBuilder(N, N, max_num_triplets=N * 6) + F_b = ti.ndarray(ti.f32, shape=N) + + fill_laplacian_matrix(K) + L = K.build() + solver = ti.linalg.SparseSolver(solver_type="LLT") + solver.analyze_pattern(L) + solver.factorize(L) + + +@ti.func +def sample(qf, u, v): + I = ti.Vector([int(u), int(v)]) + I = ti.max(0, ti.min(res - 1, I)) + return qf[I] + + +@ti.func +def lerp(vl, vr, frac): + # frac: [0.0, 1.0] + return vl + frac * (vr - vl) + + +@ti.func +def bilerp(vf, p): + u, v = p + s, t = u - 0.5, v - 0.5 + # floor + iu, iv = ti.floor(s), ti.floor(t) + # fract + fu, fv = s - iu, t - iv + a = sample(vf, iu, iv) + b = sample(vf, iu + 1, iv) + c = sample(vf, iu, iv + 1) + d = sample(vf, iu + 1, iv + 1) + return lerp(lerp(a, b, fu), lerp(c, d, fu), fv) + + +# 3rd order Runge-Kutta +@ti.func +def backtrace(vf: ti.template(), p, dt_: ti.template()): + v1 = bilerp(vf, p) + p1 = p - 0.5 * dt_ * v1 + v2 = bilerp(vf, p1) + p2 = p - 0.75 * dt_ * v2 + v3 = bilerp(vf, p2) + p -= dt_ * ((2 / 9) * v1 + (1 / 3) * v2 + (4 / 9) * v3) + return p + + +@ti.kernel +def advect(vf: ti.template(), qf: ti.template(), new_qf: ti.template()): + for i, j in vf: + p = ti.Vector([i, j]) + 0.5 + p = backtrace(vf, p, dt) + new_qf[i, j] = bilerp(qf, p) * dye_decay + + +@ti.kernel +def apply_impulse(vf: ti.template(), dyef: ti.template(), imp_data: ti.types.ndarray()): + g_dir = -ti.Vector([0, 9.8]) * 300 + for i, j in vf: + omx, omy = imp_data[2], imp_data[3] + mdir = ti.Vector([imp_data[0], imp_data[1]]) + dx, dy = (i + 0.5 - omx), (j + 0.5 - omy) + d2 = dx * dx + dy * dy + # dv = F * dt + factor = ti.exp(-d2 / force_radius) + + dc = dyef[i, j] + a = dc.norm() + + momentum = (mdir * f_strength * factor + g_dir * a / (1 + a)) * dt + + v = vf[i, j] + vf[i, j] = v + momentum + # add dye + if mdir.norm() > 0.5: + dc += ti.exp(-d2 * (4 / (res / 15) ** 2)) * ti.Vector([imp_data[4], imp_data[5], imp_data[6]]) + + dyef[i, j] = dc + + +@ti.kernel +def divergence(vf: ti.template()): + for i, j in vf: + vl = sample(vf, i - 1, j) + vr = sample(vf, i + 1, j) + vb = sample(vf, i, j - 1) + vt = sample(vf, i, j + 1) + vc = sample(vf, i, j) + if i == 0: + vl.x = -vc.x + if i == res - 1: + vr.x = -vc.x + if j == 0: + vb.y = -vc.y + if j == res - 1: + vt.y = -vc.y + velocity_divs[i, j] = (vr.x - vl.x + vt.y - vb.y) * 0.5 + + +@ti.kernel +def vorticity(vf: ti.template()): + for i, j in vf: + vl = sample(vf, i - 1, j) + vr = sample(vf, i + 1, j) + vb = sample(vf, i, j - 1) + vt = sample(vf, i, j + 1) + velocity_curls[i, j] = (vr.y - vl.y - vt.x + vb.x) * 0.5 + + +@ti.kernel +def pressure_jacobi(pf: ti.template(), new_pf: ti.template()): + for i, j in pf: + pl = sample(pf, i - 1, j) + pr = sample(pf, i + 1, j) + pb = sample(pf, i, j - 1) + pt = sample(pf, i, j + 1) + div = velocity_divs[i, j] + new_pf[i, j] = (pl + pr + pb + pt - div) * 0.25 + + +@ti.kernel +def subtract_gradient(vf: ti.template(), pf: ti.template()): + for i, j in vf: + pl = sample(pf, i - 1, j) + pr = sample(pf, i + 1, j) + pb = sample(pf, i, j - 1) + pt = sample(pf, i, j + 1) + vf[i, j] -= 0.5 * ti.Vector([pr - pl, pt - pb]) + + +@ti.kernel +def enhance_vorticity(vf: ti.template(), cf: ti.template()): + # anti-physics visual enhancement... + for i, j in vf: + cl = sample(cf, i - 1, j) + cr = sample(cf, i + 1, j) + cb = sample(cf, i, j - 1) + ct = sample(cf, i, j + 1) + cc = sample(cf, i, j) + force = ti.Vector([abs(ct) - abs(cb), abs(cl) - abs(cr)]).normalized(1e-3) + force *= curl_strength * cc + vf[i, j] = ti.min(ti.max(vf[i, j] + force * dt, -1e3), 1e3) + + +@ti.kernel +def copy_divergence(div_in: ti.template(), div_out: ti.types.ndarray()): + for I in ti.grouped(div_in): + div_out[I[0] * res + I[1]] = -div_in[I] + + +@ti.kernel +def apply_pressure(p_in: ti.types.ndarray(), p_out: ti.template()): + for I in ti.grouped(p_out): + p_out[I] = p_in[I[0] * res + I[1]] + + +def solve_pressure_sp_mat(): + copy_divergence(velocity_divs, F_b) + x = solver.solve(F_b) + apply_pressure(x, pressures_pair.cur) + + +def solve_pressure_jacobi(): + for _ in range(p_jacobi_iters): + pressure_jacobi(pressures_pair.cur, pressures_pair.nxt) + pressures_pair.swap() + + +def step(mouse_data): + advect(velocities_pair.cur, velocities_pair.cur, velocities_pair.nxt) + advect(velocities_pair.cur, dyes_pair.cur, dyes_pair.nxt) + velocities_pair.swap() + dyes_pair.swap() + + apply_impulse(velocities_pair.cur, dyes_pair.cur, mouse_data) + + divergence(velocities_pair.cur) + + if curl_strength: + vorticity(velocities_pair.cur) + enhance_vorticity(velocities_pair.cur, velocity_curls) + + if use_sparse_matrix: + solve_pressure_sp_mat() + else: + solve_pressure_jacobi() + + subtract_gradient(velocities_pair.cur, pressures_pair.cur) + + if debug: + divergence(velocities_pair.cur) + div_s = np.sum(velocity_divs.to_numpy()) + print(f"divergence={div_s}") + + +class MouseDataGen: + def __init__(self): + self.prev_mouse = None + self.prev_color = None + + def __call__(self, gui): + # [0:2]: normalized delta direction + # [2:4]: current mouse xy + # [4:7]: color + mouse_data = np.zeros(8, dtype=np.float32) + if gui.is_pressed(ti.GUI.LMB): + mxy = np.array(gui.get_cursor_pos(), dtype=np.float32) * res + if self.prev_mouse is None: + self.prev_mouse = mxy + # Set lower bound to 0.3 to prevent too dark colors + self.prev_color = (np.random.rand(3) * 0.7) + 0.3 + else: + mdir = mxy - self.prev_mouse + mdir = mdir / (np.linalg.norm(mdir) + 1e-5) + mouse_data[0], mouse_data[1] = mdir[0], mdir[1] + mouse_data[2], mouse_data[3] = mxy[0], mxy[1] + mouse_data[4:7] = self.prev_color + self.prev_mouse = mxy + else: + self.prev_mouse = None + self.prev_color = None + return mouse_data + + +def reset(): + velocities_pair.cur.fill(0) + pressures_pair.cur.fill(0) + dyes_pair.cur.fill(0) + + +def main(): + global debug, curl_strength + visualize_d = True # visualize dye (default) + visualize_v = False # visualize velocity + visualize_c = False # visualize curl + + paused = False + + gui = ti.GUI("Stable Fluid", (res, res)) + md_gen = MouseDataGen() + + while gui.running: + if gui.get_event(ti.GUI.PRESS): + e = gui.event + if e.key == ti.GUI.ESCAPE: + break + elif e.key == "r": + paused = False + reset() + elif e.key == "s": + if curl_strength: + curl_strength = 0 + else: + curl_strength = 7 + elif e.key == "v": + visualize_v = True + visualize_c = False + visualize_d = False + elif e.key == "d": + visualize_d = True + visualize_v = False + visualize_c = False + elif e.key == "c": + visualize_c = True + visualize_d = False + visualize_v = False + elif e.key == "p": + paused = not paused + elif e.key == "d": + debug = not debug + + # Debug divergence: + # print(max((abs(velocity_divs.to_numpy().reshape(-1))))) + + if not paused: + mouse_data = md_gen(gui) + step(mouse_data) + if visualize_c: + vorticity(velocities_pair.cur) + gui.set_image(velocity_curls.to_numpy() * 0.03 + 0.5) + elif visualize_d: + gui.set_image(dyes_pair.cur) + elif visualize_v: + gui.set_image(velocities_pair.cur.to_numpy() * 0.01 + 0.5) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/two_stream_instability.py b/modules/createEVENT/TaichiEvent/two_stream_instability.py new file mode 100644 index 000000000..3a618d065 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/two_stream_instability.py @@ -0,0 +1,98 @@ +# Authored by Luhuai Jiao +# This is a 1D simulation of "two-stream instability" in Plasma Physicis. +# Some settings of the grids and particles are taken from "Introduction to Computational Plasma Physics"(ISBN: 9787030563675) + +import taichi as ti + +ti.init(arch=ti.gpu) # Try to run on GPU +PI = 3.141592653589793 +L = 8 * PI # simulation domain size +dt = 0.1 # time step +substepping = 8 +ng = 32 # number of grids +np = 16384 # numer of particles +vb = 1.0 # beam-velocity, one is vb, the other is -vb +vt = 0.3 # thermal velocity +wp = 1 # Plasma frequency +qm = -1 # charge-mass ratio +q = wp * wp / (qm * np / L) # charge of a particle +rho_back = -q * np / L # background charge density +dx = L / ng # grid spacing +inv_dx = 1.0 / dx +x = ti.Vector.field(1, ti.f32, np) # position +v = ti.Vector.field(1, ti.f32, np) # velocity +rho = ti.Vector.field(1, ti.f32, ng) # charge density +e = ti.Vector.field(1, ti.f32, ng) # electric fields +# to show x-vx on the screen +v_x_pos1 = ti.Vector.field(2, ti.f32, int(np / 2)) +v_x_pos2 = ti.Vector.field(2, ti.f32, int(np / 2)) + + +@ti.kernel +def initialize(): + for p in x: + x[p].x = (p + 1) * L / np + v[p].x = vt * ti.randn() + (-1) ** p * vb # two streams + + +@ti.kernel +def substep(): + for p in x: + x[p] += v[p] * dt + if x[p].x >= L: # periodic boundary condition + x[p] += -L + if x[p].x < 0: + x[p] += L + rho.fill(rho_back) # fill rho with background charge density + for p in x: # Particle state update and scatter to grid (P2G) + base = (x[p] * inv_dx - 0.5).cast(int) + fx = x[p] * inv_dx - 0.5 - base.cast(float) + rho[base] += (1.0 - fx) * q * inv_dx + if base[0] < ng - 1: + rho[base + 1] += fx * q * inv_dx + e.fill(0.0) + ti.loop_config(serialize=True) + for i in range(ng): # compute electric fields + if i == 0: + e[i] = rho[i] * dx * 0.5 + else: + e[i] = e[i - 1] + (rho[i - 1] + rho[i]) * dx * 0.5 + s = 0.0 + for i in e: + s += e[i].x + for i in e: + e[i] += -s / ng + for p in v: # G2P + base = (x[p] * inv_dx - 0.5).cast(int) + fx = x[p] * inv_dx - 0.5 - base.cast(float) + a = e[base] * (1.0 - fx) * qm + if base[0] < ng - 1: + a += e[base + 1] * fx * qm # compute electric force + v[p] += a * dt + + +@ti.kernel +def vx_pos(): # to show x-vx on the screen + for p in x: + if p % 2: + v_x_pos1[int((p - 1) / 2)].x = x[p].x / L + v_x_pos1[int((p - 1) / 2)].y = (v[p].x) / 10 + 0.5 + else: + v_x_pos2[int(p / 2)].x = x[p].x / L + v_x_pos2[int(p / 2)].y = (v[p].x) / 10 + 0.5 + + +def main(): + initialize() + gui = ti.GUI("Shortest PIC", (800, 800)) + while not gui.get_event(ti.GUI.ESCAPE, ti.GUI.EXIT): + for s in range(substepping): + substep() + vx_pos() + gui.circles(v_x_pos1.to_numpy(), color=0x0000FF, radius=2) + gui.circles(v_x_pos2.to_numpy(), color=0xFF0000, radius=2) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/vortex_rings.py b/modules/createEVENT/TaichiEvent/vortex_rings.py new file mode 100644 index 000000000..c87a34366 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/vortex_rings.py @@ -0,0 +1,97 @@ +# C++ reference and tutorial (Chinese): https://zhuanlan.zhihu.com/p/26882619 +import math + +import numpy as np + +import taichi as ti + +ti.init(arch=ti.gpu) + +eps = 0.01 +dt = 0.1 + +n_vortex = 4 +n_tracer = 200000 + +pos = ti.Vector.field(2, ti.f32, shape=n_vortex) +new_pos = ti.Vector.field(2, ti.f32, shape=n_vortex) +vort = ti.field(ti.f32, shape=n_vortex) + +tracer = ti.Vector.field(2, ti.f32, shape=n_tracer) + + +@ti.func +def compute_u_single(p, i): + r2 = (p - pos[i]).norm() ** 2 + uv = ti.Vector([pos[i].y - p.y, p.x - pos[i].x]) + return vort[i] * uv / (r2 * math.pi) * 0.5 * (1.0 - ti.exp(-r2 / eps**2)) + + +@ti.func +def compute_u_full(p): + u = ti.Vector([0.0, 0.0]) + for i in range(n_vortex): + u += compute_u_single(p, i) + return u + + +@ti.kernel +def integrate_vortex(): + for i in range(n_vortex): + v = ti.Vector([0.0, 0.0]) + for j in range(n_vortex): + if i != j: + v += compute_u_single(pos[i], j) + new_pos[i] = pos[i] + dt * v + + for i in range(n_vortex): + pos[i] = new_pos[i] + + +@ti.kernel +def advect(): + for i in range(n_tracer): + # Ralston's third-order method + p = tracer[i] + v1 = compute_u_full(p) + v2 = compute_u_full(p + v1 * dt * 0.5) + v3 = compute_u_full(p + v2 * dt * 0.75) + tracer[i] += (2 / 9 * v1 + 1 / 3 * v2 + 4 / 9 * v3) * dt + + +pos[0] = [0, 1] +pos[1] = [0, -1] +pos[2] = [0, 0.3] +pos[3] = [0, -0.3] +vort[0] = 1 +vort[1] = -1 +vort[2] = 1 +vort[3] = -1 + + +@ti.kernel +def init_tracers(): + for i in range(n_tracer): + tracer[i] = [ti.random() - 0.5, ti.random() * 3 - 1.5] + + +def main(): + init_tracers() + gui = ti.GUI("Vortex Rings", (1024, 512), background_color=0xFFFFFF) + + while gui.running: + for i in range(4): # substeps + advect() + integrate_vortex() + + gui.circles( + tracer.to_numpy() * np.array([[0.05, 0.1]]) + np.array([[0.0, 0.5]]), + radius=0.5, + color=0x0, + ) + + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/createEVENT/TaichiEvent/waterwave.py b/modules/createEVENT/TaichiEvent/waterwave.py new file mode 100644 index 000000000..73a238ba3 --- /dev/null +++ b/modules/createEVENT/TaichiEvent/waterwave.py @@ -0,0 +1,100 @@ +# Water wave effect partially based on shallow water equations +# https://en.wikipedia.org/wiki/Shallow_water_equations#Non-conservative_form + +import taichi as ti + +ti.init(arch=ti.gpu) + +light_color = 1 +gravity = 2.0 # larger gravity makes wave propagates faster +damping = 0.2 # larger damping makes wave vanishes faster when propagating +dx = 0.02 +dt = 0.01 +shape = 512, 512 +pixels = ti.field(dtype=float, shape=shape) +background = ti.field(dtype=float, shape=shape) +height = ti.field(dtype=float, shape=shape) +velocity = ti.field(dtype=float, shape=shape) + + +@ti.kernel +def reset(): + for i, j in height: + t = i // 16 + j // 16 + if t % 2 == 0: + background[i, j] = 0.0 + else: + background[i, j] = 0.25 + height[i, j] = 0 + velocity[i, j] = 0 + + +@ti.func +def laplacian(i, j): + return (-4 * height[i, j] + height[i, j - 1] + height[i, j + 1] + height[i + 1, j] + height[i - 1, j]) / ( + 4 * dx**2 + ) + + +@ti.func +def gradient(i, j): + return ti.Vector( + [ + (height[i + 1, j] if i < shape[0] - 1 else 0) - (height[i - 1, j] if i > 1 else 0), + (height[i, j + 1] if j < shape[1] - 1 else 0) - (height[i, j - 1] if j > 1 else 0), + ] + ) * (0.5 / dx) + + +@ti.kernel +def create_wave(amplitude: ti.f32, x: ti.f32, y: ti.f32): + for i, j in ti.ndrange((1, shape[0] - 1), (1, shape[1] - 1)): + r2 = (i - x) ** 2 + (j - y) ** 2 + height[i, j] = height[i, j] + amplitude * ti.exp(-0.02 * r2) + + +@ti.kernel +def update(): + for i, j in ti.ndrange((1, shape[0] - 1), (1, shape[1] - 1)): + acceleration = gravity * laplacian(i, j) - damping * velocity[i, j] + velocity[i, j] = velocity[i, j] + acceleration * dt + + for i, j in ti.ndrange((1, shape[0] - 1), (1, shape[1] - 1)): + height[i, j] = height[i, j] + velocity[i, j] * dt + + +@ti.kernel +def visualize_wave(): + # visualizes the wave using a fresnel-like shading + # a brighter color indicates a steeper wave + # (closer to grazing angle when looked from above) + for i, j in pixels: + g = gradient(i, j) + cos_i = 1 / ti.sqrt(1 + g.norm_sqr()) + brightness = pow(1 - cos_i, 2) + color = background[i, j] + pixels[i, j] = (1 - brightness) * color + brightness * light_color + + +def main(): + print("[Hint] click on the window to create waves") + + reset() + gui = ti.GUI("Water Wave", shape) + while gui.running: + for e in gui.get_events(ti.GUI.PRESS): + if e.key in [ti.GUI.ESCAPE, ti.GUI.EXIT]: + gui.running = False + elif e.key == "r": + reset() + elif e.key == ti.GUI.LMB: + x, y = e.pos + create_wave(3, x * shape[0], y * shape[1]) + update() + visualize_wave() + gui.set_image(pixels) + gui.show() + + +if __name__ == "__main__": + main() diff --git a/modules/performRegionalEventSimulation/regionalWindField/WindFieldModel.h b/modules/performRegionalEventSimulation/regionalWindField/WindFieldModel.h index 8a4cede10..b4e5de271 100644 --- a/modules/performRegionalEventSimulation/regionalWindField/WindFieldModel.h +++ b/modules/performRegionalEventSimulation/regionalWindField/WindFieldModel.h @@ -1,5 +1,5 @@ #ifndef WIND_FIELD_MODEL_H_ -#define WIND_FILED_MODEL_H_ +#define WIND_FIELD_MODEL_H_ #include @@ -89,4 +89,4 @@ class WindFieldModel }; -#endif +#endif // WIND_FIELD_MODEL_H_ diff --git a/pyproject.toml b/pyproject.toml index c6dbc62a2..3501f12c6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,7 @@ [tool.ruff] line-length = 85 +# Example Taichi programs for users to test are held here until formal examples are made/documented for HydroUQ, etc. Don't need to lint/format - JustinBonus +extend-exclude = ["modules/createEVENT/TaichiEvent"] [tool.ruff.lint] # Enable all known categories