diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h index b832333377f37e33a99eb0b9fdb7cb15d1233d0b..850ac68352dcd38426cdef36442ede93eda2270b 100644 --- a/src/hydro/Default/hydro.h +++ b/src/hydro/Default/hydro.h @@ -72,9 +72,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part( p->rho = 0.f; p->rho_dh = 0.f; p->density.div_v = 0.f; - p->density.curl_v[0] = 0.f; - p->density.curl_v[1] = 0.f; - p->density.curl_v[2] = 0.f; + p->density.rot_v[0] = 0.f; + p->density.rot_v[1] = 0.f; + p->density.rot_v[2] = 0.f; } /** @@ -110,6 +110,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( /* Compute the derivative term */ 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( /* Pre-compute some stuff for the balsara switch. */ 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] + - p->density.curl_v[1] * p->density.curl_v[1] + - p->density.curl_v[2] * p->density.curl_v[2]) / + const float normRot_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] + + p->density.rot_v[1] * p->density.rot_v[1] + + p->density.rot_v[2] * p->density.rot_v[2]) / p->rho * ih4; /* Compute this particle's sound speed. */ @@ -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); /* 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 */ const float tau = h / (2.f * const_viscosity_length * p->force.c); diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h index e67d6422a8a3019beee7942f3381fabe5951640c..db32a9a7b54177e73c514d7cebd1ed12acae30ce 100644 --- a/src/hydro/Default/hydro_iact.h +++ b/src/hydro/Default/hydro_iact.h @@ -80,7 +80,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( pi->density.wcount_dh -= xi * 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. */ h_inv = 1.0 / hj; @@ -93,7 +93,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( pj->density.wcount_dh -= xj * 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( vector dx[3], dv[3]; vector vi[3], vj[3]; 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; #if VEC_SIZE == 8 @@ -178,14 +178,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density( wcounti.v = wi.v; wcounti_dh.v = xi.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_dh.v = mi.v * (vec_set1(3.0f) * wj.v + xj.v * wj_dx.v); wcountj.v = wj.v; wcountj_dh.v = xj.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++) { pi[k]->rho += rhoi.f[k]; @@ -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_dh -= wcounti_dh.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_dh -= rhoj_dh.f[k]; pj[k]->density.wcount += wcountj.f[k]; pj[k]->density.wcount_dh -= wcountj_dh.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 @@ -255,7 +255,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( pi->density.wcount_dh -= xi * 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, vector dx[3], dv[3]; vector vi[3], vj[3]; vector dvdr; - vector curlvr[3], curl_vi[3]; + vector curlvr[3], rot_vi[3]; int k, j; #if VEC_SIZE == 8 @@ -331,7 +331,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj, wcounti.v = wi.v; wcounti_dh.v = xi.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++) { pi[k]->rho += rhoi.f[k]; @@ -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_dh -= wcounti_dh.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 diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h index a4096d5c30d525307d5327559ef6b007c6931486..946145587262a1502f9e2ee847f10500c1a6b986 100644 --- a/src/hydro/Default/hydro_part.h +++ b/src/hydro/Default/hydro_part.h @@ -80,7 +80,7 @@ struct part { float wcount_dh; /* Particle velocity curl. */ - float curl_v[3]; + float rot_v[3]; /* Particle number density. */ float wcount;