如何在GLSL计算着色器中模拟类肥皂泡/细胞的球体碰撞？

需求说明 #

在GPU上运行包含1000个粒子的N体模拟系统，粒子间可相互作用，需要实现类似肥皂泡或细胞分裂的效果：

• 碰撞无需刚性，仅需球体间相互推挤
• 整体结构不过度扩张
• 工作组大小设置为10×10×10

现有GLSL计算着色器代码 #

#[compute]
#version 450

layout(local_size_x = 10, local_size_y = 10, local_size_z = 10) in;

layout(set = 0, binding = 0, std430) restrict buffer MyPositionXBuffer {
    float data[];
}
my_position_x_buffer;
layout(set = 0, binding = 1, std430) restrict buffer MyPositionYBuffer {
    float data[];
}
my_position_y_buffer;
layout(set = 0, binding = 2, std430) restrict buffer MyPositionZBuffer {
    float data[];
}
my_position_z_buffer;

layout(set = 0, binding = 3, std430) restrict buffer MyVelocityXBuffer {
    float data[];
}
my_velocity_x_buffer;
layout(set = 0, binding = 4, std430) restrict buffer MyVelocityYBuffer {
    float data[];
}
my_velocity_y_buffer;
layout(set = 0, binding = 5, std430) restrict buffer MyVelocityZBuffer {
    float data[];
}
my_velocity_z_buffer;


float webby_force(float dist) {
    //return clamp(pow(dist, 16.0), 0.0, 10.0);
    float force_reversal_dist = 0.5;
    return -clamp(pow((dist - force_reversal_dist) * 100.0, 3.0), -100.0, 100.0);
}


float attraction(float dist) {
    //return pow(dist - 0.9, 3.0);
    
    dist = clamp(dist, 0.0, 1.0);
    float value = pow(asin(dist - 0.1), 5.0);
    value = clamp(value, -1.0, 1.0);
    return value;
}


float simple_force(float dist) {
    float diff = dist - 0.1;
    float mult = 50.0;
    float value = pow(diff * mult, 5.0) * max((-abs(diff) * 2.0) + 0.05, 0.0);
    value = clamp(value, -1.0, 1.0);
    return -value;
}


float repulsion(float dist) {
    float diff = dist - 0.1;
    return clamp(pow(diff * 100.0, 2.0), -5.0, 5.0);
}


float is_colliding(float dist) {
    return max(-sign(dist - 0.1), 0.0);
}


float neighbor_is_not_self(uint loc, uint loc_neighbor) {
    return abs(sign(loc - loc_neighbor));
}


void main() {
    uint location = ((gl_WorkGroupID.x * gl_NumWorkGroups.x * gl_NumWorkGroups.x) + (gl_WorkGroupID.y * gl_NumWorkGroups.x) + (gl_WorkGroupID.z));
    vec3 position = vec3(my_position_x_buffer.data[location], my_position_y_buffer.data[location], my_position_z_buffer.data[location]);
    vec3 velocity = vec3(my_velocity_x_buffer.data[location], my_velocity_y_buffer.data[location], my_velocity_z_buffer.data[location]);
    
    vec3 neighbor_pos = vec3(my_position_x_buffer.data[gl_LocalInvocationIndex], my_position_y_buffer.data[gl_LocalInvocationIndex], my_position_z_buffer.data[gl_LocalInvocationIndex]);
    uint i = 1000;
    vec3 tinyvec = vec3(my_position_x_buffer.data[i], my_position_y_buffer.data[i], my_position_z_buffer.data[i]);
    vec3 neighbor_vec = normalize(neighbor_pos - position + tinyvec);
    vec3 neighbor_vec_unnorm = neighbor_pos - position + tinyvec;
    
    memoryBarrierShared();
    
    float distance_multiplier = 1.0;
    
    // collision detection
    float neighbor_not_self_flag = float(neighbor_is_not_self(location, gl_LocalInvocationIndex));
    float neighbor_dist = distance(position, neighbor_pos);
    float collision_flag = is_colliding(neighbor_dist * distance_multiplier) * neighbor_not_self_flag;
    float no_collision_flag = 1.0 - collision_flag;
    
    // calculated forces
    float repulsive_factor = repulsion(neighbor_dist * distance_multiplier);
    float webby_factor = webby_force(neighbor_dist * distance_multiplier);
    float simple_force_factor = simple_force(neighbor_dist * distance_multiplier);
    
    memoryBarrierShared();
    
    float force_multiplier = 1.0;
    
    // force integration
    vec3 neighbor_force = neighbor_vec_unnorm * force_multiplier * simple_force_factor * webby_factor;
    vec3 center_force = normalize(position + tinyvec) * length(position) * 0.00001;
    
    memoryBarrierShared();
    
    float damping = 0.9;
    
    velocity += neighbor_force + center_force;
    velocity *= damping;
    
    memoryBarrierShared();
    
    vec3 newposition = position - velocity;
    
    memoryBarrierShared();

    my_position_x_buffer.data[location] = newposition.x;
    my_position_y_buffer.data[location] = newposition.y;
    my_position_z_buffer.data[location] = newposition.z;
    my_velocity_x_buffer.data[location] = velocity.x;
    my_velocity_y_buffer.data[location] = velocity.y;
    my_velocity_z_buffer.data[location] = velocity.z;
}

疑问与请求 #

已尝试多种力计算函数，甚至组合simple_force_factor与webby_factor，但模拟效果始终不符合预期。请解答：

1. 代码中存在哪些核心问题？
2. 当前实现偏离正确方向的程度如何？
3. 代码中的memoryBarrier有多少是必要的？
4. 针对需求给出具体的改进建议。

内容的提问来源于stack exchange，提问作者pgcomai