Commit cdffd1b9 authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Moved GizmoMFM particle variables around and put some of them in unions.

parent 422560d9
......@@ -57,7 +57,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
vrel[2] = p->v[2] - xp->v_full[2];
float vmax =
sqrtf(vrel[0] * vrel[0] + vrel[1] * vrel[1] + vrel[2] * vrel[2]) +
sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
sqrtf(hydro_gamma * p->P / p->rho);
vmax = max(vmax, p->timestepvars.vmax);
// MATTHIEU: Bert is this correct? Do we need more cosmology terms here?
......@@ -152,7 +152,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
time the density loop is repeated, and the whole point of storing wcorr
is to have a way of remembering that we need more neighbours for this
particle */
p->density.wcorr = 1.0f;
p->geometry.wcorr = 1.0f;
}
/**
......@@ -268,7 +268,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
hydro_dimension_inv * sqrtf(condition_number_E * condition_number_Einv);
if (condition_number > const_gizmo_max_condition_number &&
p->density.wcorr > const_gizmo_min_wcorr) {
p->geometry.wcorr > const_gizmo_min_wcorr) {
#ifdef GIZMO_PATHOLOGICAL_ERROR
error("Condition number larger than %g (%g)!",
const_gizmo_max_condition_number, condition_number);
......@@ -278,7 +278,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
condition_number, const_gizmo_max_condition_number, p->id);
#endif
/* add a correction to the number of neighbours for this particle */
p->density.wcorr *= const_gizmo_w_correction_factor;
p->geometry.wcorr *= const_gizmo_w_correction_factor;
}
/* compute primitive variables */
......@@ -300,7 +300,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
momentum[0] = p->conserved.momentum[0];
momentum[1] = p->conserved.momentum[1];
momentum[2] = p->conserved.momentum[2];
p->primitives.rho = m / volume;
p->rho = m / volume;
if (m == 0.0f) {
p->v[0] = 0.0f;
p->v[1] = 0.0f;
......@@ -315,7 +315,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
#ifdef EOS_ISOTHERMAL_GAS
/* although the pressure is not formally used anywhere if an isothermal eos
has been selected, we still make sure it is set to the correct value */
p->primitives.P = gas_pressure_from_internal_energy(p->primitives.rho, 0.0f);
p->P = gas_pressure_from_internal_energy(p->rho, 0.0f);
#else
float energy = p->conserved.energy;
......@@ -328,17 +328,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* energy contains the total thermal energy, we want the specific energy.
this is why we divide by the volume, and not by the density */
p->primitives.P = hydro_gamma_minus_one * energy / volume;
p->P = hydro_gamma_minus_one * energy / volume;
#endif
/* sanity checks */
gizmo_check_physical_quantities("density", "pressure", p->primitives.rho,
p->v[0], p->v[1], p->v[2], p->primitives.P);
gizmo_check_physical_quantities("density", "pressure", p->rho, p->v[0],
p->v[1], p->v[2], p->P);
/* Add a correction factor to wcount (to force a neighbour number increase if
the geometry matrix is close to singular) */
p->density.wcount *= p->density.wcorr;
p->density.wcount_dh *= p->density.wcorr;
p->density.wcount *= p->geometry.wcorr;
p->density.wcount_dh *= p->geometry.wcorr;
}
/**
......@@ -453,10 +453,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const struct cosmology* cosmo) {
/* Initialise values that are used in the force loop */
p->conserved.flux.momentum[0] = 0.0f;
p->conserved.flux.momentum[1] = 0.0f;
p->conserved.flux.momentum[2] = 0.0f;
p->conserved.flux.energy = 0.0f;
p->flux.momentum[0] = 0.0f;
p->flux.momentum[1] = 0.0f;
p->flux.momentum[2] = 0.0f;
p->flux.energy = 0.0f;
}
/**
......@@ -533,22 +533,20 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
if (p->conserved.mass > 0.0f) {
const float m_inv = 1.0f / p->conserved.mass;
p->v[0] += p->conserved.flux.momentum[0] * dt_drift * m_inv;
p->v[1] += p->conserved.flux.momentum[1] * dt_drift * m_inv;
p->v[2] += p->conserved.flux.momentum[2] * dt_drift * m_inv;
p->v[0] += p->flux.momentum[0] * dt_drift * m_inv;
p->v[1] += p->flux.momentum[1] * dt_drift * m_inv;
p->v[2] += p->flux.momentum[2] * dt_drift * m_inv;
#if !defined(EOS_ISOTHERMAL_GAS)
#ifdef GIZMO_TOTAL_ENERGY
const float Etot =
p->conserved.energy + p->conserved.flux.energy * dt_therm;
const float Etot = p->conserved.energy + p->flux.energy * dt_therm;
const float v2 =
(p->v[0] * p->v[0] + p->v[1] * p->v[1] + p->v[2] * p->v[2]);
const float u = (Etot * m_inv - 0.5f * v2);
#else
const float u =
(p->conserved.energy + p->conserved.flux.energy * dt_therm) * m_inv;
const float u = (p->conserved.energy + p->flux.energy * dt_therm) * m_inv;
#endif
p->primitives.P = hydro_gamma_minus_one * u * p->primitives.rho;
p->P = hydro_gamma_minus_one * u * p->rho;
#endif
}
......@@ -564,8 +562,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
}
#endif
gizmo_check_physical_quantities("density", "pressure", p->primitives.rho,
p->v[0], p->v[1], p->v[2], p->primitives.P);
gizmo_check_physical_quantities("density", "pressure", p->rho, p->v[0],
p->v[1], p->v[2], p->P);
}
/**
......@@ -609,15 +607,15 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
float a_grav[3];
/* Update conserved variables (note: the mass does not change). */
p->conserved.momentum[0] += p->conserved.flux.momentum[0] * dt;
p->conserved.momentum[1] += p->conserved.flux.momentum[1] * dt;
p->conserved.momentum[2] += p->conserved.flux.momentum[2] * dt;
p->conserved.momentum[0] += p->flux.momentum[0] * dt;
p->conserved.momentum[1] += p->flux.momentum[1] * dt;
p->conserved.momentum[2] += p->flux.momentum[2] * dt;
#if defined(EOS_ISOTHERMAL_GAS)
/* We use the EoS equation in a sneaky way here just to get the constant u */
p->conserved.energy =
p->conserved.mass * gas_internal_energy_from_entropy(0.0f, 0.0f);
#else
p->conserved.energy += p->conserved.flux.energy * dt;
p->conserved.energy += p->flux.energy * dt;
#endif
/* Apply the minimal energy limit */
......@@ -625,7 +623,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
hydro_props->minimal_internal_energy * cosmo->a_factor_internal_energy;
if (p->conserved.energy < min_energy * p->conserved.mass) {
p->conserved.energy = min_energy * p->conserved.mass;
p->conserved.flux.energy = 0.0f;
p->flux.energy = 0.0f;
}
gizmo_check_physical_quantities(
......@@ -639,7 +637,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
error(
"Negative energy after conserved variables update (energy: %g, "
"denergy: %g)!",
p->conserved.energy, p->conserved.flux.energy);
p->conserved.energy, p->flux.energy);
}
#endif
......@@ -662,7 +660,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* Set the velocities: */
/* We first set the particle velocity */
if (p->conserved.mass > 0.0f && p->primitives.rho > 0.0f) {
if (p->conserved.mass > 0.0f && p->rho > 0.0f) {
const float inverse_mass = 1.0f / p->conserved.mass;
......@@ -685,7 +683,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
}
/* reset wcorr */
p->density.wcorr = 1.0f;
p->geometry.wcorr = 1.0f;
}
/**
......@@ -696,9 +694,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_internal_energy(const struct part* restrict p) {
if (p->primitives.rho > 0.0f) {
return gas_internal_energy_from_pressure(p->primitives.rho,
p->primitives.P);
if (p->rho > 0.0f) {
return gas_internal_energy_from_pressure(p->rho, p->P);
} else {
return 0.0f;
}
......@@ -726,8 +723,8 @@ hydro_get_physical_internal_energy(const struct part* restrict p,
__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
const struct part* restrict p) {
if (p->primitives.rho > 0.0f) {
return gas_entropy_from_pressure(p->primitives.rho, p->primitives.P);
if (p->rho > 0.0f) {
return gas_entropy_from_pressure(p->rho, p->P);
} else {
return 0.0f;
}
......@@ -755,8 +752,8 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_soundspeed(const struct part* restrict p) {
if (p->primitives.rho > 0.0f) {
return gas_soundspeed_from_pressure(p->primitives.rho, p->primitives.P);
if (p->rho > 0.0f) {
return gas_soundspeed_from_pressure(p->rho, p->P);
} else {
return 0.0f;
}
......@@ -783,7 +780,7 @@ hydro_get_physical_soundspeed(const struct part* restrict p,
__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
const struct part* restrict p) {
return p->primitives.P;
return p->P;
}
/**
......@@ -795,7 +792,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
const struct part* restrict p, const struct cosmology* cosmo) {
return cosmo->a_factor_pressure * p->primitives.P;
return cosmo->a_factor_pressure * p->P;
}
/**
......@@ -836,12 +833,9 @@ __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
if (p->conserved.mass > 0.0f) {
const float inverse_mass = 1.0f / p->conserved.mass;
v[0] =
p->v[0] + p->conserved.flux.momentum[0] * dt_kick_hydro * inverse_mass;
v[1] =
p->v[1] + p->conserved.flux.momentum[1] * dt_kick_hydro * inverse_mass;
v[2] =
p->v[2] + p->conserved.flux.momentum[2] * dt_kick_hydro * inverse_mass;
v[0] = p->v[0] + p->flux.momentum[0] * dt_kick_hydro * inverse_mass;
v[1] = p->v[1] + p->flux.momentum[1] * dt_kick_hydro * inverse_mass;
v[2] = p->v[2] + p->flux.momentum[2] * dt_kick_hydro * inverse_mass;
} else {
v[0] = p->v[0];
v[1] = p->v[1];
......@@ -862,7 +856,7 @@ __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
const struct part* restrict p) {
return p->primitives.rho;
return p->rho;
}
/**
......@@ -874,7 +868,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
__attribute__((always_inline)) INLINE static float hydro_get_physical_density(
const struct part* restrict p, const struct cosmology* cosmo) {
return cosmo->a3_inv * p->primitives.rho;
return cosmo->a3_inv * p->rho;
}
/**
......@@ -900,7 +894,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
(p->conserved.momentum[0] * p->v[0] + p->conserved.momentum[1] * p->v[1] +
p->conserved.momentum[2] * p->v[2]);
#endif
p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u;
p->P = hydro_gamma_minus_one * p->rho * u;
}
/**
......@@ -915,7 +909,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
__attribute__((always_inline)) INLINE static void hydro_set_entropy(
struct part* restrict p, float S) {
p->conserved.energy = S * pow_gamma_minus_one(p->primitives.rho) *
p->conserved.energy = S * pow_gamma_minus_one(p->rho) *
hydro_one_over_gamma_minus_one * p->conserved.mass;
#ifdef GIZMO_TOTAL_ENERGY
/* add the kinetic energy */
......@@ -924,7 +918,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
(p->conserved.momentum[0] * p->v[0] + p->conserved.momentum[1] * p->v[1] +
p->conserved.momentum[2] * p->v[2]);
#endif
p->primitives.P = S * pow_gamma(p->primitives.rho);
p->P = S * pow_gamma(p->rho);
}
/**
......@@ -949,7 +943,7 @@ hydro_set_init_internal_energy(struct part* p, float u_init) {
(p->conserved.momentum[0] * p->v[0] + p->conserved.momentum[1] * p->v[1] +
p->conserved.momentum[2] * p->v[2]);
#endif
p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u_init;
p->P = hydro_gamma_minus_one * p->rho * u_init;
}
#endif /* SWIFT_GIZMO_MFM_HYDRO_H */
......@@ -27,7 +27,6 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"a=[%.3e,%.3e,%.3e], "
"h=%.3e, "
"time_bin=%d, "
"primitives={"
"rho=%.3e, "
"P=%.3e, "
"gradients={"
......@@ -38,7 +37,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"rho=[%.3e,%.3e], "
"v=[[%.3e,%.3e],[%.3e,%.3e],[%.3e,%.3e]], "
"P=[%.3e,%.3e], "
"maxr=%.3e}}, "
"maxr=%.3e}, "
"conserved={"
"momentum=[%.3e,%.3e,%.3e], "
"mass=%.3e, "
......@@ -52,23 +51,18 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"wcount_dh=%.3e, "
"wcount=%.3e}\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
p->a_hydro[1], p->a_hydro[2], p->h, p->time_bin, p->primitives.rho,
p->primitives.P, p->primitives.gradients.rho[0],
p->primitives.gradients.rho[1], p->primitives.gradients.rho[2],
p->primitives.gradients.v[0][0], p->primitives.gradients.v[0][1],
p->primitives.gradients.v[0][2], p->primitives.gradients.v[1][0],
p->primitives.gradients.v[1][1], p->primitives.gradients.v[1][2],
p->primitives.gradients.v[2][0], p->primitives.gradients.v[2][1],
p->primitives.gradients.v[2][2], p->primitives.gradients.P[0],
p->primitives.gradients.P[1], p->primitives.gradients.P[2],
p->primitives.limiter.rho[0], p->primitives.limiter.rho[1],
p->primitives.limiter.v[0][0], p->primitives.limiter.v[0][1],
p->primitives.limiter.v[1][0], p->primitives.limiter.v[1][1],
p->primitives.limiter.v[2][0], p->primitives.limiter.v[2][1],
p->primitives.limiter.P[0], p->primitives.limiter.P[1],
p->primitives.limiter.maxr, p->conserved.momentum[0],
p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.mass,
p->conserved.energy, p->geometry.volume, p->geometry.matrix_E[0][0],
p->a_hydro[1], p->a_hydro[2], p->h, p->time_bin, p->rho, p->P,
p->gradients.rho[0], p->gradients.rho[1], p->gradients.rho[2],
p->gradients.v[0][0], p->gradients.v[0][1], p->gradients.v[0][2],
p->gradients.v[1][0], p->gradients.v[1][1], p->gradients.v[1][2],
p->gradients.v[2][0], p->gradients.v[2][1], p->gradients.v[2][2],
p->gradients.P[0], p->gradients.P[1], p->gradients.P[2],
p->limiter.rho[0], p->limiter.rho[1], p->limiter.v[0][0],
p->limiter.v[0][1], p->limiter.v[1][0], p->limiter.v[1][1],
p->limiter.v[2][0], p->limiter.v[2][1], p->limiter.P[0], p->limiter.P[1],
p->limiter.maxr, p->conserved.momentum[0], p->conserved.momentum[1],
p->conserved.momentum[2], p->conserved.mass, p->conserved.energy,
p->geometry.volume, p->geometry.matrix_E[0][0],
p->geometry.matrix_E[0][1], p->geometry.matrix_E[0][2],
p->geometry.matrix_E[1][0], p->geometry.matrix_E[1][1],
p->geometry.matrix_E[1][2], p->geometry.matrix_E[2][0],
......
......@@ -102,38 +102,28 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
const float xij_j[3] = {xfac * dx[0], xfac * dx[1], xfac * dx[2]};
float dWi[5];
dWi[0] = pi->primitives.gradients.rho[0] * xij_i[0] +
pi->primitives.gradients.rho[1] * xij_i[1] +
pi->primitives.gradients.rho[2] * xij_i[2];
dWi[1] = pi->primitives.gradients.v[0][0] * xij_i[0] +
pi->primitives.gradients.v[0][1] * xij_i[1] +
pi->primitives.gradients.v[0][2] * xij_i[2];
dWi[2] = pi->primitives.gradients.v[1][0] * xij_i[0] +
pi->primitives.gradients.v[1][1] * xij_i[1] +
pi->primitives.gradients.v[1][2] * xij_i[2];
dWi[3] = pi->primitives.gradients.v[2][0] * xij_i[0] +
pi->primitives.gradients.v[2][1] * xij_i[1] +
pi->primitives.gradients.v[2][2] * xij_i[2];
dWi[4] = pi->primitives.gradients.P[0] * xij_i[0] +
pi->primitives.gradients.P[1] * xij_i[1] +
pi->primitives.gradients.P[2] * xij_i[2];
dWi[0] = pi->gradients.rho[0] * xij_i[0] + pi->gradients.rho[1] * xij_i[1] +
pi->gradients.rho[2] * xij_i[2];
dWi[1] = pi->gradients.v[0][0] * xij_i[0] + pi->gradients.v[0][1] * xij_i[1] +
pi->gradients.v[0][2] * xij_i[2];
dWi[2] = pi->gradients.v[1][0] * xij_i[0] + pi->gradients.v[1][1] * xij_i[1] +
pi->gradients.v[1][2] * xij_i[2];
dWi[3] = pi->gradients.v[2][0] * xij_i[0] + pi->gradients.v[2][1] * xij_i[1] +
pi->gradients.v[2][2] * xij_i[2];
dWi[4] = pi->gradients.P[0] * xij_i[0] + pi->gradients.P[1] * xij_i[1] +
pi->gradients.P[2] * xij_i[2];
float dWj[5];
dWj[0] = pj->primitives.gradients.rho[0] * xij_j[0] +
pj->primitives.gradients.rho[1] * xij_j[1] +
pj->primitives.gradients.rho[2] * xij_j[2];
dWj[1] = pj->primitives.gradients.v[0][0] * xij_j[0] +
pj->primitives.gradients.v[0][1] * xij_j[1] +
pj->primitives.gradients.v[0][2] * xij_j[2];
dWj[2] = pj->primitives.gradients.v[1][0] * xij_j[0] +
pj->primitives.gradients.v[1][1] * xij_j[1] +
pj->primitives.gradients.v[1][2] * xij_j[2];
dWj[3] = pj->primitives.gradients.v[2][0] * xij_j[0] +
pj->primitives.gradients.v[2][1] * xij_j[1] +
pj->primitives.gradients.v[2][2] * xij_j[2];
dWj[4] = pj->primitives.gradients.P[0] * xij_j[0] +
pj->primitives.gradients.P[1] * xij_j[1] +
pj->primitives.gradients.P[2] * xij_j[2];
dWj[0] = pj->gradients.rho[0] * xij_j[0] + pj->gradients.rho[1] * xij_j[1] +
pj->gradients.rho[2] * xij_j[2];
dWj[1] = pj->gradients.v[0][0] * xij_j[0] + pj->gradients.v[0][1] * xij_j[1] +
pj->gradients.v[0][2] * xij_j[2];
dWj[2] = pj->gradients.v[1][0] * xij_j[0] + pj->gradients.v[1][1] * xij_j[1] +
pj->gradients.v[1][2] * xij_j[2];
dWj[3] = pj->gradients.v[2][0] * xij_j[0] + pj->gradients.v[2][1] * xij_j[1] +
pj->gradients.v[2][2] * xij_j[2];
dWj[4] = pj->gradients.P[0] * xij_j[0] + pj->gradients.P[1] * xij_j[1] +
pj->gradients.P[2] * xij_j[2];
/* Apply the slope limiter at this interface */
hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
......
This diff is collapsed.
......@@ -28,24 +28,24 @@
__attribute__((always_inline)) INLINE static void hydro_gradients_init(
struct part *p) {
p->primitives.gradients.rho[0] = 0.0f;
p->primitives.gradients.rho[1] = 0.0f;
p->primitives.gradients.rho[2] = 0.0f;
p->gradients.rho[0] = 0.0f;
p->gradients.rho[1] = 0.0f;
p->gradients.rho[2] = 0.0f;
p->primitives.gradients.v[0][0] = 0.0f;
p->primitives.gradients.v[0][1] = 0.0f;
p->primitives.gradients.v[0][2] = 0.0f;
p->gradients.v[0][0] = 0.0f;
p->gradients.v[0][1] = 0.0f;
p->gradients.v[0][2] = 0.0f;
p->primitives.gradients.v[1][0] = 0.0f;
p->primitives.gradients.v[1][1] = 0.0f;
p->primitives.gradients.v[1][2] = 0.0f;
p->primitives.gradients.v[2][0] = 0.0f;
p->primitives.gradients.v[2][1] = 0.0f;
p->primitives.gradients.v[2][2] = 0.0f;
p->gradients.v[1][0] = 0.0f;
p->gradients.v[1][1] = 0.0f;
p->gradients.v[1][2] = 0.0f;
p->gradients.v[2][0] = 0.0f;
p->gradients.v[2][1] = 0.0f;
p->gradients.v[2][2] = 0.0f;
p->primitives.gradients.P[0] = 0.0f;
p->primitives.gradients.P[1] = 0.0f;
p->primitives.gradients.P[2] = 0.0f;
p->gradients.P[0] = 0.0f;
p->gradients.P[1] = 0.0f;
p->gradients.P[2] = 0.0f;
hydro_slope_limit_cell_init(p);
}
......@@ -73,40 +73,25 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
kernel_deval(xi, &wi, &wi_dx);
/* very basic gradient estimate */
pi->primitives.gradients.rho[0] -=
wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.rho[1] -=
wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.rho[2] -=
wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pi->primitives.gradients.v[0][0] -=
wi_dx * dx[0] * (pi->v[0] - pj->v[0]) * r_inv;
pi->primitives.gradients.v[0][1] -=
wi_dx * dx[1] * (pi->v[0] - pj->v[0]) * r_inv;
pi->primitives.gradients.v[0][2] -=
wi_dx * dx[2] * (pi->v[0] - pj->v[0]) * r_inv;
pi->primitives.gradients.v[1][0] -=
wi_dx * dx[0] * (pi->v[1] - pj->v[1]) * r_inv;
pi->primitives.gradients.v[1][1] -=
wi_dx * dx[1] * (pi->v[1] - pj->v[1]) * r_inv;
pi->primitives.gradients.v[1][2] -=
wi_dx * dx[2] * (pi->v[1] - pj->v[1]) * r_inv;
pi->primitives.gradients.v[2][0] -=
wi_dx * dx[0] * (pi->v[2] - pj->v[2]) * r_inv;
pi->primitives.gradients.v[2][1] -=
wi_dx * dx[1] * (pi->v[2] - pj->v[2]) * r_inv;
pi->primitives.gradients.v[2][2] -=
wi_dx * dx[2] * (pi->v[2] - pj->v[2]) * r_inv;
pi->primitives.gradients.P[0] -=
wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
pi->primitives.gradients.P[1] -=
wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
pi->primitives.gradients.P[2] -=
wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
pi->gradients.rho[0] -= wi_dx * dx[0] * (pi->rho - pj->rho) * r_inv;
pi->gradients.rho[1] -= wi_dx * dx[1] * (pi->rho - pj->rho) * r_inv;
pi->gradients.rho[2] -= wi_dx * dx[2] * (pi->rho - pj->rho) * r_inv;
pi->gradients.v[0][0] -= wi_dx * dx[0] * (pi->v[0] - pj->v[0]) * r_inv;
pi->gradients.v[0][1] -= wi_dx * dx[1] * (pi->v[0] - pj->v[0]) * r_inv;
pi->gradients.v[0][2] -= wi_dx * dx[2] * (pi->v[0] - pj->v[0]) * r_inv;
pi->gradients.v[1][0] -= wi_dx * dx[0] * (pi->v[1] - pj->v[1]) * r_inv;
pi->gradients.v[1][1] -= wi_dx * dx[1] * (pi->v[1] - pj->v[1]) * r_inv;
pi->gradients.v[1][2] -= wi_dx * dx[2] * (pi->v[1] - pj->v[1]) * r_inv;
pi->gradients.v[2][0] -= wi_dx * dx[0] * (pi->v[2] - pj->v[2]) * r_inv;
pi->gradients.v[2][1] -= wi_dx * dx[1] * (pi->v[2] - pj->v[2]) * r_inv;
pi->gradients.v[2][2] -= wi_dx * dx[2] * (pi->v[2] - pj->v[2]) * r_inv;
pi->gradients.P[0] -= wi_dx * dx[0] * (pi->P - pj->P) * r_inv;
pi->gradients.P[1] -= wi_dx * dx[1] * (pi->P - pj->P) * r_inv;
pi->gradients.P[2] -= wi_dx * dx[2] * (pi->P - pj->P) * r_inv;
hydro_slope_limit_cell_collect(pi, pj, r);
......@@ -116,39 +101,24 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
kernel_deval(xj, &wj, &wj_dx);
/* signs are the same as before, since we swap i and j twice */
pj->primitives.gradients.rho[0] -=
wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pj->primitives.gradients.rho[1] -=
wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pj->primitives.gradients.rho[2] -=
wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
pj->primitives.gradients.v[0][0] -=
wj_dx * dx[0] * (pi->v[0] - pj->v[0]) * r_inv;
pj->primitives.gradients.v[0][1] -=
wj_dx * dx[1] * (pi->v[0] - pj->v[0]) * r_inv;
pj->primitives.gradients.v[0][2] -=
wj_dx * dx[2] * (pi->v[0] - pj->v[0]) * r_inv;
pj->primitives.gradients.v[1][0] -=
wj_dx * dx[0] * (pi->v[1] - pj->v[1]) * r_inv;
pj->primitives.gradients.v[1][1] -=
wj_dx * dx[1] * (pi->v[1] - pj->v[1]) * r_inv;
pj->primitives.gradients.v[1][2] -=
wj_dx * dx[2] * (pi->v[1] - pj->v[1]) * r_inv;
pj->primitives.gradients.v[2][0] -=
wj_dx * dx[0] * (pi->v[2] - pj->v[2]) * r_inv;
pj->primitives.gradients.v[2][1] -=
wj_dx * dx[1] * (pi->v[2] - pj->v[2]) * r_inv;
pj->primitives.gradients.v[2][2] -=
wj_dx * dx[2] * (pi->v[2] - pj->v[2]) * r_inv;
pj->primitives.gradients.P[0] -=
wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
pj->primitives.gradients.P[1] -=
wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
pj->primitives.gradients.P[2] -=
wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
pj->gradients.rho[0] -= wj_dx * dx[0] * (pi->rho - pj->rho) * r_inv;
pj->gradients.rho[1] -= wj_dx * dx[1] * (pi->rho - pj->rho) * r_inv;
pj->gradients.rho[2] -= wj_dx * dx[2] * (pi->rho - pj->rho) * r_inv;
pj->gradients.v[0][0] -= wj_dx * dx[0] * (pi->v[0] - pj->v[0]) * r_inv;