Commit fa504958 authored by James Willis's avatar James Willis
Browse files

Merge branch 'gadget2-part-update' into intrinsic-vectorisation

parents 77444460 46ea1a8b
...@@ -638,11 +638,11 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset) { ...@@ -638,11 +638,11 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset) {
*/ */
void cell_init_parts(struct cell *c, void *data) { void cell_init_parts(struct cell *c, void *data) {
struct part *p = c->parts; struct part *restrict p = c->parts;
struct xpart *xp = c->xparts; struct xpart *restrict xp = c->xparts;
const int count = c->count; const size_t count = c->count;
for (int i = 0; i < count; ++i) { for (size_t i = 0; i < count; ++i) {
p[i].ti_begin = 0; p[i].ti_begin = 0;
p[i].ti_end = 0; p[i].ti_end = 0;
xp[i].v_full[0] = p[i].v[0]; xp[i].v_full[0] = p[i].v[0];
...@@ -665,13 +665,14 @@ void cell_init_parts(struct cell *c, void *data) { ...@@ -665,13 +665,14 @@ void cell_init_parts(struct cell *c, void *data) {
*/ */
void cell_init_gparts(struct cell *c, void *data) { void cell_init_gparts(struct cell *c, void *data) {
struct gpart *gp = c->gparts; struct gpart *restrict gp = c->gparts;
const int gcount = c->gcount; const size_t gcount = c->gcount;
for (int i = 0; i < gcount; ++i) { for (size_t i = 0; i < gcount; ++i) {
gp[i].ti_begin = 0; gp[i].ti_begin = 0;
gp[i].ti_end = 0; gp[i].ti_end = 0;
gravity_first_init_gpart(&gp[i]); gravity_first_init_gpart(&gp[i]);
gravity_init_gpart(&gp[i]);
} }
c->ti_end_min = 0; c->ti_end_min = 0;
c->ti_end_max = 0; c->ti_end_max = 0;
......
...@@ -82,7 +82,7 @@ __attribute__((always_inline)) INLINE static void gravity_first_init_gpart( ...@@ -82,7 +82,7 @@ __attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
* *
* @param gp The particle to act upon * @param gp The particle to act upon
*/ */
__attribute__((always_inline)) INLINE static void gravity_init_part( __attribute__((always_inline)) INLINE static void gravity_init_gpart(
struct gpart* gp) { struct gpart* gp) {
/* Zero the acceleration */ /* Zero the acceleration */
......
...@@ -136,15 +136,12 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -136,15 +136,12 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Some smoothing length multiples. */ /* Some smoothing length multiples. */
const float h = p->h; const float h = p->h;
const float ih = 1.0f / h; const float ih = 1.0f / h;
const float ih2 = ih * ih;
const float ih4 = ih2 * ih2;
/* 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);
const float normRot_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] + 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[1] * p->density.rot_v[1] +
p->density.rot_v[2] * p->density.rot_v[2]) / p->density.rot_v[2] * p->density.rot_v[2]);
p->rho * ih4;
/* Compute this particle's sound speed. */ /* Compute this particle's sound speed. */
const float u = p->u; const float u = p->u;
......
...@@ -42,7 +42,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list, ...@@ -42,7 +42,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY, list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
UNIT_CONV_LENGTH, parts, h); UNIT_CONV_LENGTH, parts, h);
list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY, list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
UNIT_CONV_ENERGY, parts, u); UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY, list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
UNIT_CONV_NO_UNITS, parts, id); UNIT_CONV_NO_UNITS, parts, id);
list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL, list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
...@@ -72,7 +72,7 @@ void hydro_write_particles(struct part* parts, struct io_props* list, ...@@ -72,7 +72,7 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass); io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH, list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
parts, h); parts, h);
list[4] = io_make_output_field("Entropy", FLOAT, 1, list[4] = io_make_output_field("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u); UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1, list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id); UNIT_CONV_NO_UNITS, parts, id);
......
...@@ -69,43 +69,43 @@ struct part { ...@@ -69,43 +69,43 @@ struct part {
float alpha; float alpha;
/* Store density/force specific stuff. */ /* Store density/force specific stuff. */
// union { union {
struct { struct {
/* Particle velocity divergence. */ /* Particle velocity divergence. */
float div_v; float div_v;
/* Derivative of particle number density. */ /* Derivative of particle number density. */
float wcount_dh; float wcount_dh;
/* Particle velocity curl. */ /* Particle velocity curl. */
float rot_v[3]; float rot_v[3];
/* Particle number density. */ /* Particle number density. */
float wcount; float wcount;
} density; } density;
struct { struct {
/* Balsara switch */ /* Balsara switch */
float balsara; float balsara;
/* Aggregate quantities. */ /* Aggregate quantities. */
float POrho2; float POrho2;
/* Change in particle energy over time. */ /* Change in particle energy over time. */
float u_dt; float u_dt;
/* Signal velocity */ /* Signal velocity */
float v_sig; float v_sig;
/* Sound speed */ /* Sound speed */
float c; float c;
} force; } force;
//}; };
/* Particle mass. */ /* Particle mass. */
float mass; float mass;
......
...@@ -65,7 +65,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part( ...@@ -65,7 +65,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->density.wcount_dh = 0.f; p->density.wcount_dh = 0.f;
p->rho = 0.f; p->rho = 0.f;
p->rho_dh = 0.f; p->rho_dh = 0.f;
p->div_v = 0.f; p->density.div_v = 0.f;
p->density.rot_v[0] = 0.f; p->density.rot_v[0] = 0.f;
p->density.rot_v[1] = 0.f; p->density.rot_v[1] = 0.f;
p->density.rot_v[2] = 0.f; p->density.rot_v[2] = 0.f;
...@@ -111,7 +111,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( ...@@ -111,7 +111,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.rot_v[2] *= ih4 * irho; p->density.rot_v[2] *= ih4 * irho;
/* Finish calculation of the velocity divergence */ /* Finish calculation of the velocity divergence */
p->div_v *= ih4 * irho; p->density.div_v *= ih4 * irho;
} }
/** /**
...@@ -128,17 +128,31 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -128,17 +128,31 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp, int ti_current, struct part *restrict p, struct xpart *restrict xp, int ti_current,
double timeBase) { double timeBase) {
const float fac_mu = 1.f; /* Will change with cosmological integration */
/* Compute the norm of the curl */ /* Compute the norm of the curl */
p->force.curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] + const float curl_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[1] * p->density.rot_v[1] +
p->density.rot_v[2] * p->density.rot_v[2]); p->density.rot_v[2] * p->density.rot_v[2]);
/* Compute the pressure */ /* Compute the pressure */
const float dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase; const float dt = (ti_current - (p->ti_begin + p->ti_end) * 0.5f) * timeBase;
p->force.pressure = (p->entropy + p->entropy_dt * dt) * pow_gamma(p->rho); const float pressure = (p->entropy + p->force.entropy_dt * dt) * pow_gamma(p->rho);
/* Divide the pressure by the density and density gradient */
const float P_over_rho = pressure / (p->rho * p->rho) * p->rho_dh;
/* Compute the sound speed */ /* Compute the sound speed */
p->force.soundspeed = sqrtf(hydro_gamma * p->force.pressure / p->rho); const float soundspeed = sqrtf(hydro_gamma * pressure / p->rho);
/* Compute the Balsara switch */
float balsara = fabsf(p->density.div_v) /
(fabsf(p->density.div_v) + curl_v + 0.0001f * p->force.soundspeed / fac_mu / p->h);
/* Update variables. */
p->force.P_over_rho = P_over_rho;
p->force.soundspeed = soundspeed;
p->force.balsara = balsara;
} }
/** /**
...@@ -160,7 +174,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( ...@@ -160,7 +174,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
p->h_dt = 0.0f; p->h_dt = 0.0f;
/* Reset the time derivatives. */ /* Reset the time derivatives. */
p->entropy_dt = 0.0f; p->force.entropy_dt = 0.0f;
/* Reset maximal signal velocity */ /* Reset maximal signal velocity */
p->force.v_sig = 0.0f; p->force.v_sig = 0.0f;
...@@ -181,11 +195,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( ...@@ -181,11 +195,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Drift the pressure */ /* Drift the pressure */
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase; const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
p->force.pressure = const float pressure =
(p->entropy + p->entropy_dt * dt_entr) * pow_gamma(p->rho); (p->entropy + p->force.entropy_dt * dt_entr) * pow_gamma(p->rho);
/* Divide the pressure by the density and density gradient */
const float P_over_rho = pressure / (p->rho * p->rho) * p->rho_dh;
/* Compute the new sound speed */ /* Compute the new sound speed */
p->force.soundspeed = sqrtf(hydro_gamma * p->force.pressure / p->rho); const float soundspeed = sqrtf(hydro_gamma * pressure / p->rho);
/* Update variables */
p->force.P_over_rho = P_over_rho;
p->force.soundspeed = soundspeed;
} }
/** /**
...@@ -198,7 +219,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( ...@@ -198,7 +219,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
__attribute__((always_inline)) INLINE static void hydro_end_force( __attribute__((always_inline)) INLINE static void hydro_end_force(
struct part *restrict p) { struct part *restrict p) {
p->entropy_dt *= hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho); p->force.entropy_dt *= hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
} }
/** /**
...@@ -214,15 +235,15 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -214,15 +235,15 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
float half_dt) { float half_dt) {
/* Do not decrease the entropy (temperature) by more than a factor of 2*/ /* Do not decrease the entropy (temperature) by more than a factor of 2*/
const float entropy_change = p->entropy_dt * dt; const float entropy_change = p->force.entropy_dt * dt;
if (entropy_change > -0.5f * p->entropy) if (entropy_change > -0.5f * p->entropy)
p->entropy += entropy_change; p->entropy += entropy_change;
else else
p->entropy *= 0.5f; p->entropy *= 0.5f;
/* Do not 'overcool' when timestep increases */ /* Do not 'overcool' when timestep increases */
if (p->entropy + 0.5f * p->entropy_dt * dt < 0.5f * p->entropy) if (p->entropy + 0.5f * p->force.entropy_dt * dt < 0.5f * p->entropy)
p->entropy_dt = -0.5f * p->entropy / dt; p->force.entropy_dt = -0.5f * p->entropy / dt;
} }
/** /**
...@@ -248,7 +269,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( ...@@ -248,7 +269,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy( __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const struct part *restrict p, float dt) { const struct part *restrict p, float dt) {
const float entropy = p->entropy + p->entropy_dt * dt; const float entropy = p->entropy + p->force.entropy_dt * dt;
return entropy * pow_gamma_minus_one(p->rho) * hydro_one_over_gamma_minus_one; return entropy * pow_gamma_minus_one(p->rho) * hydro_one_over_gamma_minus_one;
} }
...@@ -23,15 +23,15 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle( ...@@ -23,15 +23,15 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"x=[%.3e,%.3e,%.3e], " "x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n " "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"h=%.3e, " "h=%.3e, "
"wcount=%d, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, " "wcount=%d, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, P_over_rho=%.3e, "
"S=%.3e, " "S=%.3e, "
"dS/dt=%.3e, c=%.3e\n" "dS/dt=%.3e, c=%.3e\n"
"divV=%.3e, curlV=%.3e, rotV=[%.3e,%.3e,%.3e] \n " "divV=%.3e, rotV=[%.3e,%.3e,%.3e], balsara=%.3e \n "
"v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n", "v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0], p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2], xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
p->h, (int)p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh, p->h, (int)p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh,
p->rho, p->force.pressure, p->entropy, p->entropy_dt, p->force.soundspeed, p->rho, p->force.P_over_rho, p->entropy, p->force.entropy_dt, p->force.soundspeed,
p->div_v, p->force.curl_v, p->density.rot_v[0], p->density.rot_v[1], p->density.div_v, p->density.rot_v[0], p->density.rot_v[1],
p->density.rot_v[2], p->force.v_sig, p->h_dt, p->ti_begin, p->ti_end); p->density.rot_v[2], p->force.balsara, p->force.v_sig, p->h_dt, p->ti_begin, p->ti_end);
} }
...@@ -86,8 +86,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( ...@@ -86,8 +86,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
dv[2] = pi->v[2] - pj->v[2]; dv[2] = pi->v[2] - pj->v[2];
const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]; const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
pi->div_v -= faci * dvdr; pi->density.div_v -= faci * dvdr;
pj->div_v -= facj * dvdr; pj->density.div_v -= facj * dvdr;
/* Compute dv cross r */ /* Compute dv cross r */
curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1]; curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
...@@ -209,13 +209,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density( ...@@ -209,13 +209,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
pi[k]->rho_dh -= rhoi_dh.f[k]; pi[k]->rho_dh -= rhoi_dh.f[k];
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]->div_v -= div_vi.f[k]; pi[k]->density.div_v -= div_vi.f[k];
for (j = 0; j < 3; j++) pi[k]->density.rot_v[j] += curl_vi[j].f[k]; for (j = 0; j < 3; j++) pi[k]->density.rot_v[j] += curl_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]->div_v -= div_vj.f[k]; pj[k]->density.div_v -= div_vj.f[k];
for (j = 0; j < 3; j++) pj[k]->density.rot_v[j] += curl_vj[j].f[k]; for (j = 0; j < 3; j++) pj[k]->density.rot_v[j] += curl_vj[j].f[k];
} }
...@@ -263,7 +263,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( ...@@ -263,7 +263,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
dv[1] = pi->v[1] - pj->v[1]; dv[1] = pi->v[1] - pj->v[1];
dv[2] = pi->v[2] - pj->v[2]; dv[2] = pi->v[2] - pj->v[2];
const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]; const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
pi->div_v -= fac * dvdr; pi->density.div_v -= fac * dvdr;
/* Compute dv cross r */ /* Compute dv cross r */
curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1]; curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
...@@ -363,7 +363,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj, ...@@ -363,7 +363,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
pi[k]->rho_dh -= rhoi_dh.f[k]; pi[k]->rho_dh -= rhoi_dh.f[k];
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]->div_v -= div_vi.f[k]; pi[k]->density.div_v -= div_vi.f[k];
for (j = 0; j < 3; j++) pi[k]->density.rot_v[j] += curl_vi[j].f[k]; for (j = 0; j < 3; j++) pi[k]->density.rot_v[j] += curl_vi[j].f[k];
} }
...@@ -393,8 +393,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -393,8 +393,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float mj = pj->mass; const float mj = pj->mass;
const float rhoi = pi->rho; const float rhoi = pi->rho;
const float rhoj = pj->rho; const float rhoj = pj->rho;
const float pressurei = pi->force.pressure;
const float pressurej = pj->force.pressure;
/* Get the kernel for hi. */ /* Get the kernel for hi. */
const float hi_inv = 1.0f / hi; const float hi_inv = 1.0f / hi;
...@@ -411,8 +409,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -411,8 +409,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float wj_dr = hj2_inv * hj2_inv * wj_dx; const float wj_dr = hj2_inv * hj2_inv * wj_dx;
/* Compute gradient terms */ /* Compute gradient terms */
const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh; const float P_over_rho_i = pi->force.P_over_rho;
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh; const float P_over_rho_j = pj->force.P_over_rho;
/* Compute sound speeds */ /* Compute sound speeds */
const float ci = pi->force.soundspeed; const float ci = pi->force.soundspeed;
...@@ -424,13 +422,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -424,13 +422,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
(pi->v[2] - pj->v[2]) * dx[2]; (pi->v[2] - pj->v[2]) * dx[2];
/* Balsara term */ /* Balsara term */
const float balsara_i = const float balsara_i = pi->force.balsara;
fabsf(pi->div_v) / const float balsara_j = pj->force.balsara;
(fabsf(pi->div_v) + pi->force.curl_v + 0.0001f * ci / fac_mu / hi);
const float balsara_j =
fabsf(pj->div_v) /
(fabsf(pj->div_v) + pj->force.curl_v + 0.0001f * cj / fac_mu / hj);
/* Are the particles moving towards each others ? */ /* Are the particles moving towards each others ? */
const float omega_ij = fminf(dvdr, 0.f); const float omega_ij = fminf(dvdr, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
...@@ -468,8 +462,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -468,8 +462,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->force.v_sig = fmaxf(pj->force.v_sig, v_sig); pj->force.v_sig = fmaxf(pj->force.v_sig, v_sig);
/* Change in entropy */ /* Change in entropy */
pi->entropy_dt += 0.5f * mj * visc_term * dvdr; pi->force.entropy_dt += 0.5f * mj * visc_term * dvdr;
pj->entropy_dt -= 0.5f * mi * visc_term * dvdr; pj->force.entropy_dt -= 0.5f * mi * visc_term * dvdr;
} }
/** /**
...@@ -501,8 +495,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -501,8 +495,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float mj = pj->mass; const float mj = pj->mass;
const float rhoi = pi->rho; const float rhoi = pi->rho;
const float rhoj = pj->rho; const float rhoj = pj->rho;
const float pressurei = pi->force.pressure;
const float pressurej = pj->force.pressure;
/* Get the kernel for hi. */ /* Get the kernel for hi. */
const float hi_inv = 1.0f / hi; const float hi_inv = 1.0f / hi;
...@@ -519,8 +511,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -519,8 +511,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float wj_dr = hj2_inv * hj2_inv * wj_dx; const float wj_dr = hj2_inv * hj2_inv * wj_dx;
/* Compute gradient terms */ /* Compute gradient terms */
const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh; const float P_over_rho_i = pi->force.P_over_rho;
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh; const float P_over_rho_j = pj->force.P_over_rho;
/* Compute sound speeds */ /* Compute sound speeds */
const float ci = pi->force.soundspeed; const float ci = pi->force.soundspeed;
...@@ -532,12 +524,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -532,12 +524,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
(pi->v[2] - pj->v[2]) * dx[2]; (pi->v[2] - pj->v[2]) * dx[2];
/* Balsara term */ /* Balsara term */
const float balsara_i = const float balsara_i = pi->force.balsara;
fabsf(pi->div_v) / const float balsara_j = pj->force.balsara;
(fabsf(pi->div_v) + pi->force.curl_v + 0.0001f * ci / fac_mu / hi);
const float balsara_j =
fabsf(pj->div_v) /
(fabsf(pj->div_v) + pj->force.curl_v + 0.0001f * cj / fac_mu / hj);
/* Are the particles moving towards each others ? */ /* Are the particles moving towards each others ? */
const float omega_ij = fminf(dvdr, 0.f); const float omega_ij = fminf(dvdr, 0.f);
...@@ -570,7 +558,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -570,7 +558,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig); pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
/* Change in entropy */ /* Change in entropy */
pi->entropy_dt += 0.5f * mj * visc_term * dvdr; pi->force.entropy_dt += 0.5f * mj * visc_term * dvdr;
} }
/** /**
......