Elastic_trampoline_files#59
Conversation
| 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 |
There was a problem hiding this comment.
please write the vectorized versions (slice-based, similar to functions in kernels), Let's skip using the sparse matrices if possible. Also, collision forces should be a part of the main package.
| 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) |
There was a problem hiding this comment.
can't we use compute_vorticity_from_velocity for this? maybe a good idea is to rename that function as compute_vorticity_from_velocity_forcing?
| 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 | ||
| ) |
There was a problem hiding this comment.
please rename the variables here referring to Z and R, and not X and Y, for better clarity.
| fx, fy = collision_force( | ||
| 1 * G, ball_phi2, ball_phi1, 2 * moll_zone, grid_size_z, grid_size_r, eps, dx | ||
| ) |
There was a problem hiding this comment.
x and y? shouldnt it be z and r?
| 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) |
| 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 |
There was a problem hiding this comment.
these things should be refactored later in a cleaner fashion, for now, let's keep it, but maybe put a TODO note here saying need to be refactored in a cleaner way.
| diff = 0 | ||
| bad_phi = 0 * Z | ||
| phi_orig = 0 * Z | ||
| total_flux_double = 0 * R_double |
There was a problem hiding this comment.
remove variables like this that are not being used. For higher grid resolutions this redundant mem allocation can actually build up quickly and slow down things.
Resolves #51
trampoline.mp4