Skip to content

Conversation

@Gauravu2
Copy link
Collaborator

@Gauravu2 Gauravu2 commented Sep 8, 2022

Resolves #51

trampoline.mp4

@Gauravu2 Gauravu2 requested a review from bhosale2 September 8, 2022 22:16
@Gauravu2 Gauravu2 self-assigned this Sep 8, 2022
@bhosale2 bhosale2 added enhancement New feature or request prio:high High priority labels Sep 8, 2022
Comment on lines +5 to +48
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
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Comment on lines +5 to +30
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)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Comment on lines +5 to +34
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
)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please rename the variables here referring to Z and R, and not X and Y, for better clarity.

Comment on lines +364 to +366
fx, fy = collision_force(
1 * G, ball_phi2, ball_phi1, 2 * moll_zone, grid_size_z, grid_size_r, eps, dx
)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why 0 * fy?

Comment on lines +376 to +385
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
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

@Gauravu2 Gauravu2 removed their assignment Sep 12, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request prio:high High priority

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Rigid sphere falling on an elastic trampoline

2 participants