diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h
index 3ac40845236f2c22c38997893d3523add5bb4d45..28545e7c434800504fbad5629f0275189354d38f 100644
--- a/src/hydro/GizmoMFM/hydro.h
+++ b/src/hydro/GizmoMFM/hydro.h
@@ -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 */
diff --git a/src/hydro/GizmoMFM/hydro_debug.h b/src/hydro/GizmoMFM/hydro_debug.h
index 6dfe739f05f2423a6e08410d98d4b74652fe9578..e8b0914bd3cf6a99210399c6fc654e526319009f 100644
--- a/src/hydro/GizmoMFM/hydro_debug.h
+++ b/src/hydro/GizmoMFM/hydro_debug.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],
diff --git a/src/hydro/GizmoMFM/hydro_gradients.h b/src/hydro/GizmoMFM/hydro_gradients.h
index 964a2adcfe09b95c2a221af540e5e3ff0830dd67..063553b888a6d27fc7c24ec6a7213b14b0b359d8 100644
--- a/src/hydro/GizmoMFM/hydro_gradients.h
+++ b/src/hydro/GizmoMFM/hydro_gradients.h
@@ -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);
diff --git a/src/hydro/GizmoMFM/hydro_gradients_gizmo.h b/src/hydro/GizmoMFM/hydro_gradients_gizmo.h
index ca82b88b4972dee8c56c1e10b65ffada23d82593..089c4f5a510b8e5859f7128a96938de78a2f7075 100644
--- a/src/hydro/GizmoMFM/hydro_gradients_gizmo.h
+++ b/src/hydro/GizmoMFM/hydro_gradients_gizmo.h
@@ -28,25 +28,25 @@
__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->gradients.v[1][0] = 0.0f;
+ p->gradients.v[1][1] = 0.0f;
+ p->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[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);
}
@@ -80,111 +80,96 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
Bj[k][l] = pj->geometry.matrix_E[k][l];
}
}
- Wi[0] = pi->primitives.rho;
+ Wi[0] = pi->rho;
Wi[1] = pi->v[0];
Wi[2] = pi->v[1];
Wi[3] = pi->v[2];
- Wi[4] = pi->primitives.P;
- Wj[0] = pj->primitives.rho;
+ Wi[4] = pi->P;
+ Wj[0] = pj->rho;
Wj[1] = pj->v[0];
Wj[2] = pj->v[1];
Wj[3] = pj->v[2];
- Wj[4] = pj->primitives.P;
+ Wj[4] = pj->P;
/* Compute kernel of pi. */
const float hi_inv = 1.f / hi;
const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
- if (pi->density.wcorr > const_gizmo_min_wcorr) {
+ if (pi->geometry.wcorr > const_gizmo_min_wcorr) {
/* Compute gradients for pi */
/* there is a sign difference w.r.t. eqn. (6) because of the inverse
* definition of dx */
- pi->primitives.gradients.rho[0] +=
+ pi->gradients.rho[0] +=
(Wi[0] - Wj[0]) * wi *
(Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.rho[1] +=
+ pi->gradients.rho[1] +=
(Wi[0] - Wj[0]) * wi *
(Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.rho[2] +=
+ pi->gradients.rho[2] +=
(Wi[0] - Wj[0]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->primitives.gradients.v[0][0] +=
+ pi->gradients.v[0][0] +=
(Wi[1] - Wj[1]) * wi *
(Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.v[0][1] +=
+ pi->gradients.v[0][1] +=
(Wi[1] - Wj[1]) * wi *
(Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.v[0][2] +=
+ pi->gradients.v[0][2] +=
(Wi[1] - Wj[1]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->primitives.gradients.v[1][0] +=
+ pi->gradients.v[1][0] +=
(Wi[2] - Wj[2]) * wi *
(Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.v[1][1] +=
+ pi->gradients.v[1][1] +=
(Wi[2] - Wj[2]) * wi *
(Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.v[1][2] +=
+ pi->gradients.v[1][2] +=
(Wi[2] - Wj[2]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->primitives.gradients.v[2][0] +=
+ pi->gradients.v[2][0] +=
(Wi[3] - Wj[3]) * wi *
(Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.v[2][1] +=
+ pi->gradients.v[2][1] +=
(Wi[3] - Wj[3]) * wi *
(Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.v[2][2] +=
+ pi->gradients.v[2][2] +=
(Wi[3] - Wj[3]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->primitives.gradients.P[0] +=
+ pi->gradients.P[0] +=
(Wi[4] - Wj[4]) * wi *
(Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.P[1] +=
+ pi->gradients.P[1] +=
(Wi[4] - Wj[4]) * wi *
(Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.P[2] +=
+ pi->gradients.P[2] +=
(Wi[4] - Wj[4]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
} else {
/* The gradient matrix was not well-behaved, switch to SPH gradients */
- 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);
@@ -194,95 +179,80 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
- if (pj->density.wcorr > const_gizmo_min_wcorr) {
+ if (pj->geometry.wcorr > const_gizmo_min_wcorr) {
/* Compute gradients for pj */
/* there is no sign difference w.r.t. eqn. (6) because dx is now what we
* want
* it to be */
- pj->primitives.gradients.rho[0] +=
+ pj->gradients.rho[0] +=
(Wi[0] - Wj[0]) * wj *
(Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.rho[1] +=
+ pj->gradients.rho[1] +=
(Wi[0] - Wj[0]) * wj *
(Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.rho[2] +=
+ pj->gradients.rho[2] +=
(Wi[0] - Wj[0]) * wj *
(Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
- pj->primitives.gradients.v[0][0] +=
+ pj->gradients.v[0][0] +=
(Wi[1] - Wj[1]) * wj *
(Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.v[0][1] +=
+ pj->gradients.v[0][1] +=
(Wi[1] - Wj[1]) * wj *
(Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.v[0][2] +=
+ pj->gradients.v[0][2] +=
(Wi[1] - Wj[1]) * wj *
(Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
- pj->primitives.gradients.v[1][0] +=
+ pj->gradients.v[1][0] +=
(Wi[2] - Wj[2]) * wj *
(Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.v[1][1] +=
+ pj->gradients.v[1][1] +=
(Wi[2] - Wj[2]) * wj *
(Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.v[1][2] +=
+ pj->gradients.v[1][2] +=
(Wi[2] - Wj[2]) * wj *
(Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
- pj->primitives.gradients.v[2][0] +=
+ pj->gradients.v[2][0] +=
(Wi[3] - Wj[3]) * wj *
(Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.v[2][1] +=
+ pj->gradients.v[2][1] +=
(Wi[3] - Wj[3]) * wj *
(Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.v[2][2] +=
+ pj->gradients.v[2][2] +=
(Wi[3] - Wj[3]) * wj *
(Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
- pj->primitives.gradients.P[0] +=
+ pj->gradients.P[0] +=
(Wi[4] - Wj[4]) * wj *
(Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.P[1] +=
+ pj->gradients.P[1] +=
(Wi[4] - Wj[4]) * wj *
(Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.P[2] +=
+ pj->gradients.P[2] +=
(Wi[4] - Wj[4]) * wj *
(Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
} else {
/* SPH gradients */
- 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;
+ pj->gradients.v[0][1] -= wj_dx * dx[1] * (pi->v[0] - pj->v[0]) * r_inv;
+ pj->gradients.v[0][2] -= wj_dx * dx[2] * (pi->v[0] - pj->v[0]) * r_inv;
+
+ pj->gradients.v[1][0] -= wj_dx * dx[0] * (pi->v[1] - pj->v[1]) * r_inv;
+ pj->gradients.v[1][1] -= wj_dx * dx[1] * (pi->v[1] - pj->v[1]) * r_inv;
+ pj->gradients.v[1][2] -= wj_dx * dx[2] * (pi->v[1] - pj->v[1]) * r_inv;
+ pj->gradients.v[2][0] -= wj_dx * dx[0] * (pi->v[2] - pj->v[2]) * r_inv;
+ pj->gradients.v[2][1] -= wj_dx * dx[1] * (pi->v[2] - pj->v[2]) * r_inv;
+ pj->gradients.v[2][2] -= wj_dx * dx[2] * (pi->v[2] - pj->v[2]) * r_inv;
+
+ pj->gradients.P[0] -= wj_dx * dx[0] * (pi->P - pj->P) * r_inv;
+ pj->gradients.P[1] -= wj_dx * dx[1] * (pi->P - pj->P) * r_inv;
+ pj->gradients.P[2] -= wj_dx * dx[2] * (pi->P - pj->P) * r_inv;
}
hydro_slope_limit_cell_collect(pj, pi, r);
@@ -315,16 +285,16 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
Bi[k][l] = pi->geometry.matrix_E[k][l];
}
}
- Wi[0] = pi->primitives.rho;
+ Wi[0] = pi->rho;
Wi[1] = pi->v[0];
Wi[2] = pi->v[1];
Wi[3] = pi->v[2];
- Wi[4] = pi->primitives.P;
- Wj[0] = pj->primitives.rho;
+ Wi[4] = pi->P;
+ Wj[0] = pj->rho;
Wj[1] = pj->v[0];
Wj[2] = pj->v[1];
Wj[3] = pj->v[2];
- Wj[4] = pj->primitives.P;
+ Wj[4] = pj->P;
/* Compute kernel of pi. */
float wi, wi_dx;
@@ -332,94 +302,79 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
- if (pi->density.wcorr > const_gizmo_min_wcorr) {
+ if (pi->geometry.wcorr > const_gizmo_min_wcorr) {
/* Compute gradients for pi */
/* there is a sign difference w.r.t. eqn. (6) because of the inverse
* definition of dx */
- pi->primitives.gradients.rho[0] +=
+ pi->gradients.rho[0] +=
(Wi[0] - Wj[0]) * wi *
(Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.rho[1] +=
+ pi->gradients.rho[1] +=
(Wi[0] - Wj[0]) * wi *
(Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.rho[2] +=
+ pi->gradients.rho[2] +=
(Wi[0] - Wj[0]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->primitives.gradients.v[0][0] +=
+ pi->gradients.v[0][0] +=
(Wi[1] - Wj[1]) * wi *
(Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.v[0][1] +=
+ pi->gradients.v[0][1] +=
(Wi[1] - Wj[1]) * wi *
(Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.v[0][2] +=
+ pi->gradients.v[0][2] +=
(Wi[1] - Wj[1]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->primitives.gradients.v[1][0] +=
+ pi->gradients.v[1][0] +=
(Wi[2] - Wj[2]) * wi *
(Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.v[1][1] +=
+ pi->gradients.v[1][1] +=
(Wi[2] - Wj[2]) * wi *
(Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.v[1][2] +=
+ pi->gradients.v[1][2] +=
(Wi[2] - Wj[2]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->primitives.gradients.v[2][0] +=
+ pi->gradients.v[2][0] +=
(Wi[3] - Wj[3]) * wi *
(Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.v[2][1] +=
+ pi->gradients.v[2][1] +=
(Wi[3] - Wj[3]) * wi *
(Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.v[2][2] +=
+ pi->gradients.v[2][2] +=
(Wi[3] - Wj[3]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->primitives.gradients.P[0] +=
+ pi->gradients.P[0] +=
(Wi[4] - Wj[4]) * wi *
(Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.P[1] +=
+ pi->gradients.P[1] +=
(Wi[4] - Wj[4]) * wi *
(Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.P[2] +=
+ pi->gradients.P[2] +=
(Wi[4] - Wj[4]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
} else {
/* Gradient matrix is not well-behaved, switch to SPH gradients */
- 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);
@@ -440,46 +395,46 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
const float ihdim = pow_dimension(h_inv);
const float ihdimp1 = pow_dimension_plus_one(h_inv);
- if (p->density.wcorr > const_gizmo_min_wcorr) {
- p->primitives.gradients.rho[0] *= ihdim;
- p->primitives.gradients.rho[1] *= ihdim;
- p->primitives.gradients.rho[2] *= ihdim;
-
- p->primitives.gradients.v[0][0] *= ihdim;
- p->primitives.gradients.v[0][1] *= ihdim;
- p->primitives.gradients.v[0][2] *= ihdim;
- p->primitives.gradients.v[1][0] *= ihdim;
- p->primitives.gradients.v[1][1] *= ihdim;
- p->primitives.gradients.v[1][2] *= ihdim;
- p->primitives.gradients.v[2][0] *= ihdim;
- p->primitives.gradients.v[2][1] *= ihdim;
- p->primitives.gradients.v[2][2] *= ihdim;
-
- p->primitives.gradients.P[0] *= ihdim;
- p->primitives.gradients.P[1] *= ihdim;
- p->primitives.gradients.P[2] *= ihdim;
+ if (p->geometry.wcorr > const_gizmo_min_wcorr) {
+ p->gradients.rho[0] *= ihdim;
+ p->gradients.rho[1] *= ihdim;
+ p->gradients.rho[2] *= ihdim;
+
+ p->gradients.v[0][0] *= ihdim;
+ p->gradients.v[0][1] *= ihdim;
+ p->gradients.v[0][2] *= ihdim;
+ p->gradients.v[1][0] *= ihdim;
+ p->gradients.v[1][1] *= ihdim;
+ p->gradients.v[1][2] *= ihdim;
+ p->gradients.v[2][0] *= ihdim;
+ p->gradients.v[2][1] *= ihdim;
+ p->gradients.v[2][2] *= ihdim;
+
+ p->gradients.P[0] *= ihdim;
+ p->gradients.P[1] *= ihdim;
+ p->gradients.P[2] *= ihdim;
} else {
/* finalize gradients by multiplying with volume */
- p->primitives.gradients.rho[0] *= ihdimp1 * volume;
- p->primitives.gradients.rho[1] *= ihdimp1 * volume;
- p->primitives.gradients.rho[2] *= ihdimp1 * volume;
-
- p->primitives.gradients.v[0][0] *= ihdimp1 * volume;
- p->primitives.gradients.v[0][1] *= ihdimp1 * volume;
- p->primitives.gradients.v[0][2] *= ihdimp1 * volume;
-
- p->primitives.gradients.v[1][0] *= ihdimp1 * volume;
- p->primitives.gradients.v[1][1] *= ihdimp1 * volume;
- p->primitives.gradients.v[1][2] *= ihdimp1 * volume;
- p->primitives.gradients.v[2][0] *= ihdimp1 * volume;
- p->primitives.gradients.v[2][1] *= ihdimp1 * volume;
- p->primitives.gradients.v[2][2] *= ihdimp1 * volume;
-
- p->primitives.gradients.P[0] *= ihdimp1 * volume;
- p->primitives.gradients.P[1] *= ihdimp1 * volume;
- p->primitives.gradients.P[2] *= ihdimp1 * volume;
+ p->gradients.rho[0] *= ihdimp1 * volume;
+ p->gradients.rho[1] *= ihdimp1 * volume;
+ p->gradients.rho[2] *= ihdimp1 * volume;
+
+ p->gradients.v[0][0] *= ihdimp1 * volume;
+ p->gradients.v[0][1] *= ihdimp1 * volume;
+ p->gradients.v[0][2] *= ihdimp1 * volume;
+
+ p->gradients.v[1][0] *= ihdimp1 * volume;
+ p->gradients.v[1][1] *= ihdimp1 * volume;
+ p->gradients.v[1][2] *= ihdimp1 * volume;
+ p->gradients.v[2][0] *= ihdimp1 * volume;
+ p->gradients.v[2][1] *= ihdimp1 * volume;
+ p->gradients.v[2][2] *= ihdimp1 * volume;
+
+ p->gradients.P[0] *= ihdimp1 * volume;
+ p->gradients.P[1] *= ihdimp1 * volume;
+ p->gradients.P[2] *= ihdimp1 * volume;
}
hydro_slope_limit_cell(p);
diff --git a/src/hydro/GizmoMFM/hydro_gradients_sph.h b/src/hydro/GizmoMFM/hydro_gradients_sph.h
index 1960feef989f28835323faca1eb996e0e3272243..b4b6593a74dcea38a81d7e7b65eb0c5209216f3d 100644
--- a/src/hydro/GizmoMFM/hydro_gradients_sph.h
+++ b/src/hydro/GizmoMFM/hydro_gradients_sph.h
@@ -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;
+ pj->gradients.v[0][1] -= wj_dx * dx[1] * (pi->v[0] - pj->v[0]) * r_inv;
+ pj->gradients.v[0][2] -= wj_dx * dx[2] * (pi->v[0] - pj->v[0]) * r_inv;
+
+ pj->gradients.v[1][0] -= wj_dx * dx[0] * (pi->v[1] - pj->v[1]) * r_inv;
+ pj->gradients.v[1][1] -= wj_dx * dx[1] * (pi->v[1] - pj->v[1]) * r_inv;
+ pj->gradients.v[1][2] -= wj_dx * dx[2] * (pi->v[1] - pj->v[1]) * r_inv;
+ pj->gradients.v[2][0] -= wj_dx * dx[0] * (pi->v[2] - pj->v[2]) * r_inv;
+ pj->gradients.v[2][1] -= wj_dx * dx[1] * (pi->v[2] - pj->v[2]) * r_inv;
+ pj->gradients.v[2][2] -= wj_dx * dx[2] * (pi->v[2] - pj->v[2]) * r_inv;
+
+ pj->gradients.P[0] -= wj_dx * dx[0] * (pi->P - pj->P) * r_inv;
+ pj->gradients.P[1] -= wj_dx * dx[1] * (pi->P - pj->P) * r_inv;
+ pj->gradients.P[2] -= wj_dx * dx[2] * (pi->P - pj->P) * r_inv;
hydro_slope_limit_cell_collect(pj, pi, r);
}
@@ -178,40 +148,25 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
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);
}
@@ -230,25 +185,25 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
const float volume = p->geometry.volume;
/* finalize gradients by multiplying with volume */
- p->primitives.gradients.rho[0] *= ihdimp1 * volume;
- p->primitives.gradients.rho[1] *= ihdimp1 * volume;
- p->primitives.gradients.rho[2] *= ihdimp1 * volume;
+ p->gradients.rho[0] *= ihdimp1 * volume;
+ p->gradients.rho[1] *= ihdimp1 * volume;
+ p->gradients.rho[2] *= ihdimp1 * volume;
- p->primitives.gradients.v[0][0] *= ihdimp1 * volume;
- p->primitives.gradients.v[0][1] *= ihdimp1 * volume;
- p->primitives.gradients.v[0][2] *= ihdimp1 * volume;
+ p->gradients.v[0][0] *= ihdimp1 * volume;
+ p->gradients.v[0][1] *= ihdimp1 * volume;
+ p->gradients.v[0][2] *= ihdimp1 * volume;
- p->primitives.gradients.v[1][0] *= ihdimp1 * volume;
- p->primitives.gradients.v[1][1] *= ihdimp1 * volume;
- p->primitives.gradients.v[1][2] *= ihdimp1 * volume;
+ p->gradients.v[1][0] *= ihdimp1 * volume;
+ p->gradients.v[1][1] *= ihdimp1 * volume;
+ p->gradients.v[1][2] *= ihdimp1 * volume;
- p->primitives.gradients.v[2][0] *= ihdimp1 * volume;
- p->primitives.gradients.v[2][1] *= ihdimp1 * volume;
- p->primitives.gradients.v[2][2] *= ihdimp1 * volume;
+ p->gradients.v[2][0] *= ihdimp1 * volume;
+ p->gradients.v[2][1] *= ihdimp1 * volume;
+ p->gradients.v[2][2] *= ihdimp1 * volume;
- p->primitives.gradients.P[0] *= ihdimp1 * volume;
- p->primitives.gradients.P[1] *= ihdimp1 * volume;
- p->primitives.gradients.P[2] *= ihdimp1 * volume;
+ p->gradients.P[0] *= ihdimp1 * volume;
+ p->gradients.P[1] *= ihdimp1 * volume;
+ p->gradients.P[2] *= ihdimp1 * volume;
hydro_slope_limit_cell(p);
}
diff --git a/src/hydro/GizmoMFM/hydro_iact.h b/src/hydro/GizmoMFM/hydro_iact.h
index 40cd749d7b54db9a4a7db44b806689c561f353d1..947dd3de01d80289cde4d136600a1991339fcf65 100644
--- a/src/hydro/GizmoMFM/hydro_iact.h
+++ b/src/hydro/GizmoMFM/hydro_iact.h
@@ -238,16 +238,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
const float Vi = pi->geometry.volume;
const float Vj = pj->geometry.volume;
float Wi[5], Wj[5];
- Wi[0] = pi->primitives.rho;
+ Wi[0] = pi->rho;
Wi[1] = pi->v[0];
Wi[2] = pi->v[1];
Wi[3] = pi->v[2];
- Wi[4] = pi->primitives.P;
- Wj[0] = pj->primitives.rho;
+ Wi[4] = pi->P;
+ Wj[0] = pj->rho;
Wj[1] = pj->v[0];
Wj[2] = pj->v[1];
Wj[3] = pj->v[2];
- Wj[4] = pj->primitives.P;
+ Wj[4] = pj->P;
/* calculate the maximal signal velocity */
float vmax;
@@ -298,17 +298,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
const float wi_dr = hidp1 * wi_dx;
const float wj_dr = hjdp1 * wj_dx;
dvdr *= r_inv;
- if (pj->primitives.rho > 0.)
- pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
- if (mode == 1 && pi->primitives.rho > 0.)
- pj->force.h_dt -= pi->conserved.mass * dvdr / pi->primitives.rho * wj_dr;
+ if (pj->rho > 0.)
+ pi->force.h_dt -= pj->conserved.mass * dvdr / pj->rho * wi_dr;
+ if (mode == 1 && pi->rho > 0.)
+ pj->force.h_dt -= pi->conserved.mass * dvdr / pi->rho * wj_dr;
/* Compute (square of) area */
/* eqn. (7) */
float Anorm2 = 0.0f;
float A[3];
- if (pi->density.wcorr > const_gizmo_min_wcorr &&
- pj->density.wcorr > const_gizmo_min_wcorr) {
+ if (pi->geometry.wcorr > const_gizmo_min_wcorr &&
+ pj->geometry.wcorr > const_gizmo_min_wcorr) {
/* in principle, we use Vi and Vj as weights for the left and right
contributions to the generalized surface vector.
However, if Vi and Vj are very different (because they have very
@@ -418,15 +418,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* We shamelessly exploit the fact that the mass flux is zero and omit all
terms involving it */
/* eqn. (16) */
- pi->conserved.flux.momentum[0] -= totflux[1];
- pi->conserved.flux.momentum[1] -= totflux[2];
- pi->conserved.flux.momentum[2] -= totflux[3];
- pi->conserved.flux.energy -= totflux[4];
+ pi->flux.momentum[0] -= totflux[1];
+ pi->flux.momentum[1] -= totflux[2];
+ pi->flux.momentum[2] -= totflux[3];
+ pi->flux.energy -= totflux[4];
#ifndef GIZMO_TOTAL_ENERGY
- pi->conserved.flux.energy += totflux[1] * pi->v[0];
- pi->conserved.flux.energy += totflux[2] * pi->v[1];
- pi->conserved.flux.energy += totflux[3] * pi->v[2];
+ pi->flux.energy += totflux[1] * pi->v[0];
+ pi->flux.energy += totflux[2] * pi->v[1];
+ pi->flux.energy += totflux[3] * pi->v[2];
#endif
/* Note that this used to be much more complicated in early implementations of
@@ -435,15 +435,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
* conservation anymore and just assume the current fluxes are representative
* for the flux over the entire time step. */
if (mode == 1) {
- pj->conserved.flux.momentum[0] += totflux[1];
- pj->conserved.flux.momentum[1] += totflux[2];
- pj->conserved.flux.momentum[2] += totflux[3];
- pj->conserved.flux.energy += totflux[4];
+ pj->flux.momentum[0] += totflux[1];
+ pj->flux.momentum[1] += totflux[2];
+ pj->flux.momentum[2] += totflux[3];
+ pj->flux.energy += totflux[4];
#ifndef GIZMO_TOTAL_ENERGY
- pj->conserved.flux.energy -= totflux[1] * pj->v[0];
- pj->conserved.flux.energy -= totflux[2] * pj->v[1];
- pj->conserved.flux.energy -= totflux[3] * pj->v[2];
+ pj->flux.energy -= totflux[1] * pj->v[0];
+ pj->flux.energy -= totflux[2] * pj->v[1];
+ pj->flux.energy -= totflux[3] * pj->v[2];
#endif
}
}
diff --git a/src/hydro/GizmoMFM/hydro_io.h b/src/hydro/GizmoMFM/hydro_io.h
index 59d579f70cd4aedc728dbf42038eff78d4c507d5..1f956edf3fdc31990c6aba254603ea69a98238eb 100644
--- a/src/hydro/GizmoMFM/hydro_io.h
+++ b/src/hydro/GizmoMFM/hydro_io.h
@@ -63,7 +63,7 @@ INLINE static void hydro_read_particles(struct part* parts,
list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
UNIT_CONV_ACCELERATION, parts, a_hydro);
list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
- UNIT_CONV_DENSITY, parts, primitives.rho);
+ UNIT_CONV_DENSITY, parts, rho);
}
/**
@@ -203,12 +203,12 @@ INLINE static void hydro_write_particles(const struct part* parts,
parts, xparts, convert_u);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
- list[6] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
- primitives.rho);
+ list[6] =
+ io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
list[7] = io_make_output_field_convert_part(
"Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, xparts, convert_A);
- list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
- parts, primitives.P);
+ list[8] =
+ io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, P);
list[9] = io_make_output_field_convert_part(
"TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, xparts, convert_Etot);
diff --git a/src/hydro/GizmoMFM/hydro_part.h b/src/hydro/GizmoMFM/hydro_part.h
index 467872efa4464b6cecfb4dbc5cdb1aab9f36e203..0055d7d86a35746a8ba90015b3a6986f8ddb5f9f 100644
--- a/src/hydro/GizmoMFM/hydro_part.h
+++ b/src/hydro/GizmoMFM/hydro_part.h
@@ -63,28 +63,23 @@ struct part {
/* Particle smoothing length. */
float h;
- /* The primitive hydrodynamical variables. */
- struct {
-
- /* Density. */
- float rho;
+ /* Density. */
+ float rho;
- /* Pressure. */
- float P;
+ /* Pressure. */
+ float P;
- /* Gradients of the primitive variables. */
+ union {
+ /* Quantities used during the volume (=density) loop. */
struct {
- /* Density gradients. */
- float rho[3];
-
- /* Fluid velocity gradients. */
- float v[3][3];
+ /* Derivative of particle number density. */
+ float wcount_dh;
- /* Pressure gradients. */
- float P[3];
+ /* Particle number density. */
+ float wcount;
- } gradients;
+ } density;
/* Quantities needed by the slope limiter. */
struct {
@@ -103,7 +98,56 @@ struct part {
} limiter;
- } primitives;
+ struct {
+ /* Fluxes. */
+ struct {
+
+ /* No mass flux, since it is always zero. */
+
+ /* Momentum flux. */
+ float momentum[3];
+
+ /* Energy flux. */
+ float energy;
+
+ } flux;
+
+ /* Variables used for timestep calculation. */
+ struct {
+
+ /* Maximum signal velocity among all the neighbours of the particle. The
+ * signal velocity encodes information about the relative fluid
+ * velocities
+ * AND particle velocities of the neighbour and this particle, as well
+ * as
+ * the sound speed of both particles. */
+ float vmax;
+
+ } timestepvars;
+
+ /* Quantities used during the force loop. */
+ struct {
+
+ /* Needed to drift the primitive variables. */
+ float h_dt;
+
+ } force;
+ };
+ };
+
+ /* Gradients of the primitive variables. */
+ struct {
+
+ /* Density gradients. */
+ float rho[3];
+
+ /* Fluid velocity gradients. */
+ float v[3][3];
+
+ /* Pressure gradients. */
+ float P[3];
+
+ } gradients;
/* The conserved hydrodynamical variables. */
struct {
@@ -117,19 +161,6 @@ struct part {
/* Fluid thermal energy (not per unit mass!). */
float energy;
- /* Fluxes. */
- struct {
-
- /* No mass flux, since it is always zero. */
-
- /* Momentum flux. */
- float momentum[3];
-
- /* Energy flux. */
- float energy;
-
- } flux;
-
} conserved;
/* Geometrical quantities used for hydro. */
@@ -145,40 +176,10 @@ struct part {
/* Centroid of the "cell". */
float centroid[3];
- } geometry;
-
- /* Variables used for timestep calculation. */
- struct {
-
- /* Maximum signal velocity among all the neighbours of the particle. The
- * signal velocity encodes information about the relative fluid velocities
- * AND particle velocities of the neighbour and this particle, as well as
- * the sound speed of both particles. */
- float vmax;
-
- } timestepvars;
-
- /* Quantities used during the volume (=density) loop. */
- struct {
-
- /* Derivative of particle number density. */
- float wcount_dh;
-
- /* Particle number density. */
- float wcount;
-
/* Correction factor for wcount. */
float wcorr;
- } density;
-
- /* Quantities used during the force loop. */
- struct {
-
- /* Needed to drift the primitive variables. */
- float h_dt;
-
- } force;
+ } geometry;
/* Chemistry information */
struct chemistry_part_data chemistry_data;
diff --git a/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h b/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h
index c03a626c688ccf51e8d8cc02dabf261ef52b7806..a33d03dee8d54f6044b2080d095e34e7df15cfa0 100644
--- a/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h
+++ b/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h
@@ -29,18 +29,18 @@
__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
struct part* p) {
- p->primitives.limiter.rho[0] = FLT_MAX;
- p->primitives.limiter.rho[1] = -FLT_MAX;
- p->primitives.limiter.v[0][0] = FLT_MAX;
- p->primitives.limiter.v[0][1] = -FLT_MAX;
- p->primitives.limiter.v[1][0] = FLT_MAX;
- p->primitives.limiter.v[1][1] = -FLT_MAX;
- p->primitives.limiter.v[2][0] = FLT_MAX;
- p->primitives.limiter.v[2][1] = -FLT_MAX;
- p->primitives.limiter.P[0] = FLT_MAX;
- p->primitives.limiter.P[1] = -FLT_MAX;
-
- p->primitives.limiter.maxr = -FLT_MAX;
+ p->limiter.rho[0] = FLT_MAX;
+ p->limiter.rho[1] = -FLT_MAX;
+ p->limiter.v[0][0] = FLT_MAX;
+ p->limiter.v[0][1] = -FLT_MAX;
+ p->limiter.v[1][0] = FLT_MAX;
+ p->limiter.v[1][1] = -FLT_MAX;
+ p->limiter.v[2][0] = FLT_MAX;
+ p->limiter.v[2][1] = -FLT_MAX;
+ p->limiter.P[0] = FLT_MAX;
+ p->limiter.P[1] = -FLT_MAX;
+
+ p->limiter.maxr = -FLT_MAX;
}
/**
@@ -56,30 +56,20 @@ hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
/* basic slope limiter: collect the maximal and the minimal value for the
* primitive variables among the ngbs */
- pi->primitives.limiter.rho[0] =
- min(pj->primitives.rho, pi->primitives.limiter.rho[0]);
- pi->primitives.limiter.rho[1] =
- max(pj->primitives.rho, pi->primitives.limiter.rho[1]);
-
- pi->primitives.limiter.v[0][0] =
- min(pj->v[0], pi->primitives.limiter.v[0][0]);
- pi->primitives.limiter.v[0][1] =
- max(pj->v[0], pi->primitives.limiter.v[0][1]);
- pi->primitives.limiter.v[1][0] =
- min(pj->v[1], pi->primitives.limiter.v[1][0]);
- pi->primitives.limiter.v[1][1] =
- max(pj->v[1], pi->primitives.limiter.v[1][1]);
- pi->primitives.limiter.v[2][0] =
- min(pj->v[2], pi->primitives.limiter.v[2][0]);
- pi->primitives.limiter.v[2][1] =
- max(pj->v[2], pi->primitives.limiter.v[2][1]);
-
- pi->primitives.limiter.P[0] =
- min(pj->primitives.P, pi->primitives.limiter.P[0]);
- pi->primitives.limiter.P[1] =
- max(pj->primitives.P, pi->primitives.limiter.P[1]);
-
- pi->primitives.limiter.maxr = max(r, pi->primitives.limiter.maxr);
+ pi->limiter.rho[0] = min(pj->rho, pi->limiter.rho[0]);
+ pi->limiter.rho[1] = max(pj->rho, pi->limiter.rho[1]);
+
+ pi->limiter.v[0][0] = min(pj->v[0], pi->limiter.v[0][0]);
+ pi->limiter.v[0][1] = max(pj->v[0], pi->limiter.v[0][1]);
+ pi->limiter.v[1][0] = min(pj->v[1], pi->limiter.v[1][0]);
+ pi->limiter.v[1][1] = max(pj->v[1], pi->limiter.v[1][1]);
+ pi->limiter.v[2][0] = min(pj->v[2], pi->limiter.v[2][0]);
+ pi->limiter.v[2][1] = max(pj->v[2], pi->limiter.v[2][1]);
+
+ pi->limiter.P[0] = min(pj->P, pi->limiter.P[0]);
+ pi->limiter.P[1] = max(pj->P, pi->limiter.P[1]);
+
+ pi->limiter.maxr = max(r, pi->limiter.maxr);
}
/**
@@ -92,94 +82,94 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
float gradtrue, gradrho[3], gradv[3][3], gradP[3];
- gradrho[0] = p->primitives.gradients.rho[0];
- gradrho[1] = p->primitives.gradients.rho[1];
- gradrho[2] = p->primitives.gradients.rho[2];
+ gradrho[0] = p->gradients.rho[0];
+ gradrho[1] = p->gradients.rho[1];
+ gradrho[2] = p->gradients.rho[2];
- gradv[0][0] = p->primitives.gradients.v[0][0];
- gradv[0][1] = p->primitives.gradients.v[0][1];
- gradv[0][2] = p->primitives.gradients.v[0][2];
+ gradv[0][0] = p->gradients.v[0][0];
+ gradv[0][1] = p->gradients.v[0][1];
+ gradv[0][2] = p->gradients.v[0][2];
- gradv[1][0] = p->primitives.gradients.v[1][0];
- gradv[1][1] = p->primitives.gradients.v[1][1];
- gradv[1][2] = p->primitives.gradients.v[1][2];
+ gradv[1][0] = p->gradients.v[1][0];
+ gradv[1][1] = p->gradients.v[1][1];
+ gradv[1][2] = p->gradients.v[1][2];
- gradv[2][0] = p->primitives.gradients.v[2][0];
- gradv[2][1] = p->primitives.gradients.v[2][1];
- gradv[2][2] = p->primitives.gradients.v[2][2];
+ gradv[2][0] = p->gradients.v[2][0];
+ gradv[2][1] = p->gradients.v[2][1];
+ gradv[2][2] = p->gradients.v[2][2];
- gradP[0] = p->primitives.gradients.P[0];
- gradP[1] = p->primitives.gradients.P[1];
- gradP[2] = p->primitives.gradients.P[2];
+ gradP[0] = p->gradients.P[0];
+ gradP[1] = p->gradients.P[1];
+ gradP[2] = p->gradients.P[2];
gradtrue = sqrtf(gradrho[0] * gradrho[0] + gradrho[1] * gradrho[1] +
gradrho[2] * gradrho[2]);
if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- const float gradmax = p->primitives.limiter.rho[1] - p->primitives.rho;
- const float gradmin = p->primitives.rho - p->primitives.limiter.rho[0];
+ gradtrue *= p->limiter.maxr;
+ const float gradmax = p->limiter.rho[1] - p->rho;
+ const float gradmin = p->rho - p->limiter.rho[0];
const float gradtrue_inv = 1.f / gradtrue;
const float alpha =
min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
- p->primitives.gradients.rho[0] *= alpha;
- p->primitives.gradients.rho[1] *= alpha;
- p->primitives.gradients.rho[2] *= alpha;
+ p->gradients.rho[0] *= alpha;
+ p->gradients.rho[1] *= alpha;
+ p->gradients.rho[2] *= alpha;
}
gradtrue = sqrtf(gradv[0][0] * gradv[0][0] + gradv[0][1] * gradv[0][1] +
gradv[0][2] * gradv[0][2]);
if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- const float gradmax = p->primitives.limiter.v[0][1] - p->v[0];
- const float gradmin = p->v[0] - p->primitives.limiter.v[0][0];
+ gradtrue *= p->limiter.maxr;
+ const float gradmax = p->limiter.v[0][1] - p->v[0];
+ const float gradmin = p->v[0] - p->limiter.v[0][0];
const float gradtrue_inv = 1.f / gradtrue;
const float alpha =
min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
- p->primitives.gradients.v[0][0] *= alpha;
- p->primitives.gradients.v[0][1] *= alpha;
- p->primitives.gradients.v[0][2] *= alpha;
+ p->gradients.v[0][0] *= alpha;
+ p->gradients.v[0][1] *= alpha;
+ p->gradients.v[0][2] *= alpha;
}
gradtrue = sqrtf(gradv[1][0] * gradv[1][0] + gradv[1][1] * gradv[1][1] +
gradv[1][2] * gradv[1][2]);
if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- const float gradmax = p->primitives.limiter.v[1][1] - p->v[1];
- const float gradmin = p->v[1] - p->primitives.limiter.v[1][0];
+ gradtrue *= p->limiter.maxr;
+ const float gradmax = p->limiter.v[1][1] - p->v[1];
+ const float gradmin = p->v[1] - p->limiter.v[1][0];
const float gradtrue_inv = 1.f / gradtrue;
const float alpha =
min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
- p->primitives.gradients.v[1][0] *= alpha;
- p->primitives.gradients.v[1][1] *= alpha;
- p->primitives.gradients.v[1][2] *= alpha;
+ p->gradients.v[1][0] *= alpha;
+ p->gradients.v[1][1] *= alpha;
+ p->gradients.v[1][2] *= alpha;
}
gradtrue = sqrtf(gradv[2][0] * gradv[2][0] + gradv[2][1] * gradv[2][1] +
gradv[2][2] * gradv[2][2]);
if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- const float gradmax = p->primitives.limiter.v[2][1] - p->v[2];
- const float gradmin = p->v[2] - p->primitives.limiter.v[2][0];
+ gradtrue *= p->limiter.maxr;
+ const float gradmax = p->limiter.v[2][1] - p->v[2];
+ const float gradmin = p->v[2] - p->limiter.v[2][0];
const float gradtrue_inv = 1.f / gradtrue;
const float alpha =
min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
- p->primitives.gradients.v[2][0] *= alpha;
- p->primitives.gradients.v[2][1] *= alpha;
- p->primitives.gradients.v[2][2] *= alpha;
+ p->gradients.v[2][0] *= alpha;
+ p->gradients.v[2][1] *= alpha;
+ p->gradients.v[2][2] *= alpha;
}
gradtrue =
sqrtf(gradP[0] * gradP[0] + gradP[1] * gradP[1] + gradP[2] * gradP[2]);
if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- const float gradmax = p->primitives.limiter.P[1] - p->primitives.P;
- const float gradmin = p->primitives.P - p->primitives.limiter.P[0];
+ gradtrue *= p->limiter.maxr;
+ const float gradmax = p->limiter.P[1] - p->P;
+ const float gradmin = p->P - p->limiter.P[0];
const float gradtrue_inv = 1.f / gradtrue;
const float alpha =
min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
- p->primitives.gradients.P[0] *= alpha;
- p->primitives.gradients.P[1] *= alpha;
- p->primitives.gradients.P[2] *= alpha;
+ p->gradients.P[0] *= alpha;
+ p->gradients.P[1] *= alpha;
+ p->gradients.P[2] *= alpha;
}
}