Commit e804a9bb authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Correct calculation of rot_v and div_v in default case.

parent 967a0af8
...@@ -72,9 +72,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part( ...@@ -72,9 +72,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->rho = 0.f; p->rho = 0.f;
p->rho_dh = 0.f; p->rho_dh = 0.f;
p->density.div_v = 0.f; p->density.div_v = 0.f;
p->density.curl_v[0] = 0.f; p->density.rot_v[0] = 0.f;
p->density.curl_v[1] = 0.f; p->density.rot_v[1] = 0.f;
p->density.curl_v[2] = 0.f; p->density.rot_v[2] = 0.f;
} }
/** /**
...@@ -110,6 +110,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( ...@@ -110,6 +110,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Compute the derivative term */ /* Compute the derivative term */
p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho); p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho);
/* Finish calculation of the velocity curl components */
p->density.rot_v[0] *= ih4 * irho;
p->density.rot_v[1] *= ih4 * irho;
p->density.rot_v[2] *= ih4 * irho;
/* Finish calculation of the velocity divergence */
p->density.div_v *= ih4 * irho;
} }
/** /**
...@@ -133,9 +141,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -133,9 +141,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Pre-compute some stuff for the balsara switch. */ /* Pre-compute some stuff for the balsara switch. */
const float normDiv_v = fabs(p->density.div_v / p->rho * ih4); const float normDiv_v = fabs(p->density.div_v / p->rho * ih4);
const float normCurl_v = sqrtf(p->density.curl_v[0] * p->density.curl_v[0] + const float normRot_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
p->density.curl_v[1] * p->density.curl_v[1] + p->density.rot_v[1] * p->density.rot_v[1] +
p->density.curl_v[2] * p->density.curl_v[2]) / p->density.rot_v[2] * p->density.rot_v[2]) /
p->rho * ih4; p->rho * ih4;
/* Compute this particle's sound speed. */ /* Compute this particle's sound speed. */
...@@ -147,7 +155,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -147,7 +155,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->force.POrho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega); p->force.POrho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega);
/* Balsara switch */ /* Balsara switch */
p->force.balsara = normDiv_v / (normDiv_v + normCurl_v + 0.0001f * fc * ih); p->force.balsara = normDiv_v / (normDiv_v + normRot_v + 0.0001f * fc * ih);
/* Viscosity parameter decay time */ /* Viscosity parameter decay time */
const float tau = h / (2.f * const_viscosity_length * p->force.c); const float tau = h / (2.f * const_viscosity_length * p->force.c);
......
...@@ -80,7 +80,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( ...@@ -80,7 +80,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pi->density.wcount_dh -= xi * wi_dx; pi->density.wcount_dh -= xi * wi_dx;
pi->density.div_v += mj * dvdr * wi_dx; pi->density.div_v += mj * dvdr * wi_dx;
for (k = 0; k < 3; k++) pi->density.curl_v[k] += mj * curlvr[k] * wi_dx; for (k = 0; k < 3; k++) pi->density.rot_v[k] += mj * curlvr[k] * wi_dx;
/* Compute density of pj. */ /* Compute density of pj. */
h_inv = 1.0 / hj; h_inv = 1.0 / hj;
...@@ -93,7 +93,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( ...@@ -93,7 +93,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pj->density.wcount_dh -= xj * wj_dx; pj->density.wcount_dh -= xj * wj_dx;
pj->density.div_v += mi * dvdr * wj_dx; pj->density.div_v += mi * dvdr * wj_dx;
for (k = 0; k < 3; k++) pj->density.curl_v[k] += mi * curlvr[k] * wj_dx; for (k = 0; k < 3; k++) pj->density.rot_v[k] += mi * curlvr[k] * wj_dx;
} }
/** /**
...@@ -111,7 +111,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density( ...@@ -111,7 +111,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
vector dx[3], dv[3]; vector dx[3], dv[3];
vector vi[3], vj[3]; vector vi[3], vj[3];
vector dvdr, div_vi, div_vj; vector dvdr, div_vi, div_vj;
vector curlvr[3], curl_vi[3], curl_vj[3]; vector curlvr[3], rot_vi[3], rot_vj[3];
int k, j; int k, j;
#if VEC_SIZE == 8 #if VEC_SIZE == 8
...@@ -178,14 +178,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density( ...@@ -178,14 +178,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
wcounti.v = wi.v; wcounti.v = wi.v;
wcounti_dh.v = xi.v * wi_dx.v; wcounti_dh.v = xi.v * wi_dx.v;
div_vi.v = mj.v * dvdr.v * wi_dx.v; div_vi.v = mj.v * dvdr.v * wi_dx.v;
for (k = 0; k < 3; k++) curl_vi[k].v = mj.v * curlvr[k].v * wi_dx.v; for (k = 0; k < 3; k++) rot_vi[k].v = mj.v * curlvr[k].v * wi_dx.v;
rhoj.v = mi.v * wj.v; rhoj.v = mi.v * wj.v;
rhoj_dh.v = mi.v * (vec_set1(3.0f) * wj.v + xj.v * wj_dx.v); rhoj_dh.v = mi.v * (vec_set1(3.0f) * wj.v + xj.v * wj_dx.v);
wcountj.v = wj.v; wcountj.v = wj.v;
wcountj_dh.v = xj.v * wj_dx.v; wcountj_dh.v = xj.v * wj_dx.v;
div_vj.v = mi.v * dvdr.v * wj_dx.v; div_vj.v = mi.v * dvdr.v * wj_dx.v;
for (k = 0; k < 3; k++) curl_vj[k].v = mi.v * curlvr[k].v * wj_dx.v; for (k = 0; k < 3; k++) rot_vj[k].v = mi.v * curlvr[k].v * wj_dx.v;
for (k = 0; k < VEC_SIZE; k++) { for (k = 0; k < VEC_SIZE; k++) {
pi[k]->rho += rhoi.f[k]; pi[k]->rho += rhoi.f[k];
...@@ -193,13 +193,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density( ...@@ -193,13 +193,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
pi[k]->density.wcount += wcounti.f[k]; pi[k]->density.wcount += wcounti.f[k];
pi[k]->density.wcount_dh -= wcounti_dh.f[k]; pi[k]->density.wcount_dh -= wcounti_dh.f[k];
pi[k]->density.div_v += div_vi.f[k]; pi[k]->density.div_v += div_vi.f[k];
for (j = 0; j < 3; j++) pi[k]->density.curl_v[j] += curl_vi[j].f[k]; for (j = 0; j < 3; j++) pi[k]->density.rot_v[j] += rot_vi[j].f[k];
pj[k]->rho += rhoj.f[k]; pj[k]->rho += rhoj.f[k];
pj[k]->rho_dh -= rhoj_dh.f[k]; pj[k]->rho_dh -= rhoj_dh.f[k];
pj[k]->density.wcount += wcountj.f[k]; pj[k]->density.wcount += wcountj.f[k];
pj[k]->density.wcount_dh -= wcountj_dh.f[k]; pj[k]->density.wcount_dh -= wcountj_dh.f[k];
pj[k]->density.div_v += div_vj.f[k]; pj[k]->density.div_v += div_vj.f[k];
for (j = 0; j < 3; j++) pj[k]->density.curl_v[j] += curl_vj[j].f[k]; for (j = 0; j < 3; j++) pj[k]->density.rot_v[j] += rot_vj[j].f[k];
} }
#else #else
...@@ -255,7 +255,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( ...@@ -255,7 +255,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
pi->density.wcount_dh -= xi * wi_dx; pi->density.wcount_dh -= xi * wi_dx;
pi->density.div_v += mj * dvdr * wi_dx; pi->density.div_v += mj * dvdr * wi_dx;
for (k = 0; k < 3; k++) pi->density.curl_v[k] += mj * curlvr[k] * wi_dx; for (k = 0; k < 3; k++) pi->density.rot_v[k] += mj * curlvr[k] * wi_dx;
} }
/** /**
...@@ -273,7 +273,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj, ...@@ -273,7 +273,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
vector dx[3], dv[3]; vector dx[3], dv[3];
vector vi[3], vj[3]; vector vi[3], vj[3];
vector dvdr; vector dvdr;
vector curlvr[3], curl_vi[3]; vector curlvr[3], rot_vi[3];
int k, j; int k, j;
#if VEC_SIZE == 8 #if VEC_SIZE == 8
...@@ -331,7 +331,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj, ...@@ -331,7 +331,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
wcounti.v = wi.v; wcounti.v = wi.v;
wcounti_dh.v = xi.v * wi_dx.v; wcounti_dh.v = xi.v * wi_dx.v;
div_vi.v = mj.v * dvdr.v * wi_dx.v; div_vi.v = mj.v * dvdr.v * wi_dx.v;
for (k = 0; k < 3; k++) curl_vi[k].v = mj.v * curlvr[k].v * wi_dx.v; for (k = 0; k < 3; k++) rot_vi[k].v = mj.v * curlvr[k].v * wi_dx.v;
for (k = 0; k < VEC_SIZE; k++) { for (k = 0; k < VEC_SIZE; k++) {
pi[k]->rho += rhoi.f[k]; pi[k]->rho += rhoi.f[k];
...@@ -339,7 +339,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj, ...@@ -339,7 +339,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
pi[k]->density.wcount += wcounti.f[k]; pi[k]->density.wcount += wcounti.f[k];
pi[k]->density.wcount_dh -= wcounti_dh.f[k]; pi[k]->density.wcount_dh -= wcounti_dh.f[k];
pi[k]->density.div_v += div_vi.f[k]; pi[k]->density.div_v += div_vi.f[k];
for (j = 0; j < 3; j++) pi[k]->density.curl_v[j] += curl_vi[j].f[k]; for (j = 0; j < 3; j++) pi[k]->density.rot_v[j] += rot_vi[j].f[k];
} }
#else #else
......
...@@ -80,7 +80,7 @@ struct part { ...@@ -80,7 +80,7 @@ struct part {
float wcount_dh; float wcount_dh;
/* Particle velocity curl. */ /* Particle velocity curl. */
float curl_v[3]; float rot_v[3];
/* Particle number density. */ /* Particle number density. */
float wcount; float wcount;
......
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