Commit e6a3af88 authored by Johannes Blaschke's avatar Johannes Blaschke
Browse files

working on vparticle collision rules

parent eec2fc80
......@@ -141,7 +141,7 @@ send_particles_to_right_cpu(particle_state * part, int * delta_n_particles, int
// return vdot(& delta_r, & p->del_p);
// }
void
void
sort_particles_into_cells(vector3 * gshift, particle * particles, int nparticles, cell * cells)
{
particle * firstp = & particles[0];
......@@ -163,7 +163,7 @@ sort_particles_into_cells(vector3 * gshift, particle * particles, int nparticles
#endif
}
}
void
void
sort_vparticles_into_cells(vector3 * gshift, vparticle * vparticles, int nvparticles, cell * cells)
{
vparticle * firstvp = & vparticles[0];
......@@ -176,7 +176,7 @@ sort_vparticles_into_cells(vector3 * gshift, vparticle * vparticles, int nvparti
c->nparticles++;
c->nvparticles++;
// #ifdef ENABLE_PFIELD
// c->div_v += particle_div_v_vp(vp, c);
// #endif
......@@ -311,7 +311,7 @@ calculate_random_velocity_sum(vector3 * gshift, particle * particles, int nparti
}
#ifdef ENABLE_PFIELD
void grid_assign_grad_pressure(world_state * world)
void grid_assign_grad_pressure(world_state * world)
{
particle * particles = get_particles(world);
int nparticles = get_nparticles(world);
......@@ -343,7 +343,7 @@ void grid_assign_grad_pressure(world_state * world)
for(int i = 0; i < cpu.numbercells; i++){
cell * c_ptr = & world->cells[i];
int n_avg = config.pressure_nhist;
vcpy(& c_ptr->hist_grad_p[c_ptr->index_hist_grad_p], & c_ptr->grad_pressure);
c_ptr->index_hist_grad_p ++;
if(c_ptr->index_hist_grad_p >= n_avg) c_ptr->index_hist_grad_p = 0;
......@@ -362,7 +362,7 @@ void grid_assign_grad_pressure(world_state * world)
}
void grid_diffuse_pressure(world_state * world, int last) {
static int * n_tracers_in;
if(n_tracers_in == NULL) n_tracers_in = calloc(cpu.numbercells, sizeof(int));
static double * cell__tmp_in;
......@@ -385,14 +385,14 @@ void grid_diffuse_pressure(world_state * world, int last) {
// assign pressure before bounceback
// note: bounceback can alter pressure value
p->p = 0.01 * c->pressure;
p->p = 0.01 * c->pressure;
generate_random_v(& p->v, 0.1);
vapply_velocity_periodic(& p->r, & p->v, 1.0);
if(config.boundaries[0] + config.boundaries[1] + config.boundaries[2] > 0)
check_if_particle_hit_wall(p, 0.5, config.N.x, config.N.y-1, config.N.z-1);
resolve_particle_collisions_with_obstacles(p, collision_check_trypanosome, 1.0, world, 0);
resolve_particle_collisions_with_obstacles(p, collision_check_squirmers, 1.0, world, 0);
......@@ -422,7 +422,7 @@ void grid_diffuse_pressure(world_state * world, int last) {
c_ptr->delta_pressure += c_ptr->pressure - old_pressure;
}
if(last) {
if(last) {
free(n_tracers_in);
free(n_tracers_out);
free(cell__tmp_in);
......@@ -460,8 +460,8 @@ void grid_update_local_pressure(world_state * world)
double total_weight = 0;
// in case simulations haven't caught up:
if(world->timestep < n_avg) n_avg = world->timestep + 1;
if(world->timestep < n_avg) n_avg = world->timestep + 1;
for(int m=0; m<n_avg; m++){
avg_nparticles += c_ptr->hist_nparticles[m];
// avg_div_v += c_ptr->hist_nparticles[m] * c_ptr->hist_div_v[m];
......@@ -473,7 +473,7 @@ void grid_update_local_pressure(world_state * world)
// if(total_weight > 0)
// avg_div_v /= total_weight;
// else
// else
// avg_div_v = 0;
// ***** UPDATE pressure and delta_pressure ***********************************************************************
......@@ -554,7 +554,7 @@ conserve_momentum_vp(vparticle * vparticles, int nvparticles, squirmer * squirme
if((config.squirmeron == 1) && (config.without_squirmer_stream == 0))
clear_momentumsums(msums, config.nsquirmers);
// ******* ITTERATE OVER VPARTICLES *************************************************************
vparticle * firstvp = & vparticles[0];
......@@ -693,7 +693,7 @@ conserve_momentum(particle* particles, int nparticles, cell* cells)
// particle * p = & part->p[i];
//
// if(!p->active) continue;
//
//
// int cellnumber = access_cells(p->cell.x, p->cell.y, p->cell.z);
// cell * c = & cells[cellnumber];
// // if(c->p_thermostat_done == 1) continue;
......@@ -702,13 +702,13 @@ conserve_momentum(particle* particles, int nparticles, cell* cells)
//
// if(delta_n < 0){
// double coin = generate_random_number_in_interval(0.0, 1.0);
//
//
// if(coin < exp(20.*delta_n/config.ppercell)) continue;
//
// vector3 dv; vcpy(& dv, & p->v);
// vinverse(& dv);
// vadd_ip(& c->averagev, & dv);
//
//
// if(config.angular)
// {
// vcross_add_ip(& c->sum_l, & p->r, & dv);
......@@ -761,7 +761,7 @@ conserve_momentum(particle* particles, int nparticles, cell* cells)
// coordinate_valid = 1;
// break;
// }
//
//
// if(coordinate_valid == 1){
// // fprintf(stdout, "%f, %f, %f => %f, %f, %f => %f, %f, %f\n",rmin.x, rmin.y, rmin.z, r_rand.x, r_rand.y, r_rand.z, rmax.x, rmax.y, rmax.z);
//
......@@ -803,7 +803,7 @@ void undo_gridshift(vector3 * shift, world_state * world)
cpu.upleftlimit.z += shift->z;
cpu.downrightlimit.y += shift->y;
cpu.downrightlimit.z += shift->z;
cell__clear_obstacle_cache(world);
//vzero(shift); // Just for testing purposes
......
......@@ -306,7 +306,7 @@ generate_cell_ghost_particles(vector3i cell_index, vparticle_state * vparticles,
trypanosome_node * closest_node = & int_facet.nodes[0];
vector3 * r_slice_centre = & closest_node->r_center;
vector3 * r_closest_node = closest_node->r;
vector3 * r_closest_node = closest_node->r;
vector3 * v_closest_node = closest_node->v;
vector3 * v_slice_centre = & closest_node->v_center;
// vector3 * slice_del_p = closest_node->del_p_center;
......@@ -317,12 +317,13 @@ generate_cell_ghost_particles(vector3i cell_index, vparticle_state * vparticles,
vp->mass_obstacle_node = 4*closest_node->mass; // 4 nodes per facet
// vp->del_p_obstacle_node = closest_node.del_p;
vp->del_p_obstacle_node = int_facet.facet_del_p;
vp->del_p_obstacle = & world->trypanosome->del_p; //TODO: HAXOR fix asap (incompatible with multiple tryps)
// vp->del_p_obstacle = int_facet.facet_del_p;
//TODO: HAXOR fix asap (incompatible with multiple tryps)
vp->del_p_obstacle = & world->trypanosome->del_p;
// vp->del_p_obstacle = int_facet.facet_del_p;
// save slice data
// vcpy(& vp->r_obstacle , r_slice_centre);
vcpy(& vp->r_obstacle, r_closest_node);
vcpy(& vp->r_obstacle, r_closest_node);
// vp->del_p_obstacle = slice_del_p;
// mpcd
......@@ -331,12 +332,12 @@ generate_cell_ghost_particles(vector3i cell_index, vparticle_state * vparticles,
// keep_vector3_in_box(& vp->r, & config.N);
vzero(& vp->del_p);
generate_random_v(& vp->del_p, sigma);
// vadd_ip(& vp->del_p, v_slice_centre);
vadd_ip(& vp->del_p, v_slice_centre);
// vadd_ip(& vp->del_p, v_closest_node);
vector3 avg_v; vzero(& avg_v);
for(int i=1; i<5; i++)
vmuladd_ip(& avg_v, 0.25, int_facet.nodes[i].v);
vadd_ip(& vp->del_p, & avg_v);
// vector3 avg_v; vzero(& avg_v);
// for(int i=1; i<5; i++)
// vmuladd_ip(& avg_v, 0.25, int_facet.nodes[i].v);
// vadd_ip(& vp->del_p, & avg_v);
// save metadata
vp->belonging=belonging;
......
......@@ -193,7 +193,7 @@ void stream_particles(world_state * world, double dt, const int timestep, const
for(squirmer * s = get_first_s(firsts, lasts); s!=lasts; s=get_next_s(s, lasts)){
vzero(& s->pressure_f);
}
firsts = & world->osquirmers.s[0];
lasts = & world->osquirmers.s[world->osquirmers.index];
for(squirmer * s = get_first_s(firsts, lasts); s!=lasts; s=get_next_s(s, lasts)){
......@@ -258,7 +258,7 @@ void stream_particles(world_state * world, double dt, const int timestep, const
if(config.boundaries[0] + config.boundaries[1] + config.boundaries[2] > 0)
check_if_particle_hit_wall(p, dt/2, config.N.x, config.N.y-1, config.N.z-1);
resolve_particle_collisions_with_obstacles(p, collision_check_trypanosome, dt, world, 1); //HAXOR, TODO find more elegant way to do this (perhaps get rid of this coll_check);
resolve_particle_collisions_with_obstacles(p, collision_check_squirmers, dt, world, 1); //'tis complicated why it should be this order....
......
......@@ -565,7 +565,7 @@ trypanosome_time_driven_md(world_state * world, const int n_md_step, const doubl
{ // verlet algorithm
trypanosome_network__md_update_r(net, dt);
trypanosome_network__md_update_f(net, md_time, vertex_buffer, end_node_buffer,
kw_flagellum, drag_force, gamma);
kw_flagellum, drag_force, gamma);
// assign forces due to boundaries *********************************************************************
if(config.boundaries[0] + config.boundaries[1] + config.boundaries[2] > 0)
......
......@@ -1341,7 +1341,7 @@ vertex__md_update_v(vertex * node, double dt, vertex * old_vertex)
// total force acting on node **************************************************************************
vector3 f_total;
vector3 * f = & node->f;
vector3 * f = & node->f;
vector3 * f_ext = & node->f_ext;
vadd(& f_total, f, f_ext);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment