diff --git a/examples/RigidBallSoftTrampoline/collision_force.py b/examples/RigidBallSoftTrampoline/collision_force.py new file mode 100644 index 0000000..f8e59c5 --- /dev/null +++ b/examples/RigidBallSoftTrampoline/collision_force.py @@ -0,0 +1,48 @@ +import numpy as np +import scipy.sparse as spp + + +def collision_force(k, phi_1, phi_2, blend_w, grid_size_z, grid_size_r, eps, dx): + """ + computes the collision force between 2 bodies + using their level sets + """ + + gradient = ( + spp.diags( + [-0.5, 0.5, -0.5, 0.5], + [-1, 1, grid_size_r - 1, -grid_size_r + 1], + shape=(grid_size_r, grid_size_r), + format="csr", + ) + / dx + ) + + gradient_T = ( + spp.diags( + [-0.5, 0.5, -0.5, 0.5], + [-1, 1, grid_size_z - 1, -grid_size_z + 1], + shape=(grid_size_z, grid_size_z), + format="csr", + ) + / dx + ) + + phi_diff = 0.5 * (phi_2 - phi_1) + delta_phi = phi_diff * 0 + delta_phi += ( + (np.fabs(phi_diff) < blend_w) + * 0.5 + * (1 + np.cos(np.pi * phi_diff / blend_w)) + / blend_w + ) + phi_diff_x = np.mat(phi_diff) * gradient_T + phi_diff_y = gradient * np.mat(phi_diff) + s = np.sqrt(np.square(phi_diff_x) + np.square(phi_diff_y)) + phi_diff_x /= s + eps + phi_diff_y /= s + eps + S = phi_diff / np.sqrt(phi_diff**2 + eps**2) + body_mask = np.logical_or((phi_1 > 0), (phi_2 > 0)) + force_x = body_mask * S * phi_diff_x * k * delta_phi + force_y = body_mask * S * phi_diff_y * k * delta_phi + return force_x, force_y diff --git a/examples/RigidBallSoftTrampoline/curl.py b/examples/RigidBallSoftTrampoline/curl.py new file mode 100644 index 0000000..23fcbea --- /dev/null +++ b/examples/RigidBallSoftTrampoline/curl.py @@ -0,0 +1,31 @@ +import numpy as np +import scipy.sparse as spp + + +def curl(ax, ay, grid_size_z, grid_size_r, dx): + """ + compute curl of the vector using CD + """ + + gradient = ( + spp.diags( + [-0.5, 0.5, -0.5, 0.5], + [-1, 1, grid_size_r - 1, -grid_size_r + 1], + shape=(grid_size_r, grid_size_r), + format="csr", + ) + / dx + ) + + gradient_T = ( + spp.diags( + [-0.5, 0.5, -0.5, 0.5], + [-1, 1, grid_size_z - 1, -grid_size_z + 1], + shape=(grid_size_z, grid_size_z), + format="csr", + ) + / dx + ) + + vort = np.mat(ay) * gradient_T - gradient * np.mat(ax) + return vort diff --git a/examples/RigidBallSoftTrampoline/hard_sphere_on_trampoline_colliding.py b/examples/RigidBallSoftTrampoline/hard_sphere_on_trampoline_colliding.py new file mode 100644 index 0000000..11eeee0 --- /dev/null +++ b/examples/RigidBallSoftTrampoline/hard_sphere_on_trampoline_colliding.py @@ -0,0 +1,399 @@ +import numpy as np +import matplotlib.pyplot as plt +import skfmm +import os +from pyaxisymflow.utils.custom_cmap import lab_cmp +from pyaxisymflow.kernels.brinkmann_penalize import brinkmann_penalize +from pyaxisymflow.kernels.compute_forces import compute_force_on_body +from pyaxisymflow.kernels.compute_velocity_from_psi import compute_velocity_from_psi_unb +from pyaxisymflow.kernels.smooth_Heaviside import smooth_Heaviside +from pyaxisymflow.kernels.kill_boundary_vorticity_sine import ( + kill_boundary_vorticity_sine_r, + kill_boundary_vorticity_sine_z, +) +from pyaxisymflow.kernels.FastDiagonalisationStokesSolver import ( + FastDiagonalisationStokesSolver, +) +from pyaxisymflow.kernels.diffusion_RK2 import diffusion_RK2_unb +from pyaxisymflow.kernels.advect_vorticity_via_eno3 import gen_advect_vorticity_via_eno3 +from pyaxisymflow.elasto_kernels.div_tau import update_vorticity_from_solid_stress +from pyaxisymflow.elasto_kernels.solid_sigma import solid_sigma +from pyaxisymflow.elasto_kernels.extrapolate_eta_using_least_squares_unb import ( + extrapolate_eta_with_least_squares, +) +from pyaxisymflow.elasto_kernels.advect_refmap_via_eno3 import ( + gen_advect_refmap_via_eno3, +) +from curl import curl +from collision_force import collision_force +from pyaxisymflow.kernels.compute_velocity_from_psi import compute_velocity_from_psi_unb +from pyaxisymflow.kernels.compute_vorticity_from_velocity import ( + compute_vorticity_from_velocity_unb, +) +from pyaxisymflow.kernels.compute_vorticity_from_velocity import ( + compute_vorticity_from_velocity_unb, +) +from pyaxisymflow.utils.dump_vtk import vtk_init, vtk_write +from trampoline_level_set import trampoline_level_set + +# global settings +grid_size_z = 400 +domain_AR = 0.5 +dx = 1.0 / grid_size_z +grid_size_r = int(domain_AR * grid_size_z) +CFL = 0.1 +eps = np.finfo(float).eps +num_threads = 4 +implicit_diffusion = False + + +plt.figure(figsize=(5 / domain_AR, 5)) +# Parameters +fotoTimer_limit = 0.01 +moll_zone = dx * np.sqrt(16) +brink_lam = 1e4 +extrap_zone = moll_zone + 4 * dx +reinit_band = extrap_zone +scale = 1 +r_ball = scale * 0.075 +U_0 = 1.0 +Re = 100.0 +nu = 0.001 +nondim_T = 3000 +tEnd = 10 +T_ramp = 20 * r_ball / U_0 +rho_f = 1 +G = 0.1 +# Build discrete domain +z = np.linspace(0 + dx / 2, 1 - dx / 2, grid_size_z) +r = np.linspace(0 + dx / 2, domain_AR - dx / 2, grid_size_r) +Z, R = np.meshgrid(z, r) +K_tether = G / (dx**2) +# load initial conditions +vorticity = 0 * Z +penal_vorticity = 0 * Z +penal_vorticity_t = 0 * Z +temp_vorticity = 0 * Z +psi = 0 * Z +u_z = 0 * Z +u_r = 0 * Z +u_z_upen = 0 * Z +u_r_upen = 0 * Z +u_z_upen_t = 0 * Z +u_r_upen_t = 0 * Z + + +fotoTimer = 0.0 +it = 0 +Z_cm1 = 0.60 +Z_cm2 = 0.30 +R_cm = 0.0 +t = 0.0 +T = [] +rho_s = 1.05 +g = 5 * (rho_s - rho_f) * 9.81 +# create char function +ball_phi2 = -np.sqrt((Z - Z_cm2) ** 2 + (R - R_cm) ** 2) + r_ball +ball_char_func2 = 0 * Z +smooth_Heaviside(ball_char_func2, ball_phi2, moll_zone) +inside_solid2 = ball_char_func2 > 0 + +# trampoline phi +trampoline_CM_Z = 0.6 +trampoline_diameter = scale * 0.4 +trampoline_height = scale * 0.05 +left_center_R = -trampoline_diameter / 2 +left_center_Z = trampoline_CM_Z +corner_radius = 0.5 * trampoline_height +ball_phi1 = -trampoline_level_set( + R, + Z, + left_center_R, + trampoline_diameter, + trampoline_CM_Z, + corner_radius, + corner_radius, +) +ball_char_func1 = 0 * Z +smooth_Heaviside(ball_char_func1, ball_phi1, moll_zone) +inside_solid1 = ball_char_func1 > 0 +Z_cm = trampoline_CM_Z +R_cm_t = trampoline_diameter / 2 +fixed_rad = 0.01 +tether_phi = -np.sqrt((Z - Z_cm) ** 2 + (R - R_cm_t) ** 2) + fixed_rad +tether_char_func = 0 * Z +smooth_Heaviside(tether_char_func, tether_phi, moll_zone) + +part_mass = 2 * np.pi * dx * dx * rho_s * np.sum(ball_char_func2 * R) +part_vol = 2 * np.pi * dx * dx * np.sum(ball_char_func2 * R) +part_Z_cm_old = Z_cm2 +part_Z_cm_new = Z_cm2 +part_Z_cm = Z_cm2 + +advect_refmap_via_eno3 = gen_advect_refmap_via_eno3(dx, grid_size_r, grid_size_z) +FD_stokes_solver = FastDiagonalisationStokesSolver(grid_size_r, grid_size_z, dx) +advect_vorticity_via_eno3 = gen_advect_vorticity_via_eno3( + dx, grid_size_r, grid_size_z, num_threads=num_threads +) + +diffusion_dt_limit = dx / np.sqrt(G) +if implicit_diffusion: + implicit_diffusion_stepper = ImplicitEulerDiffusionStepper( + time_step=diffusion_dt_limit, + kinematic_viscosity=nu, + grid_size_r=grid_size_r, + grid_size_z=grid_size_z, + dx=dx, + ) +r1 = np.linspace(0 + dx / 2, domain_AR - dx / 2, 2 * grid_size_r) +Z_double, R_double = np.meshgrid(z, r1) +F_un = 0 +F_pen = 0.0 +U_z_cm_part = 0 +U_z_cm_part_old = 0 +diff = 0 +bad_phi = 0 * Z +phi_orig = 0 * Z +total_flux_double = 0 * R_double +sigma_s_11 = 0 * Z +sigma_s_12 = 0 * Z +sigma_s_22 = 0 * Z +eta1z = 0 * Z +eta2z = 0 * Z +eta1r = 0 * Z +eta2r = 0 * Z +tau_z = 0 * Z +tau_r = 0 * Z +eta1 = Z.copy() +eta2 = R.copy() +eta1_double = Z_double.copy() +eta2_double = R_double.copy() +ball_phi_double = 0 * Z_double +part_Z_cm_old = Z_cm2 +part_Z_cm_new = Z_cm2 +part_Z_cm = Z_cm2 +F_un = 0 +F_pen = 0.0 +U_z_cm_part = 0 +U_z_cm_part_old = 0 +diff = 0 +T = [] +U_z_cm = [] +vtk_image_data, temp_vtk_array, writer = vtk_init(grid_size_z, grid_size_r) + +# solver loop +while t < tEnd: + + # kill vorticity at boundaries + kill_boundary_vorticity_sine_z(vorticity, Z, 3, dx) + kill_boundary_vorticity_sine_r(vorticity, R, 3, dx) + + # solve for stream function and get velocity + FD_stokes_solver.solve(solution_field=psi, rhs_field=vorticity) + compute_velocity_from_psi_unb(u_z, u_r, psi, R, dx) + + if fotoTimer >= fotoTimer_limit or t == 0: + fotoTimer = 0.0 + levels = np.linspace(-0.1, 0.1, 25) + plt.contourf( + Z, + R, + vorticity, + levels=100, + extend="both", + cmap=lab_cmp, + ) + + plt.contour( + Z, + R, + ball_char_func1, + levels=[ + 0.5, + ], + colors="grey", + ) + plt.contour( + Z, + R, + ball_char_func2, + levels=[ + 0.5, + ], + colors="grey", + ) + plt.contour( + Z, + R, + tether_char_func, + levels=[ + 0.5, + ], + colors="grey", + ) + + plt.xticks([]) + plt.yticks([]) + plt.gca().set_aspect("equal") + plt.savefig("snap_" + str("%0.4d" % (t * 100)) + ".png") + plt.clf() + plt.close("all") + vtk_write( + "axisym_avg_" + str("%0.4d" % (t * 100)) + ".vti", + vtk_image_data, + temp_vtk_array, + writer, + ["ball_char_func1", "ball_char_func2", "vorticity", "psi", "u_z", "u_r"], + [ball_char_func1, ball_char_func2, vorticity, psi, u_z, u_r], + grid_size_z, + grid_size_r, + ) + + # get dt + if implicit_diffusion: + # technically we can set any higher dt here lower than the CFL limit + dt = diffusion_dt_limit + else: + dt = min( + CFL * dx / np.sqrt(G / rho_f), + 0.9 * dx**2 / 4 / nu, + CFL / (np.amax(np.fabs(u_z) + np.fabs(u_r)) + eps), + ) + + # advect refmap + advect_refmap_via_eno3( + eta1, + eta2, + u_z, + u_r, + dt, + ) + + compute_vorticity_from_velocity_unb( + penal_vorticity_t, + tether_char_func * (eta1 - Z), + tether_char_func * (eta2 - R), + dx, + ) + vorticity[...] += K_tether * dt * penal_vorticity_t + + # penalise velocity (particle) + u_z_upen[...] = u_z.copy() + u_r_upen[...] = u_r.copy() + brinkmann_penalize( + brink_lam, dt, ball_char_func2, U_z_cm_part, 0.0, u_z_upen, u_r_upen, u_z, u_r + ) + compute_vorticity_from_velocity_unb( + penal_vorticity, u_z - u_z_upen, u_r - u_r_upen, dx + ) + vorticity[...] += penal_vorticity + + F_pen, F_un = compute_force_on_body( + R, ball_char_func2, rho_f, brink_lam, u_z, U_z_cm_part, part_vol, dt, diff + ) + + ball_phi2 = -np.sqrt((Z - part_Z_cm) ** 2 + (R - R_cm) ** 2) + r_ball + ball_char_func2 = 0 * Z + smooth_Heaviside(ball_char_func2, ball_phi2, moll_zone) + + # pin eta and phi boundary + phi_orig[...] = -trampoline_level_set( + eta2, + eta1, + left_center_R, + trampoline_diameter, + trampoline_CM_Z, + corner_radius, + corner_radius, + ) + band = np.where(ball_phi1 > -3 * dx) # assuming solid_CFL << 1 + ball_phi1[band] = phi_orig[band] + + # reinit level set + bad_phi[...] = ball_phi1 + ball_phi1 = skfmm.distance(ball_phi1, dx=dx, narrow=reinit_band) + idx = ball_phi1.mask + idx_one = np.logical_and(idx, bad_phi < 0.0) + idx_two = np.logical_and(idx, bad_phi > 0.0) + ball_phi1[idx_one] = -99.0 + ball_phi1[idx_two] = 99.0 + + # get char function + ball_char_func1 = 0 * Z + smooth_Heaviside(ball_char_func1, ball_phi1, moll_zone) + inside_solid1[...] = ball_phi1 > 0.0 + + # extrapolate eta for stresses + extrapolate_eta_with_least_squares( + inside_solid1, + ball_phi1, + eta1, + eta2, + ball_phi_double, + eta1_double, + eta2_double, + extrap_zone, + grid_size_r, + z, + ) + + # compute solid stresses and blend + solid_sigma( + sigma_s_11, + sigma_s_12, + sigma_s_22, + G, + dx, + eta1, + eta2, + eta1z, + eta1r, + eta2z, + eta2r, + ) + + sigma_s_11[...] = ball_char_func1 * sigma_s_11 + sigma_s_12[...] = ball_char_func1 * sigma_s_12 + sigma_s_22[...] = ball_char_func1 * sigma_s_22 + + update_vorticity_from_solid_stress( + vorticity, tau_z, tau_r, sigma_s_11, sigma_s_12, sigma_s_22, R, dt, dx + ) + + # Calculate collision forces + fx, fy = collision_force( + 1 * G, ball_phi2, ball_phi1, 2 * moll_zone, grid_size_z, grid_size_r, eps, dx + ) + vorticity[...] += dt * curl(fx, 0 * fy, grid_size_z, grid_size_r, dx) + + # diffuse vorticity + diffusion_RK2_unb(vorticity, temp_vorticity, R, nu, dt, dx) + + # advect vorticity + advect_vorticity_via_eno3(vorticity, u_z, u_r, dt) + + # update particle location and velocity + F_total = 2 * np.pi * dx * dx * F_pen + F_un - np.sum(R * ball_char_func2 * fx) + U_z_cm_part_old = U_z_cm_part + U_z_cm_part += dt * ((g / rho_s + F_total / part_mass)) + diff = dt * (g / rho_s + F_total / part_mass) + part_Z_cm_new = part_Z_cm + Z_cm = part_Z_cm + part_Z_cm += U_z_cm_part_old * dt + ( + 0.5 * dt * dt * (g / rho_s + F_total / part_mass) + ) + part_Z_cm_old = part_Z_cm_new + + # update time + t += dt + fotoTimer += dt + it += 1 + if it % 50 == 0: + print(f"time: {t}, max vort: {np.amax(vorticity)},simulations:{U_z_cm_part}") + + +os.system("rm -f 2D_advect.mp4") +os.system( + "ffmpeg -r 20 -s 3840x2160 -f image2 -pattern_type glob -i 'snap*.png' " + "-vcodec libx264 -crf 15 -pix_fmt yuv420p 2D_advect.mp4" +) diff --git a/examples/RigidBallSoftTrampoline/trampoline_level_set.py b/examples/RigidBallSoftTrampoline/trampoline_level_set.py new file mode 100644 index 0000000..a992521 --- /dev/null +++ b/examples/RigidBallSoftTrampoline/trampoline_level_set.py @@ -0,0 +1,156 @@ +import numpy as np +from matplotlib import pyplot as plt + + +def trampoline_level_set( + inp_X, + inp_Y, + center_basal_x, + distance_between_centers_x, + center_y, + r_basal, + r_tip, +): + """Gives the level set of a Trampoline""" + # 1. Basal circle + # 2. Straight line above and below symmetry plane + # 3. Tip circle + # center_y = 0.5 + # r_basal = 0.1 + # r_tip = 0.05 + # distance_between_centers_x = 0.3 + # center_basal_x = 0.2 + center_tip_x = center_basal_x + distance_between_centers_x + phi = 0.0 * inp_X + + basal_theta = np.arccos((r_basal - r_tip) / distance_between_centers_x) + THETA = np.arctan2((inp_Y - center_y), (inp_X - center_basal_x)) + idx_one = np.abs(THETA) > basal_theta + + # In this region, the level set should be circle distance desu + phi[idx_one] = ( + np.sqrt( + (inp_X[idx_one] - center_basal_x) ** 2 + (inp_Y[idx_one] - center_y) ** 2 + ) + - r_basal + ) + + # Do for the smaller circle next + THETA = np.arctan2((inp_Y - center_y), (inp_X - center_tip_x)) + idx_two = np.abs(THETA) < basal_theta + phi[idx_two] = ( + np.sqrt((inp_X[idx_two] - center_tip_x) ** 2 + (inp_Y[idx_two] - center_y) ** 2) + - r_tip + ) + + # Now shift FOR to the point of attachment and calculate distance to line + attach_basal_x = center_basal_x + r_basal * np.cos(basal_theta) + attach_basal_y = center_y + r_basal * np.sin(basal_theta) + + # in this shifted coordinates, the line is y = -cot(basal_theta)*x + shifted_X = inp_X - attach_basal_x + shifted_Y = inp_Y - attach_basal_y + + cot_basal_theta = 1.0 / np.tan(basal_theta) + dist = (shifted_Y + cot_basal_theta * shifted_X) / np.sqrt(1 + cot_basal_theta**2) + mask_idx = (inp_Y > center_y) & np.logical_not(np.logical_or(idx_one, idx_two)) + phi[mask_idx] = dist[mask_idx] + + # Do the same for down under + attach_basal_x = center_basal_x + r_basal * np.cos(basal_theta) + attach_basal_y = center_y - r_basal * np.sin(basal_theta) + + # in this shifted coordinates, the line is y = cot(basal_theta)*x + shifted_X = inp_X - attach_basal_x + shifted_Y = inp_Y - attach_basal_y + + dist = -(shifted_Y - cot_basal_theta * shifted_X) / np.sqrt( + 1 + cot_basal_theta**2 + ) + mask_idx = (inp_Y < center_y) & np.logical_not(np.logical_or(idx_one, idx_two)) + phi[mask_idx] = dist[mask_idx] + + return phi + + +def __internal_level_set_test(inp_X, inp_Y, phi_input): + import scipy.sparse as spp + + del_x = inp_X[0, 1] - inp_X[0, 0] + del_y = inp_Y[1, 0] - inp_Y[0, 0] + print(del_x, del_y) + grid_size = inp_X.shape + mat_size = grid_size[0] * grid_size[1] + grad_x = spp.diags([-1, 1], [-1, 1], shape=(mat_size, mat_size)) / 2.0 / del_x + grad_y = ( + spp.diags([-1, 1], [-grid_size[1], grid_size[1]], shape=(mat_size, mat_size)) + / 2.0 + / del_y + ) + phi_gradient = ( + ( + grad_x + @ phi.reshape( + -1, + ) + ) + ** 2 + + ( + grad_y + @ phi.reshape( + -1, + ) + ) + ** 2 + ).reshape(grid_size[0], grid_size[1]) + return np.sqrt(phi_gradient) + + +if __name__ == "__main__": + n_points = 128 + del_x = 1.0 / float(n_points) + x = np.linspace(0.5 * del_x, 1.0 - del_x, n_points) + X, Y = np.meshgrid(x, x) + + fig = plt.figure(figsize=(8, 8)) + ax = fig.add_subplot(111) + phi = octopus_arm_level_set( + X, + Y, + center_basal_x=0.6, + distance_between_centers_x=0.25, + center_y=0.4, + r_basal=0.04, + r_tip=0.04, + ) + activation_phi = octopus_arm_level_set( + X, + Y, + center_basal_x=0.6, + distance_between_centers_x=0.25, + center_y=0.4, + r_basal=0.02, + r_tip=0.02, + ) + + # Level set of an ellipse + phi_two = np.fliplr(phi) + phi_ellipse = np.sqrt((X - 0.5) ** 2 / 0.1**2 + (Y - 0.6) ** 2 / 0.18**2) - 1 + + idx = (phi > 0.0) * (phi_two > 0.0) * (phi_ellipse > 0.0) + + phi_new = 0.0 * X + phi_new[idx] = phi[idx] + phi_two[idx] + phi_ellipse[idx] + # ax.contour(X, Y, phi_new, levels=[0], colors="white", linewidths=3.0) + # cont = ax.contourf(X, Y, phi_new, levels=30, cmap=plt.get_cmap("Reds")) + + # phi = Y + X + # ax.contour(X, Y, phi, levels=[0], colors="white", linewidths=3.0) + # cont = ax.contourf(X, Y, phi, levels=30, cmap=plt.get_cmap("Reds")) + + ax.contour(X, Y, phi, levels=[0], colors="white", linewidths=3.0) + cont = ax.contourf(X, Y, phi, levels=30, cmap=plt.get_cmap("Reds")) + ax.contour(X, Y, activation_phi, levels=[0], colors="black", linewidths=3.0) + + ax.set_aspect("equal") + plt.show()