Skip to content
Snippets Groups Projects
Commit 0dea4d65 authored by Pedro Gonnet's avatar Pedro Gonnet
Browse files

fix symmetry bug, bad time step.

Former-commit-id: ab05e570d42bc644de825ffbdc6cd50be0c492bb
parent ec55fbbd
No related branches found
No related tags found
No related merge requests found
......@@ -175,7 +175,8 @@ void engine_map_kick_first ( struct cell *c , void *data ) {
int j, k;
struct engine *e = (struct engine *)data;
float pdt, dt_step = e->dt_step, dt = e->dt, hdt = 0.5f*dt;
float dt_min, dt_max, h_max, dx, h2dx_max, dx_max, a[3], v[3], u, u_dt, h, h_dt, rho, v_old[3];
float dt_min, dt_max, h_max, dx, h2dx_max, dx_max;
float a[3], v[3], u, u_dt, h, h_dt, rho, v_old[3];
double x[3], x_old[3];
struct part *restrict p;
struct xpart *restrict xp;
......
......@@ -335,7 +335,7 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
struct cell *finger;
int i, k, redo, count = c->count;
int *pid;
float h, hg, ihg, ihg2, ihg4, h_corr;
float h, hg, ihg, ihg2, ihg4, h_corr, rho, rho_dh, u, fc;
float normDiv_v, normCurl_v;
float dt_step = r->e->dt_step;
TIMER_TIC
......@@ -381,8 +381,8 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
p->density.wcount += kernel_wroot;
/* Final operation on the density. */
p->rho = ihg * ihg2 * ( p->rho + p->mass*kernel_root );
p->rho_dh *= ihg4;
p->rho = rho = ihg * ihg2 * ( p->rho + p->mass*kernel_root );
p->rho_dh = rho_dh = p->rho_dh * ihg4;
/* Compute the smoothing length update (Newton step). */
if ( p->density.wcount_dh == 0.0f )
......@@ -417,20 +417,21 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
}
/* 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;
normDiv_v = fabs( p->density.div_v / 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] ) / rho * ihg4;
/* As of here, particle force variables will be set. Do _NOT_
try to read any particle density variables! */
/* Compute this particle's sound speed. */
p->force.c = sqrtf( const_gamma * ( const_gamma - 1.0f ) * p->u );
u = p->u;
p->force.c = fc = sqrtf( const_gamma * ( const_gamma - 1.0f ) * u );
/* Compute the P/Omega/rho2. */
p->force.POrho2 = p->u * ( const_gamma - 1.0f ) / ( p->rho + hg * p->rho_dh / 3.0f );
p->force.POrho2 = u * ( const_gamma - 1.0f ) / ( rho + hg * rho_dh / 3.0f );
/* Balsara switch */
p->force.balsara = normCurl_v / ( normDiv_v + normCurl_v + 0.0001f * p->force.c * ihg );
p->force.balsara = normCurl_v / ( normDiv_v + normCurl_v + 0.0001f * fc * ihg );
/* Reset the acceleration. */
for ( k = 0 ; k < 3 ; k++ )
......@@ -541,7 +542,7 @@ void runner_dokick2 ( struct runner *r , struct cell *c ) {
}
/* Update the particle's time step. */
p->dt = pdt = const_cfl * h / ( p->force.v_sig );
p->dt = pdt = const_cfl * kernel_gamma * h / ( p->force.v_sig );
/* Update positions and energies at the half-step. */
p->v[0] = v[0] = xp->v_old[0] + hdt * p->a[0];
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment