Skip to content

physics proposed changes #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
lornearnold wants to merge 2 commits into
base: origin
Choose a base branch
from
Open

physics proposed changes #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ecc47ac
physics proposed changes
lornearnold Mar 26, 2024
2fde6bc
Add .DS_Store to .gitignore
lornearnold Mar 26, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
.DS_Store
Cargo.lock
/target
# flatc.exe
42 changes: 23 additions & 19 deletions src/shaders/2D_LOM.wgsl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,40 +23,44 @@ struct Forces {
@group(1) @binding(6) var<storage, read_write> fixity: array<Particle_Settings>;
@group(1) @binding(7) var<storage, read_write> forces: array<Forces>;

const deltaTime: f32 = 0.0000390625;
const deltaTime: f32 = 0.0000390625; // We need to calculate this based on the model parameters
const PI = 3.141592653589793238;

// This calculation comes before the forece-displacement calculation
@compute @workgroup_size(256)
fn main(@builtin(global_invocation_id) global_id: vec3<u32>) {
let id: u32 = global_id.x;

// I think we want to use the equations of motion to update the position and the half step velocity
// I'm using velocities_buf to store the half step velocity

// Note: below, grav must be vec2<f32> with a default value of (0.0, -9.8) m/s/s
// If we made rot the third component of positions, we could use a vec3<f32> for positions and velocities
// and we could use a vec3<f32> for forces. This would allow us to use a single array for all of the data.
// mass and gravity would also need to be vec3<f32> with the third component of mass being the moment of inertia.
velocities_buf[id] = velocities[id] + (forces[id].xy/mass[id] + grav)*deltaTime/2.0;
rot_vel_buf[id] = rot_vel[id] + forces[id].rot*deltaTime/moment_of_inertia[id]/2.0;

if fixity[id].x_vel == 0 {
acc[id] = vec3(vec2((velocities_buf[id] - velocities[id]).x, acc[id].y), acc[id].z);
velocities[id] = vec2(velocities_buf[id].x, velocities[id].y);
} else {
acc[id] = vec3(vec2(0.0, acc[id].y), acc[id].z);
// I don't think the acceleration should be updated here.
velocities_buf[id].x = velocities[id].x;
}

if fixity[id].y_vel == 0 {
acc[id] = vec3(vec2(acc[id].x, (velocities_buf[id] - velocities[id]).y), acc[id].z);
velocities[id] = vec2(velocities[id].x, velocities_buf[id].y);
} else {
acc[id] = vec3(vec2(acc[id].x, 0.0), acc[id].z);
velocities_buf[id].y = velocities[id].y;
}

if fixity[id].rot_vel == 0 {
rot_vel_buf[id] = rot_vel[id];
}

velocities[id] += vec2(forces[id].x, forces[id].y)*deltaTime;

positions[id] = positions[id] + velocities[id] * deltaTime;

if fixity[id].rot_vel == 0 {
acc[id] = vec3(acc[id].xy, rot_vel_buf[id] - rot_vel[id]/deltaTime);
rot_vel[id] = rot_vel_buf[id];
}

rot_vel[id] += forces[id].rot*deltaTime;
rot_vel_buf[id] = rot_vel[id];
rot[id] = (rot[id] + rot_vel[id] * deltaTime)%(2.0*PI);
positions[id] = positions[id] + velocities_buf[id] * deltaTime;
rot[id] = rot[id] + rot_vel_buf[id] * deltaTime;

// I don't think we need to update the forces here.
// Maybe this is intened to be updating externally applied forces?
forces[id].x += forces[id].delX*deltaTime;
forces[id].y += forces[id].delY*deltaTime;
forces[id].rot += forces[id].delRot*deltaTime;
Expand Down
75 changes: 65 additions & 10 deletions src/shaders/2D_Simulation.wgsl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ struct Settings {
friction_coefficient: f32,
rotation: i32,
linear_contact_bonds: i32,
gravity_acc: f32,
gravity_acc: f32, // make this a vector with an x and y acceleration maybe a third component (value = 0) if including rotation
stiffness: f32,
bonds_tear: i32,
bond_force_limit: f32,
Expand Down Expand Up @@ -75,8 +75,13 @@ struct Material {


const deltaTime: f32 = 0.0000390625;
// deltaTime should be a function of the particle mass and stiffness dt = sqrt(min_particle_mass/max_parallel_stiffness)
// It's more complicated than that, but that is a good starting point. At the very least we should check to see what that calc would give us relative to what we're using.
const PI = 3.141592653589793238;

// This calculation comes after the laws of motion calculation
// When we start here, we have the updated position and orientation, but the velocity is one timestep back
// and the velocity_buf is half a time step back.
@compute @workgroup_size(256)
fn main(@builtin(global_invocation_id) global_id: vec3<u32>) {
// data[id*4u] = 1.0;
Expand All @@ -88,12 +93,23 @@ fn main(@builtin(global_invocation_id) global_id: vec3<u32>) {
let mat_id = material_pointers[id];
// let damping: f32 = 0.2; // Damping factor, can be adjusted

let damping: f32 = 0.2; // Damping factor, can be adjusted
// MTIF ("move to input file", see explanation below)
// This is my first time making what I expect will be a recurring comment, so I'll explain it here.
// I think we should keep all variables in a separate file of input parameters.
// This may cost us some flexibility in terms of the simulator being interactive, but we're getting to the point where we're tryin to be scientific in our approach.
// And that involves a fairly rigid, controlled environment anyway.
// We can figure out a way to bring back an interactive mode when the bond physics is working

var net_force = vec2(0.0, 0.0);
var net_moment = 0.0;
// var stress_tensor = vec3(0.0, 0.0, 0.0);

// START of NOT REVIEWED BLOCK
// Bonds
// let bond_shear_lim = 0.5;
//Bonds
let bond_shear_lim = 0.5; // MTIF
var bonded_particles = array<i32, 6u>(-1,-1,-1,-1,-1,-1);
if settings.bonds != 0 {
let start = bond_info[id].x;
Expand All @@ -119,6 +135,8 @@ fn main(@builtin(global_invocation_id) global_id: vec3<u32>) {
let ideal_rot = rot[bond_id];
let rot_disp = rot[id] - rot[bond_id];
net_moment -= (radii[id])*rot_disp/10000.0;
net_moment -= (radii[id])*rot_disp/10000.0; // MTIF

}
} else {
// Linear Bonds, w/ shear resistance
Expand Down Expand Up @@ -230,12 +248,16 @@ fn main(@builtin(global_invocation_id) global_id: vec3<u32>) {
if bonded == false {
contacts[i].bonded = -1;
}
// END of NOT REVIEWED BLOCK
let a = contacts[i].a;
let b = contacts[i].b;
let overlap = max(-distance(a, b), 0.0);
// this is an absolute approach to calculating the normal force.
// nothing wrong with it, but the tangential force needs to be calculated incementally
// I suggest we calculate the normal force incrementally as well.

var normal_stiffness = 10.0;
var shear_stiffness = 0.25;
var normal_stiffness = 10.0; // MTIF
var shear_stiffness = 0.25; // MTIF
if mat_id != -1 {
normal_stiffness = (materials[(material_pointers[b])].normal_stiffness);
shear_stiffness = (materials[(material_pointers[b])].shear_stiffness);
Expand All @@ -244,6 +266,11 @@ fn main(@builtin(global_invocation_id) global_id: vec3<u32>) {
let normal = normalize(positions[a] - positions[b]);
let tangent = vec2(-normal.y, normal.x);

// This looks dangerous to me!
// I think this is duplicating part of the work that the laws of motion should be doing.
// So if it's happening here, too, it might be happening twice.
// Or more like 1.5 times, because the values aren't being stored.
// It's possibly fine, but it's worth checking. Especially if any part of the calculation were to change...
let del_pos_a = velocities[a]*deltaTime;
let del_pos_b = velocities[b]*deltaTime;
let del_rot_a = rot_vel[a]*deltaTime*(radii[a]);//-overlap/2.0);
Expand All @@ -261,6 +288,16 @@ fn main(@builtin(global_invocation_id) global_id: vec3<u32>) {
friction_limit = settings.bond_shear_lim;
moment = false;
}
// I like some things about the code below, but it's only the contact part of the force-disp
// calculation. The bond part is separate. I think we should combine them.
// I think we should have a single function for each "contact model"
// For any given contact, this part of the code should call the appropriate function
// corresponding to the contact model, which will return the forces and moment.
// I think it will only take two functions: one for linear contacts and one for parallel bonded contacts.
// The linear contact can have a bonding ability that defaults to false for unbonded behavior.
// The parallel bond function will have its own force-displacement calculation, call the linear contact
// function (with bonding ability set to false) and then add the forces and moment together.
// With this approach, we only have one place to change the contact forces in the calculation cycle.
contacts[i].tangent_force = contacts[i].tangent_force + rel_tangent*shear_stiffness;//clamp(contacts[i].tangent_force + rel_tangent*shear_stiffness, -friction_limit, friction_limit);
net_force += settings.damping * (normal*normal_force + tangent*contacts[i].tangent_force);
// if moment {
Expand All @@ -285,18 +322,35 @@ fn distance(a: i32, b: i32) -> f32 {
fn store_forces(id: u32, mat_id: i32, net_force: vec2<f32>, net_moment: f32) {
// Apply sum of forces and gravity to velocities
var density = 1.0;
}
}

var density = 1.0; // MTIF
// Move laws of motion to the beginning of the calculation cycle
// Let's rethink this and break it down into its components and put them in the right order.
// Translational motion, then rotational motion
// The outcome will be an updated position and orientation.
// The velocity we use will be a temporary use, half step (i.e., t + dt/2) velocity

// Translational Motion
if mat_id != -1 {
density = materials[mat_id].density;
}
let mass1 = density * PI * radii[id] * radii[id];
let rot_inertia = 0.5*mass1*radii[id]*radii[id];
velocities_buf[id] = velocities[id] + net_force/mass1;
rot_vel_buf[id] = rot_vel[id] + net_moment/rot_inertia;
let mass1 = density * PI * radii[id] * radii[id]; // make it a function mass(density,radius)
//let half_step_velocity = velocities[id] + 0.5 * (resultant_force/mass1 - grav) * deltaTime
// Rotational Motion
let rot_inertia = 0.5 * mass1 * radii[id] * radii[id]; // make it a function rot_inertia(density,radius)
//let half_step_ang_vel = rot_vel[id] + 0.5 * net_moment/rot_inertia) * deltaTime;

// Not sure these belong here. They seem to be part of the laws of motion calculation.
// But they're not entirely consistent with the other set of laws of motion in 2D_LOM.wgsl,
// so something's going to be inconsistent.
rot_vel_buf[id] =
if settings.gravity == 1 && settings.planet_mode == 1 {
let delta = (vec2(0.0, 0.0) - positions[id]);
velocities_buf[id] += delta/length(delta) * 9.81 * settings.gravity_acc * deltaTime;
} else if settings.gravity == 1 {
let gravity = 9.81 * settings.gravity_acc * deltaTime;
let gravity = 9.81 * settings.gravity_acc * deltaTime; // MTIF (9.81)
velocities_buf[id] += vec2(0.0, -gravity);
}
}
Expand All @@ -305,8 +359,8 @@ fn walls(id: u32) {
// BS Walls
let pos = positions[id];
let rad = radii[id];
let elasticity = 0.5;
let anti_stick_coating = 0.01;
let elasticity = 0.5; // MTIF
let anti_stick_coating = 0.01; // MTIF
let yH = settings.vert_bound;
let xW = settings.hor_bound;

Expand All @@ -333,4 +387,5 @@ fn walls(id: u32) {
// fn stress_tensor(id: u32, force: vec2<f32>, delta: vec2<f32>) -> vec3<f32> {
// var tensor = vec3(0.0, 0.0, 0.0);
// let

// }