Commit 58d81f85 authored by Pedro Gonnet's avatar Pedro Gonnet
Browse files

squeeze particles a bit, i.e. union with overlapping ephemeric data for...

squeeze particles a bit, i.e. union with overlapping ephemeric data for density and force computations.


Former-commit-id: 0d1cc866c157a776a0d4e19266097125508d2946
parent 32870962
......@@ -37,12 +37,12 @@ void printParticle ( struct part *parts , long long int id ) {
parts[i].v[0], parts[i].v[1], parts[i].v[2],
parts[i].a[0], parts[i].a[1], parts[i].a[2],
parts[i].h,
parts[i].h_dt,
parts[i].force.h_dt,
parts[i].wcount,
parts[i].mass,
parts[i].rho,
parts[i].u,
parts[i].u_dt,
parts[i].force.u_dt,
parts[i].dt
);
}
......
......@@ -96,12 +96,12 @@ void engine_prepare ( struct engine *e ) {
for ( k = 0 ; k < s->nr_parts ; k++ )
if ( s->parts[k].dt <= dt_step ) {
s->parts[k].wcount = 0.0f;
s->parts[k].wcount_dh = 0.0f;
s->parts[k].density.wcount_dh = 0.0f;
s->parts[k].rho = 0.0f;
s->parts[k].rho_dh = 0.0f;
s->parts[k].div_v = 0.0f;
for ( j = 0 ; j < 3 ; ++j)
s->parts[k].curl_v[j] = 0.0f;
s->parts[k].density.div_v = 0.0f;
for ( j = 0 ; j < 3 ; ++j)
s->parts[k].density.curl_v[j] = 0.0f;
}
// printf( "engine_prepare: re-setting particle data took %.3f ms.\n" , (double)(getticks() - tic) / CPU_TPS * 1000 );
......@@ -212,7 +212,7 @@ void engine_step ( struct engine *e , int sort_queues ) {
xp->v_old[0] = p->v[0] + hdt * p->a[0];
xp->v_old[1] = p->v[1] + hdt * p->a[1];
xp->v_old[2] = p->v[2] + hdt * p->a[2];
xp->u_old = p->u + hdt * p->u_dt;
xp->u_old = p->u + hdt * p->force.u_dt;
/* Move the particles with the velocitie at the half-step. */
p->x[0] += dt * xp->v_old[0];
......@@ -223,13 +223,13 @@ void engine_step ( struct engine *e , int sort_queues ) {
p->v[0] += dt * p->a[0];
p->v[1] += dt * p->a[1];
p->v[2] += dt * p->a[2];
p->u *= expf( p->u_dt / p->u * dt );
p->h *= expf( p->h_dt / p->h * dt );
p->u *= expf( p->force.u_dt / p->u * dt );
p->h *= expf( p->force.h_dt / p->h * dt );
/* Integrate other values if this particle will not be updated. */
if ( p->dt > dt_step ) {
p->rho *= expf( -3.0f * p->h_dt / p->h * dt );
p->POrho2 = p->u * ( const_gamma - 1.0f ) / ( p->rho + p->h * p->rho_dh / 3.0f );
p->rho *= expf( -3.0f * p->force.h_dt / p->h * dt );
p->force.POrho2 = p->u * ( const_gamma - 1.0f ) / ( p->rho + p->h * p->rho_dh / 3.0f );
}
}
......@@ -290,18 +290,18 @@ void engine_step ( struct engine *e , int sort_queues ) {
/* Scale the derivatives if they're freshly computed. */
if ( p->dt <= dt_step ) {
p->h_dt *= p->h * 0.333333333f;
p->force.h_dt *= p->h * 0.333333333f;
lcount += 1;
}
/* Update the particle's time step. */
p->dt = const_cfl * p->h / ( p->v_sig );
p->dt = const_cfl * p->h / ( p->force.v_sig );
/* Update positions and energies at the half-step. */
p->v[0] = xp->v_old[0] + hdt * p->a[0];
p->v[1] = xp->v_old[1] + hdt * p->a[1];
p->v[2] = xp->v_old[2] + hdt * p->a[2];
p->u = xp->u_old + hdt * p->u_dt;
p->u = xp->u_old + hdt * p->force.u_dt;
/* Get the smallest/largest dt. */
ldt_min = fminf( ldt_min , p->dt );
......
......@@ -58,7 +58,7 @@ struct xpart {
struct part {
/* Particle velocity. */
float v[3] __attribute__((aligned (16)));
float v[3];
/* Particle mass. */
float mass;
......@@ -66,43 +66,57 @@ struct part {
/* Particle density. */
float rho;
/* Particle velocity divergence. */
float div_v;
/* Particle velocity curl. */
float curl_v[3] __attribute__((aligned (16)));
/* Balsara switch */
float balsara;
/* Particle pressure. */
// float P;
/* Derivative of the density with respect to this particle's smoothing length. */
float rho_dh;
/* Aggregate quantities. */
float POrho2;
/* Store density/force specific stuff. */
union {
/* Change in particle energy over time. */
float u_dt;
struct {
/* Particle velocity divergence. */
float div_v;
/* Derivative of particle number density. */
float wcount_dh;
/* Particle velocity curl. */
float curl_v[3];
} density;
struct {
/* Balsara switch */
float balsara;
/* Aggregate quantities. */
float POrho2;
/* Change in smoothing length over time. */
float h_dt;
/* Change in particle energy over time. */
float u_dt;
/* Change in smoothing length over time. */
float h_dt;
/* Particle acceleration. */
float a[3] __attribute__((aligned (16)));
/* Signal velocity */
float v_sig;
/* Sound speed */
float c;
/* Sound speed */
float c;
/* Signal velocity */
float v_sig;
} force;
};
/* Particle pressure. */
// float P;
/* Derivative of the density with respect to this particle's smoothing length. */
float rho_dh;
/* Particle acceleration. */
float a[3];
/* Particle number density. */
// int icount;
float wcount;
float wcount_dh;
/* Particle internal energy. */
float u;
......
......@@ -330,7 +330,7 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
struct cell *finger;
int i, k, redo, count = c->count;
int *pid;
float ihg, ihg2, h_corr;
float ihg, ihg2, ihg4, h_corr;
float normDiv_v, normCurl_v;
float dt_step = r->e->dt_step;
TIMER_TIC
......@@ -365,15 +365,16 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
/* Is this part within the timestep? */
if ( cp->dt <= dt_step ) {
/* Some smoothing length multiples*/
/* Some smoothing length multiples. */
ihg = kernel_igamma / p->h;
ihg2 = ihg * ihg;
ihg4 = ihg2 * ihg2;
/* Adjust the computed weighted number of neighbours. */
p->wcount += kernel_wroot;
/* Compute the smoothing length update (Newton step). */
h_corr = ( const_nwneigh - p->wcount ) / p->wcount_dh;
h_corr = ( const_nwneigh - p->wcount ) / p->density.wcount_dh;
/* Truncate to the range [ -p->h/2 , p->h ]. */
h_corr = fminf( h_corr , p->h );
......@@ -387,52 +388,44 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
if ( p->wcount > const_nwneigh + const_delta_nwneigh ||
p->wcount < const_nwneigh - const_delta_nwneigh ) {
// printf( "runner_doghost: particle %lli (h=%e,h_dt=%e,depth=%i) has bad wcount=%.3f.\n" , p->id , p->h , p->h_dt , c->depth , p->wcount ); fflush(stdout);
// p->h += ( p->wcount + kernel_root - const_nwneigh ) / p->wcount_dh;
// p->h += ( p->wcount + kernel_root - const_nwneigh ) / p->density.wcount_dh;
pid[redo] = pid[i];
redo += 1;
p->wcount = 0.0;
p->wcount_dh = 0.0;
p->density.wcount_dh = 0.0;
p->rho = 0.0;
p->rho_dh = 0.0;
p->div_v = 0.0;
p->density.div_v = 0.0;
for ( k=0 ; k < 3 ; k++)
p->curl_v[k] = 0.0;
p->density.curl_v[k] = 0.0;
continue;
}
/* Final operation on the density */
/* Final operation on the density. */
p->rho = ihg * ihg2 * ( p->rho + p->mass*kernel_root );
p->rho_dh *= ihg2 * ihg2;
p->rho_dh *= ihg4;
/* Pre-compute some stuff for the balsara switch. */
normDiv_v = fabs( p->density.div_v / p->rho * ihg4 );
normCurl_v = sqrtf( p->density.curl_v[0] * p->density.curl_v[0] + p->density.curl_v[1] * p->density.curl_v[1] + p->density.curl_v[2] * p->density.curl_v[2] ) / p->rho * ihg4;
/* Compute this particle's sound speed. */
p->c = sqrtf( const_gamma * ( const_gamma - 1.0f ) * p->u );
p->force.c = sqrtf( const_gamma * ( const_gamma - 1.0f ) * p->u );
/* Compute the P/Omega/rho2. */
p->POrho2 = p->u * ( const_gamma - 1.0f ) / ( p->rho + p->h * p->rho_dh / 3.0f );
/* Final operation on the velocity divergence */
p->div_v /= p->rho;
p->div_v *= ihg2 * ihg2;
/* Final operation on the velocity curl */
for ( k=0 ; k < 3 ; k++ ){
p->curl_v[k] /= -p->rho;
p->curl_v[k] *= ihg2 * ihg2;
}
p->force.POrho2 = p->u * ( const_gamma - 1.0f ) / ( p->rho + p->h * p->rho_dh / 3.0f );
/* Balsara switch */
normDiv_v = fabs( p->div_v );
normCurl_v = sqrtf( p->curl_v[0] * p->curl_v[0] + p->curl_v[1] * p->curl_v[1] + p->curl_v[2] * p->curl_v[2] );
p->balsara = normCurl_v / ( normDiv_v + normCurl_v + 0.0001f * p->c / p->h );
p->force.balsara = normCurl_v / ( normDiv_v + normCurl_v + 0.0001f * p->force.c / p->h );
/* Reset the acceleration. */
for ( k = 0 ; k < 3 ; k++ )
p->a[k] = 0.0f;
/* Reset the time derivatives. */
p->u_dt = 0.0f;
p->h_dt = 0.0f;
p->v_sig = 0.0f;
p->force.u_dt = 0.0f;
p->force.h_dt = 0.0f;
p->force.v_sig = 0.0f;
}
......
This diff is collapsed.
Supports Markdown
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