diff --git a/examples/Planetary/DemoImpactInitCond/plot_profiles.py b/examples/Planetary/DemoImpactInitCond/plot_profiles.py
index a03fa9a7c8f065b49d297241115408903b428d25..e60b2e3201495017a86642bd498bef98c1efc6e8 100755
--- a/examples/Planetary/DemoImpactInitCond/plot_profiles.py
+++ b/examples/Planetary/DemoImpactInitCond/plot_profiles.py
@@ -22,7 +22,7 @@ Note that, for standard SPH hydro schemes, especially at low resolution, the
standard issues that arise at discontinuities in material and density lead the
SPH particles to settle at slightly different densities near discontinuties.
For more info and to explore ways to resolve these issues, check out e.g.
-Sandnes et al. (2024), Ruiz-Bonilla el al. (2022), and Kegerreis et al. (2019).
+Sandnes et al. (2025), Ruiz-Bonilla el al. (2022), and Kegerreis et al. (2019).
The overall profile of the settled SPH planet should still align with the input.
"""
diff --git a/src/hydro.h b/src/hydro.h
index c6fd08cc47f55384f4f622819f0247a51ac449c1..95860d625b5437f7e09eab855e7da5cf076eb289 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -74,7 +74,7 @@
#elif defined(REMIX_SPH)
#include "./hydro/REMIX/hydro.h"
#include "./hydro/REMIX/hydro_iact.h"
-#define SPH_IMPLEMENTATION "SPH with multiple materials and methods"
+#define SPH_IMPLEMENTATION "REMIX (Sandnes+ 2025)"
#elif defined(SPHENIX_SPH)
#include "./hydro/SPHENIX/hydro.h"
#include "./hydro/SPHENIX/hydro_iact.h"
diff --git a/src/hydro/REMIX/hydro.h b/src/hydro/REMIX/hydro.h
index 21508917e48dc59809969faf2ae7d91d3b4a9112..4498f987009ee46cc3c04b93e103d15a59658af2 100644
--- a/src/hydro/REMIX/hydro.h
+++ b/src/hydro/REMIX/hydro.h
@@ -23,7 +23,7 @@
/**
* @file Planetary/hydro.h
- * @brief REMIX implementation of SPH (Sandnes et al. 2024)
+ * @brief REMIX implementation of SPH (Sandnes et al. 2025)
*/
#include "adiabatic_index.h"
@@ -493,7 +493,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
#endif
hydro_init_part_extra_kernel(p);
- hydro_init_part_extra_viscosity(p);
+ hydro_init_part_extra_visc_difn(p);
}
/**
@@ -535,7 +535,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->n_density *= h_inv_dim;
#endif
- hydro_end_density_extra_viscosity(p);
+ hydro_end_density_extra_visc_difn(p);
hydro_end_density_extra_kernel(p);
}
@@ -589,6 +589,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
const struct cosmology *cosmo, const struct hydro_props *hydro_props,
const struct pressure_floor_props *pressure_floor) {
+ /* Set p->is_h_max = 1 if smoothing length is h_max */
if (p->h > 0.999f * hydro_props->h_max) {
p->is_h_max = 1;
} else {
@@ -596,7 +597,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
}
hydro_prepare_gradient_extra_kernel(p);
- hydro_prepare_gradient_extra_viscosity(p);
+ hydro_prepare_gradient_extra_visc_difn(p);
}
/**
@@ -625,9 +626,9 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
struct part *p) {
hydro_end_gradient_extra_kernel(p);
- hydro_end_gradient_extra_viscosity(p);
+ hydro_end_gradient_extra_visc_difn(p);
- // Set the density to be used in the force loop to be the evolved density
+ /* Set the density to be used in the force loop to be the evolved density */
p->rho = p->rho_evol;
}
@@ -676,7 +677,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
float mod_curl_v = sqrtf(curl_v[0] * curl_v[0] + curl_v[1] * curl_v[1] +
curl_v[2] * curl_v[2]);
- // Balsara switch using normalised kernel gradients
+ /* Balsara switch using normalised kernel gradients (Sandnes+2025 Eqn. 34 with
+ * velocity gradients calculated by Eqn. 35) */
float balsara;
if (div_v == 0.f) {
balsara = 0.f;
@@ -691,6 +693,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
p->force.balsara = balsara;
+ /* Set eta_crit for arificial viscosity and diffusion slope limiters */
p->force.eta_crit = 1.f / hydro_props->eta_neighbours;
}
@@ -738,6 +741,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
/* Re-set the internal energy */
p->u = xp->u_full;
+ /* Re-set the density */
p->rho = xp->rho_evol_full;
p->rho_evol = xp->rho_evol_full;
@@ -790,6 +794,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
const float floor_rho = p->mass * kernel_root * h_inv_dim;
p->rho_evol = max(p->rho_evol, floor_rho);
+ p->rho = p->rho_evol;
/* Check against absolute minimum */
const float min_u =
@@ -804,15 +809,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
else
p->h *= expf(w1);
- /* Predict density */
- const float w2 = -hydro_dimension * w1;
- if (fabsf(w2) < 0.2f)
- p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
- else
- p->rho *= expf(w2);
-
- p->rho = p->rho_evol;
-
const float floor_u = FLT_MIN;
p->u = max(p->u, floor_u);
@@ -890,12 +886,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->u_dt = 0.f;
}
+ /* Integrate the density forward in time */
const float delta_rho = p->drho_dt * dt_therm;
+ /* Do not decrease the density by more than a factor of 2*/
xp->rho_evol_full =
max(xp->rho_evol_full + delta_rho, 0.5f * xp->rho_evol_full);
- /* Minimum SPH quantities */
+ /* Minimum SPH density */
const float h = p->h;
const float h_inv = 1.0f / h; /* 1/h */
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
diff --git a/src/hydro/REMIX/hydro_iact.h b/src/hydro/REMIX/hydro_iact.h
index fa04f2a30555dfda44e47ffb5727f860ca02bcf5..9e1b65794974e2ffc43f014345d744304615af05 100644
--- a/src/hydro/REMIX/hydro_iact.h
+++ b/src/hydro/REMIX/hydro_iact.h
@@ -23,7 +23,7 @@
/**
* @file Planetary/hydro_iact.h
- * @brief REMIX implementation of SPH (Sandnes et al. 2024)
+ * @brief REMIX implementation of SPH (Sandnes et al. 2025)
*/
#include "adiabatic_index.h"
@@ -95,7 +95,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
#endif
hydro_runner_iact_density_extra_kernel(pi, pj, dx, wi, wj, wi_dx, wj_dx);
- hydro_runner_iact_density_extra_viscosity(pi, pj, dx, wi, wj, wi_dx, wj_dx);
+ hydro_runner_iact_density_extra_visc_difn(pi, pj, dx, wi, wj, wi_dx, wj_dx);
}
/**
@@ -145,7 +145,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
#endif
hydro_runner_iact_nonsym_density_extra_kernel(pi, pj, dx, wi, wi_dx);
- hydro_runner_iact_nonsym_density_extra_viscosity(pi, pj, dx, wi, wi_dx);
+ hydro_runner_iact_nonsym_density_extra_visc_difn(pi, pj, dx, wi, wi_dx);
}
/**
@@ -183,7 +183,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
kernel_deval(uj, &wj, &wj_dx);
hydro_runner_iact_gradient_extra_kernel(pi, pj, dx, wi, wj, wi_dx, wj_dx);
- hydro_runner_iact_gradient_extra_viscosity(pi, pj, dx, wi, wj, wi_dx, wj_dx);
+ hydro_runner_iact_gradient_extra_visc_difn(pi, pj, dx, wi, wj, wi_dx, wj_dx);
}
/**
@@ -223,7 +223,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
hydro_runner_iact_nonsym_gradient_extra_kernel(pi, pj, dx, wi, wj, wi_dx,
wj_dx);
- hydro_runner_iact_nonsym_gradient_extra_viscosity(pi, pj, dx, wi, wi_dx);
+ hydro_runner_iact_nonsym_gradient_extra_visc_difn(pi, pj, dx, wi, wi_dx);
}
/**
@@ -250,7 +250,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
error("Inhibited pj in interaction function!");
#endif
- /* Get r and 1/r. */
+ /* Get r. */
const float r = sqrtf(r2);
/* Recover some data */
@@ -273,18 +273,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
float wj, wj_dx;
kernel_deval(xj, &wj, &wj_dx);
- // Linear-order reproducing kernel gradient term
+ /* Linear-order reproducing kernel gradient term (Sandnes+2025 Eqn. 28) */
float Gj[3], Gi[3], G_mean[3];
hydro_set_Gi_Gj_forceloop(Gi, Gj, pi, pj, dx, wi, wj, wi_dx, wj_dx);
for (int i = 0; i < 3; i++) {
+ /* Antisymmetric kernel grad term for conservation of momentum and energy */
G_mean[i] = 0.5f * (Gi[i] - Gj[i]);
}
- // Viscous pressures
+ /* Viscous pressures (Sandnes+2025 Eqn. 41) */
float Qi, Qj;
float visc_signal_velocity, difn_signal_velocity;
hydro_set_Qi_Qj(&Qi, &Qj, &visc_signal_velocity, &difn_signal_velocity, pi,
- pj, dx, a, H);
+ pj, dx);
/* Pressure terms to be used in evolution equations */
float P_i_term = pressurei / (pi->rho * pj->rho);
@@ -300,18 +301,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
mi * (P_i_term + P_j_term + Q_i_term + Q_j_term) * G_mean[i];
}
- // v_ij dot kernel gradient term
+ /* v_ij dot kernel gradient term */
const float dvdotG = (pi->v[0] - pj->v[0]) * G_mean[0] +
(pi->v[1] - pj->v[1]) * G_mean[1] +
(pi->v[2] - pj->v[2]) * G_mean[2];
- /* Get the time derivative for u, including the viscosity */
- float du_dt_i = (P_i_term + Q_i_term) * dvdotG;
- float du_dt_j = (P_j_term + Q_j_term) * dvdotG;
-
/* Internal energy time derivative */
- pi->u_dt += du_dt_i * mj;
- pj->u_dt += du_dt_j * mi;
+ pi->u_dt += mj * (P_i_term + Q_i_term) * dvdotG;
+ pj->u_dt += mi * (P_j_term + Q_j_term) * dvdotG;
+
+ /* Density time derivative */
+ pi->drho_dt += mj * (pi->rho / pj->rho) * dvdotG;
+ pj->drho_dt += mi * (pj->rho / pi->rho) * dvdotG;
/* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdotG / rhoj;
@@ -323,10 +324,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pi->force.v_sig = max(pi->force.v_sig, v_sig);
pj->force.v_sig = max(pj->force.v_sig, v_sig);
- pi->drho_dt += mj * (pi->rho / pj->rho) * dvdotG;
- pj->drho_dt += mi * (pj->rho / pi->rho) * dvdotG;
-
+ /* Only apply diffusion and normalising term if neither particle has h=h_max */
if ((!pi->is_h_max) && (!pj->is_h_max)) {
+ /* Calculate some quantities */
float mean_rho = 0.5f * (pi->rho + pj->rho);
float mean_balsara = 0.5f * (pi->force.balsara + pj->force.balsara);
float mod_G = sqrtf(G_mean[0] * G_mean[0] + G_mean[1] * G_mean[1] +
@@ -335,6 +335,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
(pi->v[1] - pj->v[1]) * (pi->v[1] - pj->v[1]) +
(pi->v[2] - pj->v[2]) * (pi->v[2] - pj->v[2]));
+ /* Add normalising term to density evolution (Sandnes+2025 Eqn. 51) */
const float alpha_norm = const_remix_norm_alpha;
float drho_dt_norm_and_difn_i =
alpha_norm * mj * v_sig_norm * pi->force.vac_switch *
@@ -343,19 +344,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
alpha_norm * mi * v_sig_norm * pj->force.vac_switch *
(pj->m0 * pj->rho_evol - pj->rho_evol) * mod_G / mean_rho;
- // Diffusion for same materials
+ /* Only include diffusion for same-material particle pair */
if (pi->mat_id == pj->mat_id) {
- // Diffusion parameters
+ /* Diffusion parameters */
const float a_difn_rho = const_remix_difn_a_rho;
const float b_difn_rho = const_remix_difn_b_rho;
const float a_difn_u = const_remix_difn_a_u;
const float b_difn_u = const_remix_difn_b_u;
- // ...
float utilde_i, utilde_j, rhotilde_i, rhotilde_j;
hydro_set_u_rho_difn(&utilde_i, &utilde_j, &rhotilde_i, &rhotilde_j, pi,
- pj, dx, a, H);
+ pj, dx);
float v_sig_difn = difn_signal_velocity;
+
+ /* Calculate artificial diffusion of internal energy (Sandnes+2025 Eqn. 42) */
float du_dt_difn_i = -(a_difn_u + b_difn_u * mean_balsara) * mj *
v_sig_difn * (utilde_i - utilde_j) * mod_G /
mean_rho;
@@ -363,11 +365,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
v_sig_difn * (utilde_j - utilde_i) * mod_G /
mean_rho;
- // ...
+ /* Add artificial diffusion to evolution of internal energy */
pi->u_dt += du_dt_difn_i;
pj->u_dt += du_dt_difn_j;
- // ...
+ /* Calculate artificial diffusion of density (Sandnes+2025 Eqn. 43) */
drho_dt_norm_and_difn_i += -(a_difn_rho + b_difn_rho * mean_balsara) *
mj * (pi->rho / pj->rho) * v_sig_difn *
(rhotilde_i - rhotilde_j) * mod_G / mean_rho;
@@ -376,6 +378,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
(rhotilde_j - rhotilde_i) * mod_G / mean_rho;
}
+ /* Add normalising term and artificial diffusion to evolution of density */
pi->drho_dt += drho_dt_norm_and_difn_i;
pj->drho_dt += drho_dt_norm_and_difn_j;
}
@@ -412,7 +415,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
error("Inhibited pj in interaction function!");
#endif
- /* Get r and 1/r. */
+ /* Get r. */
const float r = sqrtf(r2);
/* Recover some data */
@@ -433,18 +436,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
float wj, wj_dx;
kernel_deval(xj, &wj, &wj_dx);
- // Linear-order reproducing kernel gradient term
+ /* Linear-order reproducing kernel gradient term (Sandnes+2025 Eqn. 28) */
float Gj[3], Gi[3], G_mean[3];
hydro_set_Gi_Gj_forceloop(Gi, Gj, pi, pj, dx, wi, wj, wi_dx, wj_dx);
for (int i = 0; i < 3; i++) {
+ /* Antisymmetric kernel grad term for conservation of momentum and energy */
G_mean[i] = 0.5f * (Gi[i] - Gj[i]);
}
- // Viscous pressures
+ /* Viscous pressures (Sandnes+2025 Eqn. 41) */
float Qi, Qj;
float visc_signal_velocity, difn_signal_velocity;
hydro_set_Qi_Qj(&Qi, &Qj, &visc_signal_velocity, &difn_signal_velocity, pi,
- pj, dx, a, H);
+ pj, dx);
/* Pressure terms to be used in evolution equations */
float P_i_term = pressurei / (pi->rho * pj->rho);
@@ -458,16 +462,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
mj * (P_i_term + P_j_term + Q_i_term + Q_j_term) * G_mean[i];
}
- // v_ij dot kernel gradient term
+ /* v_ij dot kernel gradient term */
const float dvdotG = (pi->v[0] - pj->v[0]) * G_mean[0] +
(pi->v[1] - pj->v[1]) * G_mean[1] +
(pi->v[2] - pj->v[2]) * G_mean[2];
- /* Get the time derivative for u, including the viscosity */
- float du_dt_i = (P_i_term + Q_i_term) * dvdotG;
-
/* Internal energy time derivative */
- pi->u_dt += du_dt_i * mj;
+ pi->u_dt += mj * (P_i_term + Q_i_term) * dvdotG;
+
+ /* Density time derivative */
+ pi->drho_dt += mj * (pi->rho / pj->rho) * dvdotG;
/* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdotG / rhoj;
@@ -477,9 +481,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Update the signal velocity. */
pi->force.v_sig = max(pi->force.v_sig, v_sig);
- pi->drho_dt += mj * (pi->rho / pj->rho) * dvdotG;
-
+ /* Only apply diffusion and normalising term if neither particle has h=h_max */
if ((!pi->is_h_max) && (!pj->is_h_max)) {
+ /* Calculate some quantities */
float mean_rho = 0.5f * (pi->rho + pj->rho);
float mean_balsara = 0.5f * (pi->force.balsara + pj->force.balsara);
float mod_G = sqrtf(G_mean[0] * G_mean[0] + G_mean[1] * G_mean[1] +
@@ -489,37 +493,41 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
(pi->v[1] - pj->v[1]) * (pi->v[1] - pj->v[1]) +
(pi->v[2] - pj->v[2]) * (pi->v[2] - pj->v[2]));
+
+ /* Add normalising term to density evolution (Sandnes+2025 Eqn. 51) */
const float alpha_norm = const_remix_norm_alpha;
float drho_dt_norm_and_difn_i =
alpha_norm * mj * v_sig_norm * pi->force.vac_switch *
(pi->m0 * pi->rho_evol - pi->rho_evol) * mod_G / mean_rho;
- // Diffusion for same materials
+ /* Only include diffusion for same-material particle pair */
if (pi->mat_id == pj->mat_id) {
- // Diffusion parameters
+ /* Diffusion parameters */
const float a_difn_rho = const_remix_difn_a_rho;
const float b_difn_rho = const_remix_difn_b_rho;
const float a_difn_u = const_remix_difn_a_u;
const float b_difn_u = const_remix_difn_b_u;
- // ...
float utilde_i, utilde_j, rhotilde_i, rhotilde_j;
hydro_set_u_rho_difn(&utilde_i, &utilde_j, &rhotilde_i, &rhotilde_j, pi,
- pj, dx, a, H);
+ pj, dx);
float v_sig_difn = difn_signal_velocity;
+
+ /* Calculate artificial diffusion of internal energy (Sandnes+2025 Eqn. 42) */
float du_dt_difn_i = -(a_difn_u + b_difn_u * mean_balsara) * mj *
v_sig_difn * (utilde_i - utilde_j) * mod_G /
mean_rho;
- // ...
+ /* Add artificial diffusion to evolution of internal energy */
pi->u_dt += du_dt_difn_i;
- // ...
+ /* Calculate artificial diffusion of density (Sandnes+2025 Eqn. 43) */
drho_dt_norm_and_difn_i += -(a_difn_rho + b_difn_rho * mean_balsara) *
mj * (pi->rho / pj->rho) * v_sig_difn *
(rhotilde_i - rhotilde_j) * mod_G / mean_rho;
}
+ /* Add normalising term and artificial diffusion to evolution of density */
pi->drho_dt += drho_dt_norm_and_difn_i;
}
diff --git a/src/hydro/REMIX/hydro_io.h b/src/hydro/REMIX/hydro_io.h
index 6fdd3cae1823a672a1e15c1eabbd09a266ebd4a6..d26f224f8e649da5b677b661d7e101d82ad42cd5 100644
--- a/src/hydro/REMIX/hydro_io.h
+++ b/src/hydro/REMIX/hydro_io.h
@@ -23,7 +23,7 @@
/**
* @file Planetary/hydro_io.h
- * @brief REMIX implementation of SPH (Sandnes et al. 2024) i/o routines
+ * @brief REMIX implementation of SPH (Sandnes et al. 2025) i/o routines
*/
#include "adiabatic_index.h"
diff --git a/src/hydro/REMIX/hydro_kernels.h b/src/hydro/REMIX/hydro_kernels.h
index 45beda88cc6f0639d765389ce593a0e170f7ce02..dbb0e1c8bd7d97ec986a87f1a71d15ab32679833 100644
--- a/src/hydro/REMIX/hydro_kernels.h
+++ b/src/hydro/REMIX/hydro_kernels.h
@@ -32,12 +32,13 @@
/**
* @brief Prepares extra kernel parameters for a particle for the density
* calculation.
+ *
+ * @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_init_part_extra_kernel(
struct part *restrict p) {
p->m0 = 0.f;
-
p->grad_m0[0] = 0.f;
p->grad_m0[1] = 0.f;
p->grad_m0[2] = 0.f;
@@ -45,6 +46,14 @@ __attribute__((always_inline)) INLINE static void hydro_init_part_extra_kernel(
/**
* @brief Extra kernel density interaction between two particles
+ *
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param wi The value of the unmodified kernel function W(r, hi) * hi^d.
+ * @param wj The value of the unmodified kernel function W(r, hj) * hj^d.
+ * @param wi_dx The norm of the gradient of wi: dW(r, hi)/dr * hi^(d+1).
+ * @param wj_dx The norm of the gradient of wj: dW(r, hj)/dr * hj^(d+1).
*/
__attribute__((always_inline)) INLINE static void
hydro_runner_iact_density_extra_kernel(struct part *restrict pi,
@@ -57,7 +66,11 @@ hydro_runner_iact_density_extra_kernel(struct part *restrict pi,
const float r = sqrtf(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
const float r_inv = r ? 1.0f / r : 0.0f;
- /* Geometric moments and gradients (unmodified kernel, for normalisation) */
+ /* Geometric moments and gradients that use an unmodified kernel (Sandnes+2025
+ * Eqn. 50 and its gradient). Used in the normalising term (Eqn. 51) and in
+ * gradient estimates (using Eqn. 30) that are used for the calculation of grad-h
+ * terms (Eqn. 31) and in the artificial viscosity (Eqn. 35) and diffusion
+ * (Eqns. 46 and 47) schemes */
pi->m0 += pj->mass * wi / pj->rho_evol;
pj->m0 += pi->mass * wj / pi->rho_evol;
for (int i = 0; i < 3; i++) {
@@ -68,6 +81,12 @@ hydro_runner_iact_density_extra_kernel(struct part *restrict pi,
/**
* @brief Extra kernel density interaction between two particles (non-symmetric)
+ *
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param wi The value of the unmodified kernel function W(r, hi) * hi^d.
+ * @param wi_dx The norm of the gradient of wi: dW(r, hi)/dr * hi^(d+1).
*/
__attribute__((always_inline)) INLINE static void
hydro_runner_iact_nonsym_density_extra_kernel(struct part *restrict pi,
@@ -79,7 +98,11 @@ hydro_runner_iact_nonsym_density_extra_kernel(struct part *restrict pi,
const float r = sqrtf(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
const float r_inv = r ? 1.0f / r : 0.0f;
- /* Geometric moments and gradients (unmodified kernel, for normalisation) */
+ /* Geometric moments and gradients that use an unmodified kernel (Sandnes+2025
+ * Eqn. 50 and its gradient). Used in the normalising term (Eqn. 51) and in
+ * gradient estimates (using Eqn. 30) that are used for the calculation of grad-h
+ * terms (Eqn. 31) and in the artificial viscosity (Eqn. 35) and diffusion
+ * (Eqns. 46 and 47) schemes */
pi->m0 += pj->mass * wi / pj->rho_evol;
for (int i = 0; i < 3; i++) {
pi->grad_m0[i] += (pj->mass / pj->rho_evol) * dx[i] * wi_dx * r_inv;
@@ -99,7 +122,11 @@ hydro_end_density_extra_kernel(struct part *restrict p) {
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
- /* Geometric moments and gradients (unmodified kernel, for normalisation) */
+ /* Geometric moments and gradients that use an unmodified kernel (Sandnes+2025
+ * Eqn. 50 and its gradient). Used in the normalising term (Eqn. 51) and in
+ * gradient estimates (using Eqn. 30) that are used for the calculation of grad-h
+ * terms (Eqn. 31) and in the artificial viscosity (Eqn. 35) and diffusion
+ * (Eqns. 46 and 47) schemes */
p->m0 += p->mass * kernel_root / p->rho_evol;
p->m0 *= h_inv_dim;
for (int i = 0; i < 3; i++) {
@@ -124,7 +151,8 @@ hydro_prepare_gradient_extra_kernel(struct part *restrict p) {
zero_sym_matrix(&p->gradient.grad_m2_bar[2]);
zero_sym_matrix(&p->gradient.grad_m2_bar_gradhterm);
- /* Geometric moments and gradients (ij-mean, for linear-order repr. kernel) */
+ /* Geometric moments and gradients that us a kernel given by 0.5 * (W_{ij} + W_{ji}).
+ * These are used to construct the linear-order repr. kernel */
p->gradient.m0_bar = 0.f;
p->gradient.grad_m0_bar_gradhterm = 0.f;
for (int i = 0; i < 3; i++) {
@@ -139,6 +167,14 @@ hydro_prepare_gradient_extra_kernel(struct part *restrict p) {
/**
* @brief Extra kernel gradient interaction between two particles
+ *
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param wi The value of the unmodified kernel function W(r, hi) * hi^d.
+ * @param wj The value of the unmodified kernel function W(r, hj) * hj^d.
+ * @param wi_dx The norm of the gradient of wi: dW(r, hi)/dr * hi^(d+1).
+ * @param wj_dx The norm of the gradient of wj: dW(r, hj)/dr * hj^(d+1).
*/
__attribute__((always_inline)) INLINE static void
hydro_runner_iact_gradient_extra_kernel(struct part *restrict pi,
@@ -147,6 +183,7 @@ hydro_runner_iact_gradient_extra_kernel(struct part *restrict pi,
const float wj, const float wi_dx,
const float wj_dx) {
+ /* Get r and 1/r. */
const float r = sqrtf(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
const float r_inv = r ? 1.0f / r : 0.0f;
@@ -160,11 +197,11 @@ hydro_runner_iact_gradient_extra_kernel(struct part *restrict pi,
const float hj_inv_dim = pow_dimension(hj_inv); /* 1/h^d */
const float hj_inv_dim_plus_one = hj_inv_dim * hj_inv; /* 1/h^(d+1) */
- // Volume elements
+ /* Volume elements */
float volume_i = pi->mass / pi->rho_evol;
float volume_j = pj->mass / pj->rho_evol;
- // Mean ij kernels and gradients
+ /* Mean ij kernels and gradients */
float wi_term = 0.5f * (wi * hi_inv_dim + wj * hj_inv_dim);
float wj_term = wi_term;
float wi_dx_term[3], wj_dx_term[3];
@@ -175,13 +212,16 @@ hydro_runner_iact_gradient_extra_kernel(struct part *restrict pi,
(wi_dx * hi_inv_dim_plus_one + wj_dx * hj_inv_dim_plus_one);
}
- // Grad-h term, dW/dh
+ /* Grad-h term, dW/dh */
float wi_dx_gradhterm =
-0.5f * (hydro_dimension * wi + (r / hi) * wi_dx) * hi_inv_dim_plus_one;
float wj_dx_gradhterm =
-0.5f * (hydro_dimension * wj + (r / hj) * wj_dx) * hj_inv_dim_plus_one;
- // Geometric moments m_0, m_1, and m_2, their gradients and grad-h terms
+ /* Geometric moments m_0, m_1, and m_2 (Sandnes+2025 Eqns. 24--26), their
+ * gradients (Sandnes+2025 Eqns. B.10--B.12, initially we only construct the
+ * first terms in Eqns. B.11 and B.12) and grad-h terms (from second term
+ * in Eqn. 29 when used in Eqns. B.10--B.12)*/
pi->gradient.m0_bar += wi_term * volume_j;
pj->gradient.m0_bar += wj_term * volume_i;
@@ -229,12 +269,21 @@ hydro_runner_iact_gradient_extra_kernel(struct part *restrict pi,
/**
* @brief Extra kernel gradient interaction between two particles
* (non-symmetric)
+ *
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param wi The value of the unmodified kernel function W(r, hi) * hi^d.
+ * @param wj The value of the unmodified kernel function W(r, hj) * hj^d.
+ * @param wi_dx The norm of the gradient of wi: dW(r, hi)/dr * hi^(d+1).
+ * @param wj_dx The norm of the gradient of wj: dW(r, hj)/dr * hj^(d+1).
*/
__attribute__((always_inline)) INLINE static void
hydro_runner_iact_nonsym_gradient_extra_kernel(
struct part *restrict pi, struct part *restrict pj, const float dx[3],
const float wi, const float wj, const float wi_dx, const float wj_dx) {
+ /* Get r and 1/r. */
const float r = sqrtf(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
const float r_inv = r ? 1.0f / r : 0.0f;
@@ -248,10 +297,10 @@ hydro_runner_iact_nonsym_gradient_extra_kernel(
const float hj_inv_dim = pow_dimension(hj_inv); /* 1/h^d */
const float hj_inv_dim_plus_one = hj_inv_dim * hj_inv; /* 1/h^(d+1) */
- // Volume elements
+ /* Volume elements */
float volume_j = pj->mass / pj->rho_evol;
- // Mean ij kernel and gradients
+ /* Mean ij kernel and gradients */
float wi_term = 0.5f * (wi * hi_inv_dim + wj * hj_inv_dim);
float wi_dx_term[3];
for (int i = 0; i < 3; i++) {
@@ -259,11 +308,14 @@ hydro_runner_iact_nonsym_gradient_extra_kernel(
(wi_dx * hi_inv_dim_plus_one + wj_dx * hj_inv_dim_plus_one);
}
- // Grad-h term, dW/dh
+ /* Grad-h term, dW/dh */
float wi_dx_gradhterm =
-0.5f * (hydro_dimension * wi + (r / hi) * wi_dx) * hi_inv_dim_plus_one;
- // Geometric moments m_0, m_1, and m_2, their gradients and grad-h terms
+ /* Geometric moments m_0, m_1, and m_2 (Sandnes+2025 Eqns. 24--26), their
+ * gradients (Sandnes+2025 Eqns. B.10--B.12, initially we only construct the
+ * first terms in Eqns. B.11 and B.12) and grad-h terms (from second term
+ * in Eqn. 29 when used in Eqns. B.10--B.12)*/
pi->gradient.m0_bar += wi_term * volume_j;
pi->gradient.grad_m0_bar_gradhterm += wi_dx_gradhterm * volume_j;
@@ -309,15 +361,16 @@ hydro_end_gradient_extra_kernel(struct part *restrict p) {
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
- // Volume elements
+ /* Volume elements */
float volume = p->mass / p->rho_evol;
- // Self contribution to geometric moments and gradients
+ /* Self contribution to geometric moments and gradients */
p->gradient.m0_bar += volume * kernel_root * h_inv_dim;
p->gradient.grad_m0_bar_gradhterm -=
0.5f * volume * hydro_dimension * kernel_root * h_inv_dim_plus_one;
- // Multiply dh/dr into grad-h terms
+ /* Multiply dh/dr (Sandnes+2025 Eqn. 31) into grad-h terms (See second term
+ * in Sandnes+2025 Eqn. 29) */
for (int i = 0; i < 3; i++) {
p->gradient.grad_m0_bar[i] +=
p->gradient.grad_m0_bar_gradhterm * p->dh_norm_kernel[i];
@@ -342,7 +395,8 @@ hydro_prepare_force_extra_kernel(struct part *restrict p) {
if (!p->is_h_max) {
- // Add second terms to complete the geometric moment gradients
+ /* Add second terms in Sandnes+2025 Eqns. B.11 and B.12 to complete the
+ * geometric moment gradients */
for (int i = 0; i < 3; i++) {
p->gradient.grad_m1_bar[i][i] += p->gradient.m0_bar;
}
@@ -359,16 +413,19 @@ hydro_prepare_force_extra_kernel(struct part *restrict p) {
p->gradient.grad_m2_bar[2].xz += p->gradient.m1_bar[0];
p->gradient.grad_m2_bar[2].yz += p->gradient.m1_bar[1];
- // Inverse of symmetric geometric moment m_2 (bar)
+ /* Inverse of symmetric geometric moment m_2 (bar) */
struct sym_matrix m2_bar_inv;
- // Make m2_bar dimensionless for calculation of inverse
+ /* Make m2_bar dimensionless for calculation of inverse */
struct sym_matrix m2_bar_over_h2;
for (int i = 0; i < 6; i++) m2_bar_over_h2.elements[i] =
p->gradient.m2_bar.elements[i] / (p->h * p->h);
sym_matrix_invert(&m2_bar_inv, &m2_bar_over_h2);
for (int i = 0; i < 6; i++) m2_bar_inv.elements[i] /= (p->h * p->h);
- // Components for constructing linear-order kernel's A and B, and gradients
+ /* Components for constructing linear-order kernel's A and B (Sandnes+2025
+ * Eqns. 22 and 23), and gradients (Eqns. B.8 and B.9) that are calculated
+ * with sym_matrix functions from combinations of geometric moments and
+ * their gradients */
float m2_bar_inv_mult_m1_bar[3];
float m2_bar_inv_mult_grad_m1_bar[3][3];
struct sym_matrix m2_bar_inv_mult_grad_m2_bar_mult_m2_bar_inv[3];
@@ -387,22 +444,23 @@ hydro_prepare_force_extra_kernel(struct part *restrict p) {
p->gradient.m1_bar);
}
- // Linear-order reproducing kernel's A and B components and gradients
+ /* Linear-order reproducing kernel's A and B components (Sandnes+2025
+ * Eqns. 22 and 23) and gradients (Eqns. B.8 and B.9) */
float A, B[3], grad_A[3], grad_B[3][3];
- // Calculate A
+ /* Calculate A (Sandnes+2025 Eqn. 22) */
A = p->gradient.m0_bar;
for (int i = 0; i < 3; i++) {
A -= m2_bar_inv_mult_m1_bar[i] * p->gradient.m1_bar[i];
}
A = 1.f / A;
- // Calculate B
+ /* Calculate B (Sandnes+2025 Eqn. 23) */
for (int i = 0; i < 3; i++) {
B[i] = -m2_bar_inv_mult_m1_bar[i];
}
- // Calculate grad A
+ /* Calculate grad A (Sandnes+2025 Eqn. B.8) */
for (int i = 0; i < 3; i++) {
grad_A[i] = p->gradient.grad_m0_bar[i];
@@ -415,7 +473,7 @@ hydro_prepare_force_extra_kernel(struct part *restrict p) {
grad_A[i] *= -A * A;
}
- // Calculate grad B
+ /* Calculate grad B (Sandnes+2025 Eqn. B.9) */
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
grad_B[j][i] =
@@ -423,7 +481,7 @@ hydro_prepare_force_extra_kernel(struct part *restrict p) {
}
}
- // Store final values
+ /* Store final values */
p->force.A = A;
for (int i = 0; i < 3; i++) {
p->force.B[i] = B[i];
@@ -433,7 +491,7 @@ hydro_prepare_force_extra_kernel(struct part *restrict p) {
}
}
- // Vacuum-boundary proximity switch
+ /* Vacuum-boundary proximity switch (Sandnes+2025 Eqn. 33) */
p->force.vac_switch = 1.f;
float hB = p->h * sqrtf(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
float offset = 0.8f;
@@ -445,7 +503,6 @@ hydro_prepare_force_extra_kernel(struct part *restrict p) {
}
} else {
- // p->is_h_max
p->force.A = 1.f;
p->force.vac_switch = 1.f;
for (int i = 0; i < 3; i++) {
@@ -461,14 +518,25 @@ hydro_prepare_force_extra_kernel(struct part *restrict p) {
/**
* @brief Set gradient terms for linear-order reproducing kernel.
*
- * Note `G` here corresponds to `d/dr tilde{mathcal{W}}` in Sandnes et al.
- * (2024)
+ * Note `G` here corresponds to `d/dr tilde{mathcal{W}}` in Sandnes et al.(2025).
+ * These are used in the REMIX equations of motion (Sandnes+2025 Eqns. 14--16).
+ *
+ * @param Gi (return) Gradient of linear-order reproducing kernel for first particle.
+ * @param Gj (return) Gradient of linear-order reproducing kernel for second particle.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param wi The value of the unmodified kernel function W(r, hi) * hi^d.
+ * @param wj The value of the unmodified kernel function W(r, hj) * hj^d.
+ * @param wi_dx The norm of the gradient of wi: dW(r, hi)/dr * hi^(d+1).
+ * @param wj_dx The norm of the gradient of wj: dW(r, hj)/dr * hj^(d+1).
*/
__attribute__((always_inline)) INLINE static void hydro_set_Gi_Gj_forceloop(
float Gi[3], float Gj[3], const struct part *restrict pi,
const struct part *restrict pj, const float dx[3], const float wi,
const float wj, const float wi_dx, const float wj_dx) {
+ /* Get r and 1/r. */
const float r = sqrtf(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
const float r_inv = r ? 1.0f / r : 0.0f;
@@ -486,10 +554,13 @@ __attribute__((always_inline)) INLINE static void hydro_set_Gi_Gj_forceloop(
const float wj_dr = hj_inv_dim_plus_one * wj_dx;
if ((!pi->is_h_max) && (!pj->is_h_max)) {
- // Mean ij kernels and gradients with grad-h terms
+ /* Mean ij kernels and gradients with grad-h terms */
float wi_term = 0.5f * (wi * hi_inv_dim + wj * hj_inv_dim);
float wj_term = wi_term;
float wi_dx_term[3], wj_dx_term[3];
+
+ /* Get linear-order reproducing kernel's A and B components (Sandnes+2025
+ * Eqns. 22 and 23) and gradients (Eqns. B.8 and B.9) */
float Ai = pi->force.A;
float Aj = pj->force.A;
float grad_Ai[3], grad_Aj[3], Bi[3], Bj[3], grad_Bi[3][3], grad_Bj[3][3];
@@ -505,6 +576,8 @@ __attribute__((always_inline)) INLINE static void hydro_set_Gi_Gj_forceloop(
}
for (int i = 0; i < 3; i++) {
+
+ /* Assemble Sandnes+2025 Eqn. 29 */
wi_dx_term[i] =
dx[i] * r_inv * 0.5f *
(wi_dx * hi_inv_dim_plus_one + wj_dx * hj_inv_dim_plus_one);
@@ -517,6 +590,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_Gi_Gj_forceloop(
wj_dx_term[i] += -0.5f * (hydro_dimension * wj + (r / pj->h) * wj_dx) *
hj_inv_dim_plus_one * pj->dh_norm_kernel[i];
+ /* Assemble Sandnes+2025 Eqn. 28 */
Gi[i] = Ai * wi_dx_term[i] + grad_Ai[i] * wi_term + Ai * Bi[i] * wi_term;
Gj[i] = Aj * wj_dx_term[i] + grad_Aj[i] * wj_term + Aj * Bj[i] * wj_term;
@@ -531,7 +605,8 @@ __attribute__((always_inline)) INLINE static void hydro_set_Gi_Gj_forceloop(
}
}
- // Standard-kernel gradients, to be used for vacuum boundary switch
+ /* Standard-kernel gradients, to be used for vacuum boundary switch (For
+ * second term in Sandnes+2025 Eqn. 32)*/
float wi_dx_term_vac[3], wj_dx_term_vac[3];
for (int i = 0; i < 3; i++) {
wi_dx_term_vac[i] = dx[i] * r_inv * wi_dx * hi_inv_dim_plus_one;
@@ -543,7 +618,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_Gi_Gj_forceloop(
hj_inv_dim_plus_one * pj->dh_norm_kernel[i];
}
- // Gradients, including vacuum boundary switch
+ /* Gradients, including vacuum boundary switch (Sandnes+2025 Eqn. 32) */
for (int i = 0; i < 3; i++) {
Gi[i] = pi->force.vac_switch * Gi[i] +
(1.f - pi->force.vac_switch) * wi_dx_term_vac[i];
@@ -552,7 +627,8 @@ __attribute__((always_inline)) INLINE static void hydro_set_Gi_Gj_forceloop(
}
} else {
- // One or both particles at h_max
+ /* If one or both particles at h_max, revert to standard kernel grads, without
+ * grad-h terms */
for (int i = 0; i < 3; i++) {
Gi[i] = wi_dr * dx[i] * r_inv;
Gj[i] = -wj_dr * dx[i] * r_inv;
diff --git a/src/hydro/REMIX/hydro_parameters.h b/src/hydro/REMIX/hydro_parameters.h
index da4c955b476f57c3aa49dbcf2c3af06cb1813ecd..ed2932814bff8fe5664a883af5c7bf708a3330fd 100644
--- a/src/hydro/REMIX/hydro_parameters.h
+++ b/src/hydro/REMIX/hydro_parameters.h
@@ -36,7 +36,7 @@
/**
* @file Planetary/hydro_parameters.h
- * @brief REMIX implementation of SPH (Sandnes et al. 2024) (default parameters)
+ * @brief REMIX implementation of SPH (Sandnes et al. 2025) (default parameters)
*
* This file defines a number of things that are used in
* hydro_properties.c as defaults for run-time parameters
diff --git a/src/hydro/REMIX/hydro_part.h b/src/hydro/REMIX/hydro_part.h
index 01eb2357c848c5b128d0274b391d5cde306d662c..4a13cc2140a2c7cdbf810dddf0218786d793226e 100644
--- a/src/hydro/REMIX/hydro_part.h
+++ b/src/hydro/REMIX/hydro_part.h
@@ -23,13 +23,13 @@
/**
* @file Planetary/hydro_part.h
- * @brief REMIX implementation of SPH (Sandnes et al. 2024)
+ * @brief REMIX implementation of SPH (Sandnes et al. 2025)
*/
#include "black_holes_struct.h"
#include "chemistry_struct.h"
#include "cooling_struct.h"
-#include "equation_of_state.h" // For enum material_id
+#include "equation_of_state.h" /* For enum material_id */
#include "feedback_struct.h"
#include "mhd_struct.h"
#include "particle_splitting_struct.h"
diff --git a/src/hydro/REMIX/hydro_visc_difn.h b/src/hydro/REMIX/hydro_visc_difn.h
index c5d95e9302b9141ea557fb2358c881c8e792bebc..985baa3c4a34a163a34912423e4691d407192920 100644
--- a/src/hydro/REMIX/hydro_visc_difn.h
+++ b/src/hydro/REMIX/hydro_visc_difn.h
@@ -30,50 +30,64 @@
#include "math.h"
/**
- * @brief Prepares extra viscosity parameters for a particle for the density
+ * @brief Prepares extra artificial viscosity and artificial diffusion parameters for a particle for the density
* calculation.
*
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void
-hydro_init_part_extra_viscosity(struct part *restrict p) {}
+hydro_init_part_extra_visc_difn(struct part *restrict p) {}
/**
- * @brief Extra density interaction between two particles
+ * @brief Extra artificial viscosity and artificial diffusion density interaction between two particles
+ *
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param wi The value of the unmodified kernel function W(r, hi) * hi^d.
+ * @param wj The value of the unmodified kernel function W(r, hj) * hj^d.
+ * @param wi_dx The norm of the gradient of wi: dW(r, hi)/dr * hi^(d+1).
+ * @param wj_dx The norm of the gradient of wj: dW(r, hj)/dr * hj^(d+1).
*/
__attribute__((always_inline)) INLINE static void
-hydro_runner_iact_density_extra_viscosity(struct part *restrict pi,
+hydro_runner_iact_density_extra_visc_difn(struct part *restrict pi,
struct part *restrict pj,
const float dx[3], const float wi,
const float wj, const float wi_dx,
const float wj_dx) {}
/**
- * @brief Extra density interaction between two particles (non-symmetric)
+ * @brief Extra artificial viscosity and artificial diffusion density interaction between two particles (non-symmetric)
+ *
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param wi The value of the unmodified kernel function W(r, hi) * hi^d.
+ * @param wi_dx The norm of the gradient of wi: dW(r, hi)/dr * hi^(d+1).
*/
__attribute__((always_inline)) INLINE static void
-hydro_runner_iact_nonsym_density_extra_viscosity(struct part *restrict pi,
+hydro_runner_iact_nonsym_density_extra_visc_difn(struct part *restrict pi,
const struct part *restrict pj,
const float dx[3],
const float wi,
const float wi_dx) {}
/**
- * @brief Finishes extra viscosity parts of the density calculation.
+ * @brief Finishes extra artificial viscosity and artificial diffusion parts of the density calculation.
*
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void
-hydro_end_density_extra_viscosity(struct part *restrict p) {}
+hydro_end_density_extra_visc_difn(struct part *restrict p) {}
/**
- * @brief Prepares extra viscosity parameters for a particle for the gradient
+ * @brief Prepares extra artificial viscosity and artificial diffusion parameters for a particle for the gradient
* calculation.
*
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void
-hydro_prepare_gradient_extra_viscosity(struct part *restrict p) {
+hydro_prepare_gradient_extra_visc_difn(struct part *restrict p) {
for (int i = 0; i < 3; i++) {
p->du_norm_kernel[i] = 0.f;
@@ -86,15 +100,24 @@ hydro_prepare_gradient_extra_viscosity(struct part *restrict p) {
}
/**
- * @brief Extra gradient interaction between two particles
+ * @brief Extra artificial viscosity and artificial diffusion gradient interaction between two particles
+ *
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param wi The value of the unmodified kernel function W(r, hi) * hi^d.
+ * @param wj The value of the unmodified kernel function W(r, hj) * hj^d.
+ * @param wi_dx The norm of the gradient of wi: dW(r, hi)/dr * hi^(d+1).
+ * @param wj_dx The norm of the gradient of wj: dW(r, hj)/dr * hj^(d+1).
*/
__attribute__((always_inline)) INLINE static void
-hydro_runner_iact_gradient_extra_viscosity(struct part *restrict pi,
+hydro_runner_iact_gradient_extra_visc_difn(struct part *restrict pi,
struct part *restrict pj,
const float dx[3], const float wi,
const float wj, const float wi_dx,
const float wj_dx) {
+ /* Get r and 1/r. */
const float r = sqrtf(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
const float r_inv = r ? 1.0f / r : 0.0f;
@@ -111,7 +134,7 @@ hydro_runner_iact_gradient_extra_viscosity(struct part *restrict pi,
float volume_i = pi->mass / pi->rho_evol;
float volume_j = pj->mass / pj->rho_evol;
- // Gradients of normalised kernel
+ /* Gradients of normalised kernel (Sandnes+2025 Eqn. 30) */
float wi_dx_term[3], wj_dx_term[3];
for (int i = 0; i < 3; i++) {
wi_dx_term[i] = dx[i] * r_inv * wi_dx * hi_inv_dim_plus_one;
@@ -124,12 +147,13 @@ hydro_runner_iact_gradient_extra_viscosity(struct part *restrict pi,
wj_dx_term[i] += -wj * hj_inv_dim * pj->grad_m0[i] / (pj->m0 * pj->m0);
}
- // Normalised kernel h, u, rho, v gradients
+ /* Gradient estimates of h (Sandnes+2025 Eqn. 31), v (Eqn. 35), u (Eqn. 46),
+ * rho (Eqn. 47) using normalised kernel */
for (int i = 0; i < 3; i++) {
pi->dh_norm_kernel[i] += (pj->h - pi->h) * wi_dx_term[i] * volume_j;
pj->dh_norm_kernel[i] += (pi->h - pj->h) * wj_dx_term[i] * volume_i;
- // Contributions only from same-material particles for u and rho diffusion
+ /* Contributions only from same-material particles for u and rho diffusion */
if (pi->mat_id == pj->mat_id) {
pi->du_norm_kernel[i] += (pj->u - pi->u) * wi_dx_term[i] * volume_j;
pj->du_norm_kernel[i] += (pi->u - pj->u) * wj_dx_term[i] * volume_i;
@@ -150,13 +174,20 @@ hydro_runner_iact_gradient_extra_viscosity(struct part *restrict pi,
}
/**
- * @brief Extra gradient interaction between two particles (non-symmetric)
+ * @brief Extra artificial viscosity and artificial diffusion gradient interaction between two particles (non-symmetric)
+ *
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param wi The value of the unmodified kernel function W(r, hi) * hi^d.
+ * @param wi_dx The norm of the gradient of wi: dW(r, hi)/dr * hi^(d+1).
*/
__attribute__((always_inline)) INLINE static void
-hydro_runner_iact_nonsym_gradient_extra_viscosity(
+hydro_runner_iact_nonsym_gradient_extra_visc_difn(
struct part *restrict pi, const struct part *restrict pj, const float dx[3],
const float wi, const float wi_dx) {
+ /* Get r and 1/r. */
const float r = sqrtf(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
const float r_inv = r ? 1.0f / r : 0.0f;
@@ -167,7 +198,7 @@ hydro_runner_iact_nonsym_gradient_extra_viscosity(
float volume_j = pj->mass / pj->rho_evol;
- // Gradients of normalised kernel
+ /* Gradients of normalised kernel (Sandnes+2025 Eqn. 30) */
float wi_dx_term[3];
for (int i = 0; i < 3; i++) {
wi_dx_term[i] = dx[i] * r_inv * wi_dx * hi_inv_dim_plus_one;
@@ -175,11 +206,12 @@ hydro_runner_iact_nonsym_gradient_extra_viscosity(
wi_dx_term[i] += -wi * hi_inv_dim * pi->grad_m0[i] / (pi->m0 * pi->m0);
}
- // Normalised kernel h, u, rho, v gradients
+ /* Gradient estimates of h (Sandnes+2025 Eqn. 31), v (Eqn. 35), u (Eqn. 46),
+ * rho (Eqn. 47) using normalised kernel */
for (int i = 0; i < 3; i++) {
pi->dh_norm_kernel[i] += (pj->h - pi->h) * wi_dx_term[i] * volume_j;
- // Contributions only from same-material particles for u and rho diffusion
+ /* Contributions only from same-material particles for u and rho diffusion */
if (pi->mat_id == pj->mat_id) {
pi->du_norm_kernel[i] += (pj->u - pi->u) * wi_dx_term[i] * volume_j;
pi->drho_norm_kernel[i] +=
@@ -194,21 +226,28 @@ hydro_runner_iact_nonsym_gradient_extra_viscosity(
}
/**
- * @brief Finishes extra viscosity parts of the gradient calculation.
+ * @brief Finishes extra artificial viscosity and artificial diffusion parts of the gradient calculation.
*
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void
-hydro_end_gradient_extra_viscosity(struct part *restrict p) {}
+hydro_end_gradient_extra_visc_difn(struct part *restrict p) {}
/**
* @brief Returns particle viscous pressures
+ *
+ * @param Qi (return) Viscous pressure for first particle.
+ * @param Qj (return) Viscous pressure for second particle.
+ * @param visc_signal_velocity (return) Signal velocity of artificial viscosity.
+ * @param difn_signal_velocity (return) Signal velocity of artificial diffusion.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param dx Comoving vector separating both particles (pi - pj).
*/
__attribute__((always_inline)) INLINE static void hydro_set_Qi_Qj(
float *Qi, float *Qj, float *visc_signal_velocity,
float *difn_signal_velocity, const struct part *restrict pi,
- const struct part *restrict pj, const float dx[3], const float a,
- const float H) {
+ const struct part *restrict pj, const float dx[3]) {
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
const float r = sqrtf(r2);
@@ -216,10 +255,11 @@ __attribute__((always_inline)) INLINE static void hydro_set_Qi_Qj(
const float hi_inv = 1.0f / pi->h;
const float hj_inv = 1.0f / pj->h;
- // Reconstructed velocities at the halfway point between particles
+ /* Reconstructed velocities at the halfway point between particles (Sandnes+2025
+ * Eqn. 36) */
float vtilde_i[3], vtilde_j[3];
- // Viscosity parameters
+ /* Viscosity parameters. These are set in hydro_parameters.h */
const float alpha = const_remix_visc_alpha;
const float beta = const_remix_visc_beta;
const float epsilon = const_remix_visc_epsilon;
@@ -229,29 +269,27 @@ __attribute__((always_inline)) INLINE static void hydro_set_Qi_Qj(
const float slope_limiter_exp_denom = const_remix_slope_limiter_exp_denom;
if ((!pi->is_h_max) && (!pj->is_h_max)) {
- /* A numerators and denominators (eq 22 in Rosswog 2020) */
+ /* A numerators and denominators (Sandnes+2025 Eqn. 38) */
float A_i_v = 0.f;
float A_j_v = 0.f;
- /* Terms in square brackets in Rosswog 2020 eq 17 */
+ /* 1/2 * (r_j - r_i) * dv/dr in second term of Sandnes+2025 Eqn. 38 */
float v_reconst_i[3] = {0};
float v_reconst_j[3] = {0};
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
- /* Get the A numerators and denominators (eq 22 in Rosswog 2020). dv
- * is from eq 18 */
+ /* Get the A numerators and denominators (Sandnes+2025 Eqn. 38). dv_norm_kernel
+ * is from Eqn. 35 */
A_i_v += pi->dv_norm_kernel[i][j] * dx[i] * dx[j];
A_j_v += pj->dv_norm_kernel[i][j] * dx[i] * dx[j];
- /* Terms in square brackets in Rosswog 2020 eq 17. Add in FIRST
- * derivative terms */
v_reconst_i[j] -= 0.5 * pi->dv_norm_kernel[i][j] * dx[i];
v_reconst_j[j] += 0.5 * pj->dv_norm_kernel[i][j] * dx[i];
}
}
- // ...
+ /* Slope limiter (Sandnes+2025 Eqn. 37) special cases */
float phi_i_v, phi_j_v;
if ((A_i_v == 0.f) && (A_j_v == 0.f)) {
phi_i_v = 1.f;
@@ -262,7 +300,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_Qi_Qj(
phi_i_v = 0.f;
phi_j_v = 0.f;
} else {
- /* Slope limiter (eq 21 in Rosswog 2020) */
+ /* Slope limiter (Sandnes+2025 Eqn. 37) */
phi_i_v = min(1.f, 4.f * A_i_v / A_j_v / (1.f + A_i_v / A_j_v) /
(1.f + A_i_v / A_j_v));
phi_i_v = max(0.f, phi_i_v);
@@ -272,7 +310,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_Qi_Qj(
phi_j_v = max(0.f, phi_j_v);
}
- // ...
+ /* exp in slope limiter (middle case in Sandnes+2025 Eqn. 37) */
float eta_ab = min(r * hi_inv, r * hj_inv);
if (eta_ab < eta_crit) {
phi_i_v *= expf(-(eta_ab - eta_crit) * (eta_ab - eta_crit) /
@@ -281,9 +319,8 @@ __attribute__((always_inline)) INLINE static void hydro_set_Qi_Qj(
slope_limiter_exp_denom);
}
- // ...
for (int i = 0; i < 3; i++) {
- /* Assemble the reconstructed velocity (eq 17 in Rosswog 2020) */
+ /* Assemble the reconstructed velocity (Sandnes+2025 Eqn. 36) */
vtilde_i[i] =
pi->v[i] + (1.f - pi->force.balsara) * phi_i_v * v_reconst_i[i];
vtilde_j[i] =
@@ -297,7 +334,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_Qi_Qj(
}
}
- /* Finally assemble eq 15 in Rosswog 2020 */
+ /* Assemble Sandnes+2025 Eqn. 40 */
float mu_i =
min(0.f, ((vtilde_i[0] - vtilde_j[0]) * dx[0] +
(vtilde_i[1] - vtilde_j[1]) * dx[1] +
@@ -312,27 +349,38 @@ __attribute__((always_inline)) INLINE static void hydro_set_Qi_Qj(
const float ci = pi->force.soundspeed;
const float cj = pj->force.soundspeed;
- /* Get viscous pressure terms (eq 14 in Rosswog 2020) */
+ /* Finally assemble viscous pressure terms (Sandnes+2025 41) */
*Qi = (a_visc + b_visc * pi->force.balsara) * 0.5f * pi->rho *
(-alpha * ci * mu_i + beta * mu_i * mu_i);
*Qj = (a_visc + b_visc * pj->force.balsara) * 0.5f * pj->rho *
(-alpha * cj * mu_j + beta * mu_j * mu_j);
- // Account for alpha being outside brackets in timestep code
+ /* Account for alpha being outside brackets in timestep code */
const float viscosity_parameter_factor = (alpha == 0.f) ? 0.f : beta / alpha;
*visc_signal_velocity = ci + cj - 2.f * viscosity_parameter_factor * min(mu_i, mu_j);
- // ...
+ /* Signal velocity used for the artificial diffusion (Sandnes+2025 Eqns. 42 and 43) */
*difn_signal_velocity =
sqrtf((vtilde_i[0] - vtilde_j[0]) * (vtilde_i[0] - vtilde_j[0]) +
(vtilde_i[1] - vtilde_j[1]) * (vtilde_i[1] - vtilde_j[1]) +
(vtilde_i[2] - vtilde_j[2]) * (vtilde_i[2] - vtilde_j[2]));
}
+/**
+ * @brief Returns midpoint reconstructions of internal energies and densities
+ *
+ * @param utilde_i (return) u reconstructed to midpoint from first particle.
+ * @param utilde_j (return) u reconstructed to midpoint from second particle.
+ * @param rhotilde_i (return) rho reconstructed to midpoint from first particle.
+ * @param rhotilde_j (return) rho reconstructed to midpoint from second particle.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ */
__attribute__((always_inline)) INLINE static void hydro_set_u_rho_difn(
float *utilde_i, float *utilde_j, float *rhotilde_i, float *rhotilde_j,
const struct part *restrict pi, const struct part *restrict pj,
- const float dx[3], const float a, const float H) {
+ const float dx[3]) {
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
const float r = sqrtf(r2);
@@ -344,28 +392,28 @@ __attribute__((always_inline)) INLINE static void hydro_set_u_rho_difn(
const float slope_limiter_exp_denom = const_remix_slope_limiter_exp_denom;
if ((!pi->is_h_max) && (!pj->is_h_max)) {
- /* A numerators and denominators (eq 22 in Rosswog 2020) */
+ /* A numerators and denominators (Sandnes+2025 Eqns. 48 and 49) */
float A_i_u = 0.f;
float A_j_u = 0.f;
float A_i_rho = 0.f;
float A_j_rho = 0.f;
- /* Terms in square brackets in Rosswog 2020 eq 17 */
+ /* 1/2 * (r_j - r_i) * du/dr in second term of Sandnes+2025 Eqn. 44 */
float u_reconst_i = 0.f;
float u_reconst_j = 0.f;
+
+ /* 1/2 * (r_j - r_i) * drho/dr in second term of Sandnes+2025 Eqn. 45 */
float rho_reconst_i = 0.f;
float rho_reconst_j = 0.f;
for (int i = 0; i < 3; i++) {
- /* Get the A numerators and denominators (eq 22 in Rosswog 2020). dv
- * is from eq 18 */
+ /* Get the A numerators and denominators (Sandnes+2025 Eqns. 48 and 49).
+ * du_norm_kernel is from Eqn. 46 and drho_norm_kernel is from Eqn. 47 */
A_i_u += pi->du_norm_kernel[i] * dx[i];
A_j_u += pj->du_norm_kernel[i] * dx[i];
A_i_rho += pi->drho_norm_kernel[i] * dx[i];
A_j_rho += pj->drho_norm_kernel[i] * dx[i];
- /* Terms in square brackets in Rosswog 2020 eq 17. Add in FIRST
- * derivative terms */
u_reconst_i -= 0.5 * pi->du_norm_kernel[i] * dx[i];
u_reconst_j += 0.5 * pj->du_norm_kernel[i] * dx[i];
rho_reconst_i -= 0.5 * pi->drho_norm_kernel[i] * dx[i];
@@ -373,7 +421,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_u_rho_difn(
}
float phi_i_u, phi_j_u, phi_i_rho, phi_j_rho;
-
+ /* Slope limiter (Sandnes+2025 Eqn. 37) special cases */
if ((A_i_u == 0.f) && (A_j_u == 0.f)) {
phi_i_u = 1.f;
phi_j_u = 1.f;
@@ -383,7 +431,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_u_rho_difn(
phi_i_u = 0.f;
phi_j_u = 0.f;
} else {
- /* Slope limiter (eq 21 in Rosswog 2020) */
+ /* Slope limiter (Sandnes+2025 Eqn. 37) */
phi_i_u = min(1.f, 4.f * A_i_u / A_j_u / (1.f + A_i_u / A_j_u) /
(1.f + A_i_u / A_j_u));
phi_i_u = max(0.f, phi_i_u);
@@ -393,6 +441,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_u_rho_difn(
phi_j_u = max(0.f, phi_j_u);
}
+ /* Slope limiter (Sandnes+2025 Eqn. 37) special cases */
if ((A_i_rho == 0.f) && (A_j_rho == 0.f)) {
phi_i_rho = 1.f;
phi_j_rho = 1.f;
@@ -402,7 +451,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_u_rho_difn(
phi_i_rho = 0.f;
phi_j_rho = 0.f;
} else {
- /* Slope limiter (eq 21 in Rosswog 2020) */
+ /* Slope limiter (Sandnes+2025 Eqn. 37) */
phi_i_rho = min(1.f, 4.f * A_i_rho / A_j_rho / (1.f + A_i_rho / A_j_rho) /
(1.f + A_i_rho / A_j_rho));
phi_i_rho = max(0.f, phi_i_rho);
@@ -412,9 +461,8 @@ __attribute__((always_inline)) INLINE static void hydro_set_u_rho_difn(
phi_j_rho = max(0.f, phi_j_rho);
}
- // ...
+ /* exp in slope limiter (middle case in Sandnes+2025 Eqn. 37) */
float eta_ab = min(r * hi_inv, r * hj_inv);
-
if (eta_ab < eta_crit) {
phi_i_u *= expf(-(eta_ab - eta_crit) * (eta_ab - eta_crit) /
slope_limiter_exp_denom);
@@ -426,14 +474,15 @@ __attribute__((always_inline)) INLINE static void hydro_set_u_rho_difn(
slope_limiter_exp_denom);
}
- /* Assemble the reconstructed velocity (eq 17 in Rosswog 2020) */
+ /* Assemble the reconstructed internal energy (Sandnes+2025 Eqn. 44) and
+ * density (Sandnes+2025 Eqn. 45) */
*utilde_i = pi->u + phi_i_u * u_reconst_i;
*utilde_j = pj->u + phi_j_u * u_reconst_j;
*rhotilde_i = pi->rho_evol + phi_i_rho * rho_reconst_i;
*rhotilde_j = pj->rho_evol + phi_j_rho * rho_reconst_j;
} else {
- /* If h=h_max don't reconstruct velocity */
+ /* If h=h_max don't reconstruct internal energy of density */
*utilde_i = 0.f;
*utilde_j = 0.f;
*rhotilde_i = 0.f;