From ea726709aced5427fccfc2437232432ecea9a331 Mon Sep 17 00:00:00 2001
From: Thomas Sandnes <thomas.d.sandnes@durham.ac.uk>
Date: Wed, 9 Apr 2025 09:05:44 +0000
Subject: [PATCH] REMIX hydro scheme
---
configure.ac | 7 +-
.../kelvin_helmholtz.yml | 2 +-
src/debug.c | 2 +
src/hydro.h | 8 +-
src/hydro/REMIX/hydro.h | 981 ++++++++++++++++++
src/hydro/REMIX/hydro_debug.h | 53 +
src/hydro/REMIX/hydro_iact.h | 539 ++++++++++
src/hydro/REMIX/hydro_io.h | 251 +++++
src/hydro/REMIX/hydro_kernels.h | 644 ++++++++++++
src/hydro/REMIX/hydro_parameters.h | 186 ++++
src/hydro/REMIX/hydro_part.h | 313 ++++++
src/hydro/REMIX/hydro_visc_difn.h | 452 ++++++++
src/hydro_csds.h | 2 +
src/hydro_io.h | 2 +
src/hydro_parameters.h | 2 +
src/part.h | 5 +
src/runner_main.c | 18 +
src/symmetric_matrix.h | 167 +++
src/tools.c | 269 +++--
tests/makeInput.py | 17 +-
tests/test125cells.c | 38 +-
tests/test27cells.c | 3 +
tests/testActivePair.c | 47 +
tests/testEOS.c | 291 +-----
tests/testEOS.sh | 27 -
tests/testInteractions.c | 3 +-
tests/testPeriodicBC.c | 3 +
tests/testRiemannTRRS.c | 6 +
28 files changed, 3920 insertions(+), 418 deletions(-)
create mode 100644 src/hydro/REMIX/hydro.h
create mode 100644 src/hydro/REMIX/hydro_debug.h
create mode 100644 src/hydro/REMIX/hydro_iact.h
create mode 100644 src/hydro/REMIX/hydro_io.h
create mode 100644 src/hydro/REMIX/hydro_kernels.h
create mode 100644 src/hydro/REMIX/hydro_parameters.h
create mode 100644 src/hydro/REMIX/hydro_part.h
create mode 100644 src/hydro/REMIX/hydro_visc_difn.h
create mode 100644 src/symmetric_matrix.h
delete mode 100755 tests/testEOS.sh
diff --git a/configure.ac b/configure.ac
index 3db764605d..4e920c2b8d 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1449,7 +1449,7 @@ if test "x$with_lustreapi" != "xno"; then
LUSTREAPI_LIBS="-llustreapi"
LUSTREAPI_INCS=""
fi
- AC_CHECK_LIB([lustreapi], [llapi_obd_statfs], [have_lustreapi="yes"],
+ AC_CHECK_LIB([lustreapi], [llapi_obd_statfs], [have_lustreapi="yes"],
[have_lustreapi="no"],[$LUSTREAPI_LIBS])
if test "$have_lustreapi" = "yes"; then
@@ -2223,7 +2223,7 @@ fi
# Hydro scheme.
AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>],
- [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, phantom, gizmo-mfv, gizmo-mfm, shadowswift, planetary, sphenix, gasoline, anarchy-pu default: sphenix@:>@]
+ [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, phantom, gizmo-mfv, gizmo-mfm, shadowswift, planetary, remix, sphenix, gasoline, anarchy-pu default: sphenix@:>@]
)],
[with_hydro="$withval"],
[with_hydro="sphenix"]
@@ -2270,6 +2270,9 @@ case "$with_hydro" in
planetary)
AC_DEFINE([PLANETARY_SPH], [1], [Planetary SPH])
;;
+ remix)
+ AC_DEFINE([REMIX_SPH], [1], [REMIX SPH])
+ ;;
sphenix)
AC_DEFINE([SPHENIX_SPH], [1], [SPHENIX SPH])
;;
diff --git a/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/kelvin_helmholtz.yml b/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/kelvin_helmholtz.yml
index e1ab89fd99..a7f6cbbd38 100644
--- a/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/kelvin_helmholtz.yml
+++ b/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/kelvin_helmholtz.yml
@@ -34,7 +34,7 @@ InitialConditions:
periodic: 1
Scheduler:
- max_top_level_cells: 50 # Maximal number of top-level cells in any dimension.
+ max_top_level_cells: 80 # Maximal number of top-level cells in any dimension.
# Parameters related to the equation of state
EoS:
diff --git a/src/debug.c b/src/debug.c
index 0b1538b5c4..2db8f69dff 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -74,6 +74,8 @@
#include "./hydro/Shadowswift/hydro_debug.h"
#elif defined(PLANETARY_SPH)
#include "./hydro/Planetary/hydro_debug.h"
+#elif defined(REMIX_SPH)
+#include "./hydro/REMIX/hydro_debug.h"
#elif defined(SPHENIX_SPH)
#include "./hydro/SPHENIX/hydro_debug.h"
#elif defined(GASOLINE_SPH)
diff --git a/src/hydro.h b/src/hydro.h
index 048be9f314..c82c402e81 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -70,6 +70,10 @@
#include "./hydro/Planetary/hydro.h"
#include "./hydro/Planetary/hydro_iact.h"
#define SPH_IMPLEMENTATION "Minimal version of SPH with multiple materials"
+#elif defined(REMIX_SPH)
+#include "./hydro/REMIX/hydro.h"
+#include "./hydro/REMIX/hydro_iact.h"
+#define SPH_IMPLEMENTATION "REMIX (Sandnes+ 2025)"
#elif defined(SPHENIX_SPH)
#include "./hydro/SPHENIX/hydro.h"
#include "./hydro/SPHENIX/hydro_iact.h"
@@ -89,9 +93,9 @@
/* Check whether this scheme implements the density checks */
#ifdef SWIFT_HYDRO_DENSITY_CHECKS
-#if !defined(SPHENIX_SPH) && !defined(PLANETARY_SPH)
+#if !defined(SPHENIX_SPH) && !defined(PLANETARY_SPH) && !defined(REMIX_SPH)
#error \
- "Can only use the hydro brute-force density checks with the SPHENIX or PLANETARY hydro schemes."
+ "Can only use the hydro brute-force density checks with the SPHENIX or PLANETARY or REMIX hydro schemes."
#endif
#endif
diff --git a/src/hydro/REMIX/hydro.h b/src/hydro/REMIX/hydro.h
new file mode 100644
index 0000000000..bb5517e78a
--- /dev/null
+++ b/src/hydro/REMIX/hydro.h
@@ -0,0 +1,981 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2024 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+ * 2024 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ * 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_REMIX_HYDRO_H
+#define SWIFT_REMIX_HYDRO_H
+
+/**
+ * @file REMIX/hydro.h
+ * @brief REMIX implementation of SPH (Sandnes et al. 2025)
+ */
+
+#include "adiabatic_index.h"
+#include "approx_math.h"
+#include "cosmology.h"
+#include "debug.h"
+#include "dimension.h"
+#include "entropy_floor.h"
+#include "equation_of_state.h"
+#include "hydro_kernels.h"
+#include "hydro_parameters.h"
+#include "hydro_properties.h"
+#include "hydro_space.h"
+#include "hydro_visc_difn.h"
+#include "kernel_hydro.h"
+#include "minmax.h"
+#include "pressure_floor.h"
+
+/**
+ * @brief Returns the comoving internal energy of a particle at the last
+ * time the particle was kicked.
+ *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy(const struct part *restrict p,
+ const struct xpart *restrict xp) {
+
+ return xp->u_full;
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle at the last
+ * time the particle was kicked.
+ *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable and converts it to
+ * physical coordinates.
+ *
+ * @param p The particle of interest.
+ * @param xp The extended data of the particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy(const struct part *restrict p,
+ const struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
+
+ return xp->u_full * cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Returns the comoving internal energy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_comoving_internal_energy(const struct part *restrict p) {
+
+ return p->u;
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_physical_internal_energy(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ return p->u * cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Returns the comoving pressure of a particle
+ *
+ * Computes the pressure based on the particle's properties.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
+ const struct part *restrict p) {
+
+ return gas_pressure_from_internal_energy(p->rho_evol, p->u, p->mat_id);
+}
+
+/**
+ * @brief Returns the physical pressure of a particle
+ *
+ * Computes the pressure based on the particle's properties and
+ * convert it to physical coordinates.
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
+ const struct part *restrict p, const struct cosmology *cosmo) {
+
+ return cosmo->a_factor_pressure *
+ gas_pressure_from_internal_energy(p->rho_evol, p->u, p->mat_id);
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
+ const struct part *restrict p, const struct xpart *restrict xp) {
+
+ return gas_entropy_from_internal_energy(p->rho_evol, xp->u_full, p->mat_id);
+}
+
+/**
+ * @brief Returns the physical entropy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable and converts it to
+ * physical coordinates.
+ *
+ * @param p The particle of interest
+ * @param xp The extended data of the particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
+ const struct part *restrict p, const struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
+
+ /* Note: no cosmological conversion required here with our choice of
+ * coordinates. */
+ return gas_entropy_from_internal_energy(p->rho_evol, xp->u_full, p->mat_id);
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_comoving_entropy(const struct part *restrict p) {
+
+ return gas_entropy_from_internal_energy(p->rho_evol, p->u, p->mat_id);
+}
+
+/**
+ * @brief Returns the physical entropy of a particle drifted to the
+ * current time.
+ *
+ * @param p The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_drifted_physical_entropy(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ /* Note: no cosmological conversion required here with our choice of
+ * coordinates. */
+ return gas_entropy_from_internal_energy(p->rho_evol, p->u, p->mat_id);
+}
+
+/**
+ * @brief Returns the comoving sound speed of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_soundspeed(const struct part *restrict p) {
+
+ return p->force.soundspeed;
+}
+
+/**
+ * @brief Returns the physical sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_soundspeed(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ return cosmo->a_factor_sound_speed * p->force.soundspeed;
+}
+
+/**
+ * @brief Returns the comoving density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
+ const struct part *restrict p) {
+
+ return p->rho;
+}
+
+/**
+ * @brief Returns the comoving density of a particle.
+ *
+ * @param p The particle of interest
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_physical_density(
+ const struct part *restrict p, const struct cosmology *cosmo) {
+
+ return cosmo->a3_inv * p->rho;
+}
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_mass(
+ const struct part *restrict p) {
+
+ return p->mass;
+}
+
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part *restrict p, float m) {
+
+ p->mass = m;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy_dt(const struct part *restrict p) {
+
+ return p->u_dt;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest
+ * @param cosmo Cosmology data structure
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_physical_internal_energy_dt(const struct part *restrict p,
+ const struct cosmology *cosmo) {
+
+ return p->u_dt * cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest.
+ * @param du_dt The new time derivative of the internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_comoving_internal_energy_dt(struct part *restrict p, float du_dt) {
+
+ p->u_dt = du_dt;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest.
+ * @param cosmo Cosmology data structure
+ * @param du_dt The new time derivative of the internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_physical_internal_energy_dt(struct part *restrict p,
+ const struct cosmology *cosmo,
+ float du_dt) {
+
+ p->u_dt = du_dt / cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Sets the physical entropy of a particle
+ *
+ * @param p The particle of interest.
+ * @param xp The extended particle data.
+ * @param cosmo Cosmology data structure
+ * @param entropy The physical entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_physical_entropy(
+ struct part *p, struct xpart *xp, const struct cosmology *cosmo,
+ const float entropy) {
+
+ /* Note there is no conversion from physical to comoving entropy */
+ const float comoving_entropy = entropy;
+ xp->u_full = gas_internal_energy_from_entropy(p->rho_evol, comoving_entropy,
+ p->mat_id);
+}
+
+/**
+ * @brief Sets the physical internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param xp The extended particle data.
+ * @param cosmo Cosmology data structure
+ * @param u The physical internal energy
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_physical_internal_energy(struct part *p, struct xpart *xp,
+ const struct cosmology *cosmo,
+ const float u) {
+
+ xp->u_full = u / cosmo->a_factor_internal_energy;
+}
+
+/**
+ * @brief Sets the drifted physical internal energy of a particle
+ *
+ * @param p The particle of interest.
+ * @param cosmo Cosmology data structure
+ * @param u The physical internal energy
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_drifted_physical_internal_energy(struct part *p,
+ const struct cosmology *cosmo,
+ const float u) {
+
+ p->u = u / cosmo->a_factor_internal_energy;
+
+ /* Now recompute the extra quantities */
+
+ /* Compute the sound speed */
+ const float pressure =
+ gas_pressure_from_internal_energy(p->rho_evol, p->u, p->mat_id);
+ const float soundspeed = hydro_get_comoving_soundspeed(p);
+
+ /* Update variables. */
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
+
+ p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
+}
+
+/**
+ * @brief Update the value of the viscosity alpha for the scheme.
+ *
+ * @param p the particle of interest
+ * @param alpha the new value for the viscosity coefficient.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
+ struct part *restrict p, float alpha) {
+ /* This scheme has fixed alpha */
+}
+
+/**
+ * @brief Update the value of the diffusive coefficients to the
+ * feedback reset value for the scheme.
+ *
+ * @param p the particle of interest
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_diffusive_feedback_reset(struct part *restrict p) {
+ /* This scheme has fixed alpha */
+}
+
+/**
+ * @brief Computes the hydro time-step of a given particle
+ *
+ * This function returns the time-step of a particle given its hydro-dynamical
+ * state. A typical time-step calculation would be the use of the CFL condition.
+ *
+ * @param p Pointer to the particle data
+ * @param xp Pointer to the extended particle data
+ * @param hydro_properties The SPH parameters
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
+ const struct part *restrict p, const struct xpart *restrict xp,
+ const struct hydro_props *restrict hydro_properties,
+ const struct cosmology *restrict cosmo) {
+
+ const float CFL_condition = hydro_properties->CFL_condition;
+
+ /* CFL condition */
+ const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
+ (cosmo->a_factor_sound_speed * p->force.v_sig);
+
+ return dt_cfl;
+}
+
+/**
+ * @brief returns the signal velocity
+ *
+ * @brief p the particle
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_signal_velocity(
+ const struct part *restrict p) {
+
+ return p->force.v_sig;
+}
+
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part *p, float dt) {}
+
+/**
+ * @brief Prepares a particle for the density calculation.
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the various density loop over neighbours. Typically, all fields of the
+ * density sub-structure of a particle get zeroed in here.
+ *
+ * @param p The particle to act upon
+ * @param hs #hydro_space containing hydro specific space information.
+ */
+__attribute__((always_inline)) INLINE static void hydro_init_part(
+ struct part *restrict p, const struct hydro_space *hs) {
+
+ p->density.wcount = 0.f;
+ p->density.wcount_dh = 0.f;
+ p->rho = 0.f;
+ p->density.rho_dh = 0.f;
+
+ hydro_init_part_extra_kernel(p);
+}
+
+/**
+ * @brief Finishes the density calculation.
+ *
+ * Multiplies the density and number of neighbours by the appropiate constants
+ * and add the self-contribution term.
+ * Additional quantities such as velocity gradients will also get the final
+ * terms added to them here.
+ *
+ * Also adds/multiplies the cosmological terms if need be.
+ *
+ * @param p The particle to act upon
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_density(
+ struct part *restrict p, const struct cosmology *cosmo) {
+
+ /* Some smoothing length multiples. */
+ 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 */
+ const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+
+ /* Final operation on the density (add self-contribution). */
+ p->rho += p->mass * kernel_root;
+ p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
+ p->density.wcount += kernel_root;
+ p->density.wcount_dh -= hydro_dimension * kernel_root;
+
+ /* Finish the calculation by inserting the missing h-factors */
+ p->rho *= h_inv_dim;
+ p->density.rho_dh *= h_inv_dim_plus_one;
+ p->density.wcount *= h_inv_dim;
+ p->density.wcount_dh *= h_inv_dim_plus_one;
+
+ hydro_end_density_extra_kernel(p);
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * In the desperate case where a particle has no neighbours (likely because
+ * of the h_max ceiling), set the particle fields to something sensible to avoid
+ * NaNs in the next calculations.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
+
+ /* Some smoothing length multiples. */
+ 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 */
+
+ /* Re-set problematic values */
+ p->rho = p->mass * kernel_root * h_inv_dim;
+ p->density.wcount = kernel_root * h_inv_dim;
+ p->density.rho_dh = 0.f;
+ p->density.wcount_dh = 0.f;
+
+ p->is_h_max = 1;
+ p->m0 = p->mass * kernel_root * h_inv_dim / p->rho_evol;
+}
+
+/**
+ * @brief Prepare a particle for the gradient calculation.
+ *
+ * This function is called after the density loop and before the gradient loop.
+ *
+ * We use it to set the physical timestep for the particle and to copy the
+ * actual velocities, which we need to boost our interfaces during the flux
+ * calculation. We also initialize the variables used for the time step
+ * calculation.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @param cosmo The cosmological model.
+ * @param hydro_props Hydrodynamic properties.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
+ struct part *restrict p, struct xpart *restrict xp,
+ 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 {
+ p->is_h_max = 0;
+ }
+
+ hydro_prepare_gradient_extra_kernel(p);
+ hydro_prepare_gradient_extra_visc_difn(p);
+}
+
+/**
+ * @brief Resets the variables that are required for a gradient calculation.
+ *
+ * This function is called after hydro_prepare_gradient.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_gradient(
+ struct part *restrict p) {}
+
+/**
+ * @brief Finishes the gradient calculation.
+ *
+ * Just a wrapper around hydro_gradients_finalize, which can be an empty method,
+ * in which case no gradients are used.
+ *
+ * This method also initializes the force loop variables.
+ *
+ * @param p The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_gradient(
+ struct part *p) {
+
+ hydro_end_gradient_extra_kernel(p);
+ hydro_end_gradient_extra_visc_difn(p);
+
+ /* Set the density to be used in the force loop to be the evolved density */
+ p->rho = p->rho_evol;
+}
+
+/**
+ * @brief Prepare a particle for the force calculation.
+ *
+ * This function is called in the ghost task to convert some quantities coming
+ * from the density loop over neighbours into quantities ready to be used in the
+ * force loop over neighbours. Quantities are typically read from the density
+ * sub-structure and written to the force sub-structure.
+ * Examples of calculations done here include the calculation of viscosity term
+ * constants, thermal conduction terms, hydro conversions, etc.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties.
+ * @param dt_alpha The time-step used to evolve non-cosmological quantities such
+ * as the artificial viscosity.
+ * @param dt_therm The time-step used to evolve hydrodynamical quantities.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+ const struct pressure_floor_props *pressure_floor, const float dt_alpha,
+ const float dt_therm) {
+
+ hydro_prepare_force_extra_kernel(p);
+
+#ifdef PLANETARY_FIXED_ENTROPY
+ /* Override the internal energy to satisfy the fixed entropy */
+ p->u = gas_internal_energy_from_entropy(p->rho_evol, p->s_fixed, p->mat_id);
+ xp->u_full = p->u;
+#endif
+
+ /* Compute the sound speed */
+ const float soundspeed =
+ gas_soundspeed_from_internal_energy(p->rho_evol, p->u, p->mat_id);
+
+ const float div_v = p->dv_norm_kernel[0][0] + p->dv_norm_kernel[1][1] +
+ p->dv_norm_kernel[2][2];
+ float curl_v[3];
+ curl_v[0] = p->dv_norm_kernel[1][2] - p->dv_norm_kernel[2][1];
+ curl_v[1] = p->dv_norm_kernel[2][0] - p->dv_norm_kernel[0][2];
+ curl_v[2] = p->dv_norm_kernel[0][1] - p->dv_norm_kernel[1][0];
+ const 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 (Sandnes+2025 Eqn. 34 with
+ * velocity gradients calculated by Eqn. 35) */
+ float balsara;
+ if (div_v == 0.f) {
+ balsara = 0.f;
+ } else {
+ balsara = fabsf(div_v) /
+ (fabsf(div_v) + mod_curl_v + 0.0001f * soundspeed / p->h);
+ }
+
+ /* Compute the pressure */
+ const float pressure =
+ gas_pressure_from_internal_energy(p->rho_evol, p->u, p->mat_id);
+ 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;
+}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * Resets all hydro acceleration and time derivative fields in preparation
+ * for the sums taking place in the various force tasks.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
+ struct part *restrict p) {
+
+ /* Reset the acceleration. */
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+
+ /* Reset the time derivatives. */
+ p->u_dt = 0.0f;
+ p->drho_dt = 0.0f;
+ p->force.h_dt = 0.0f;
+ p->force.v_sig = p->force.soundspeed;
+}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param p The particle.
+ * @param xp The extended data of this particle.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
+ struct part *restrict p, const struct xpart *restrict xp,
+ const struct cosmology *cosmo,
+ const struct pressure_floor_props *pressure_floor) {
+
+ /* Re-set the predicted velocities */
+ p->v[0] = xp->v_full[0];
+ p->v[1] = xp->v_full[1];
+ p->v[2] = xp->v_full[2];
+
+ /* 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;
+
+ /* Compute the pressure */
+ const float pressure =
+ gas_pressure_from_internal_energy(p->rho_evol, p->u, p->mat_id);
+
+ /* Compute the sound speed */
+ const float soundspeed =
+ gas_soundspeed_from_internal_energy(p->rho_evol, p->u, p->mat_id);
+
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
+
+ p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
+}
+
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * Additional hydrodynamic quantites are drifted forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * Note the different time-step sizes used for the different quantities as they
+ * include cosmological factors.
+ *
+ * @param p The particle.
+ * @param xp The extended data of the particle.
+ * @param dt_drift The drift time-step for positions.
+ * @param dt_therm The drift time-step for thermal quantities.
+ * @param dt_kick_grav The time-step for gravity quantities.
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme
+ * @param floor_props The properties of the entropy floor.
+ */
+__attribute__((always_inline)) INLINE static void hydro_predict_extra(
+ struct part *restrict p, const struct xpart *restrict xp, float dt_drift,
+ float dt_therm, float dt_kick_grav, const struct cosmology *cosmo,
+ const struct hydro_props *hydro_props,
+ const struct entropy_floor_properties *floor_props,
+ const struct pressure_floor_props *pressure_floor) {
+
+ /* Predict the internal energy and density */
+ p->u += p->u_dt * dt_therm;
+ p->rho_evol += p->drho_dt * dt_therm;
+
+ /* compute minimum SPH quantities */
+ 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 */
+ 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 =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+
+ p->u = max(p->u, min_u);
+
+ /* Predict smoothing length */
+ const float w1 = p->force.h_dt * h_inv * dt_drift;
+ if (fabsf(w1) < 0.2f)
+ p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+ else
+ p->h *= expf(w1);
+
+ const float floor_u = FLT_MIN;
+ p->u = max(p->u, floor_u);
+
+ /* Compute the new pressure */
+ const float pressure =
+ gas_pressure_from_internal_energy(p->rho_evol, p->u, p->mat_id);
+
+ /* Compute the new sound speed */
+ const float soundspeed =
+ gas_soundspeed_from_internal_energy(p->rho_evol, p->u, p->mat_id);
+
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
+
+ p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
+}
+
+/**
+ * @brief Finishes the force calculation.
+ *
+ * Multiplies the force and accelerations by the appropiate constants
+ * and add the self-contribution term. In most cases, there is little
+ * to do here.
+ *
+ * Cosmological terms are also added/multiplied here.
+ *
+ * @param p The particle to act upon
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_force(
+ struct part *restrict p, const struct cosmology *cosmo) {
+
+ p->force.h_dt *= p->h * hydro_dimension_inv;
+}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * Additional hydrodynamic quantites are kicked forward in time here. These
+ * include thermal quantities (thermal energy or total energy or entropy, ...).
+ *
+ * @param p The particle to act upon.
+ * @param xp The particle extended data to act upon.
+ * @param dt_therm The time-step for this kick (for thermodynamic quantities).
+ * @param dt_grav The time-step for this kick (for gravity quantities).
+ * @param dt_grav_mesh The time-step for this kick (mesh gravity).
+ * @param dt_hydro The time-step for this kick (for hydro quantities).
+ * @param dt_kick_corr The time-step for this kick (for gravity corrections).
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme
+ * @param floor_props The properties of the entropy floor.
+ */
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+ struct part *restrict p, struct xpart *restrict xp, float dt_therm,
+ float dt_grav, float dt_grav_mesh, float dt_hydro, float dt_kick_corr,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+ const struct entropy_floor_properties *floor_props) {
+
+ /* Integrate the internal energy forward in time */
+ const float delta_u = p->u_dt * dt_therm;
+
+ /* Do not decrease the energy by more than a factor of 2*/
+ xp->u_full = max(xp->u_full + delta_u, 0.5f * xp->u_full);
+
+ /* Check against absolute minimum */
+ const float min_u =
+ hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+ const float floor_u = FLT_MIN;
+
+ /* Take highest of both limits */
+ const float energy_min = max(min_u, floor_u);
+
+ if (xp->u_full < energy_min) {
+ xp->u_full = energy_min;
+ 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 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 */
+ const float floor_rho = p->mass * kernel_root * h_inv_dim;
+ if (xp->rho_evol_full < floor_rho) {
+ xp->rho_evol_full = floor_rho;
+ p->drho_dt = 0.f;
+ }
+}
+
+/**
+ * @brief Converts hydro quantity of a particle at the start of a run
+ *
+ * This function is called once at the end of the engine_init_particle()
+ * routine (at the start of a calculation) after the densities of
+ * particles have been computed.
+ * This can be used to convert internal energy into entropy for instance.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle to act upon
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+ const struct pressure_floor_props *pressure_floor) {
+
+ /* Compute the pressure */
+ const float pressure =
+ gas_pressure_from_internal_energy(p->rho_evol, p->u, p->mat_id);
+
+ /* Compute the sound speed */
+ const float soundspeed =
+ gas_soundspeed_from_internal_energy(p->rho_evol, p->u, p->mat_id);
+
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
+}
+
+/**
+ * @brief Initialises the particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions or assignments between the particle
+ * and extended particle fields.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+ struct part *restrict p, struct xpart *restrict xp) {
+
+ p->time_bin = 0;
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
+ xp->u_full = p->u;
+
+ p->rho_evol = p->rho;
+ xp->rho_evol_full = p->rho_evol;
+
+ hydro_reset_acceleration(p);
+ hydro_init_part(p, NULL);
+}
+
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part *p, float u_init) {
+
+ p->u = u_init;
+}
+
+/**
+ * @brief Operations performed when a particle gets removed from the
+ * simulation volume.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ * @param time The simulation time.
+ */
+__attribute__((always_inline)) INLINE static void hydro_remove_part(
+ const struct part *p, const struct xpart *xp, const double time) {
+
+ /* Print the particle info as csv to facilitate later analysis, e.g. with
+ * grep '## Removed' -A 1 --no-group-separator output.txt > removed.txt
+ */
+ printf(
+ "## Removed particle: "
+ "id, x, y, z, vx, vy, vz, m, u, P, rho, h, mat_id, time \n"
+ "%lld, %.7g, %.7g, %.7g, %.7g, %.7g, %.7g, "
+ "%.7g, %.7g, %.7g, %.7g, %.7g, %d, %.7g \n",
+ p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->mass,
+ p->u, p->force.pressure, p->rho, p->h, p->mat_id, time);
+}
+
+#endif /* SWIFT_REMIX_HYDRO_H */
diff --git a/src/hydro/REMIX/hydro_debug.h b/src/hydro/REMIX/hydro_debug.h
new file mode 100644
index 0000000000..5fbd6c6efd
--- /dev/null
+++ b/src/hydro/REMIX/hydro_debug.h
@@ -0,0 +1,53 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2024 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+ * 2024 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ * 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_REMIX_HYDRO_DEBUG_H
+#define SWIFT_REMIX_HYDRO_DEBUG_H
+
+/**
+ * @file REMIX/hydro_debug.h
+ * @brief Debugging routines for REMIX SPH.
+ */
+
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+ const struct part* p, const struct xpart* xp) {
+ warning("[PID%lld] part:", p->id);
+ warning(
+ "[PID%lld] "
+ "x=[%.6g, %.6g, %.6g], v=[%.3g, %.3g, %.3g], "
+ "a=[%.3g, %.3g, %.3g], "
+ "m=%.3g, u=%.3g, du/dt=%.3g, P=%.3g, c_s=%.3g, "
+ "v_sig=%.3g, h=%.3g, dh/dt=%.3g, wcount=%.3g, rho=%.3g, "
+ "dh_drho=%.3g, rho_evol=%.3g, time_bin=%d, wakeup=%d, mat_id=%d",
+ p->id, 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->mass, p->u, p->u_dt,
+ hydro_get_comoving_pressure(p), p->force.soundspeed, p->force.v_sig, p->h,
+ p->force.h_dt, p->density.wcount, p->rho, p->density.rho_dh, p->rho_evol,
+ p->time_bin, p->limiter_data.wakeup, p->mat_id);
+ if (xp != NULL) {
+ warning("[PID%lld] xpart:", p->id);
+ warning(
+ "[PID%lld] "
+ "v_full=[%.3g, %.3g, %.3g]",
+ p->id, xp->v_full[0], xp->v_full[1], xp->v_full[2]);
+ }
+}
+
+#endif /* SWIFT_REMIX_HYDRO_DEBUG_H */
diff --git a/src/hydro/REMIX/hydro_iact.h b/src/hydro/REMIX/hydro_iact.h
new file mode 100644
index 0000000000..a1b27788ef
--- /dev/null
+++ b/src/hydro/REMIX/hydro_iact.h
@@ -0,0 +1,539 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2024 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+ * 2024 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ * 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_REMIX_HYDRO_IACT_H
+#define SWIFT_REMIX_HYDRO_IACT_H
+
+/**
+ * @file REMIX/hydro_iact.h
+ * @brief REMIX implementation of SPH (Sandnes et al. 2025)
+ */
+
+#include "adiabatic_index.h"
+#include "const.h"
+#include "hydro_kernels.h"
+#include "hydro_parameters.h"
+#include "hydro_visc_difn.h"
+#include "math.h"
+#include "minmax.h"
+
+/**
+ * @brief Density interaction between two particles.
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_density(
+ const float r2, const float dx[3], const float hi, const float hj,
+ struct part *restrict pi, struct part *restrict pj, const float a,
+ const float H) {
+
+ float wi, wj, wi_dx, wj_dx;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (pi->time_bin >= time_bin_inhibited)
+ error("Inhibited pi in interaction function!");
+ if (pj->time_bin >= time_bin_inhibited)
+ error("Inhibited pj in interaction function!");
+#endif
+
+ /* Get r */
+ const float r = sqrtf(r2);
+
+ /* Get the masses. */
+ const float mi = pi->mass;
+ const float mj = pj->mass;
+
+ /* Compute density of pi. */
+ const float hi_inv = 1.f / hi;
+ const float ui = r * hi_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ pi->rho += mj * wi;
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+ /* Compute density of pj. */
+ const float hj_inv = 1.f / hj;
+ const float uj = r * hj_inv;
+ kernel_deval(uj, &wj, &wj_dx);
+
+ pj->rho += mi * wj;
+ pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
+ pj->density.wcount += wj;
+ pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
+
+ hydro_runner_iact_density_extra_kernel(pi, pj, dx, wi, wj, wi_dx, wj_dx);
+}
+
+/**
+ * @brief Density interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
+ const float r2, const float dx[3], const float hi, const float hj,
+ struct part *restrict pi, const struct part *restrict pj, const float a,
+ const float H) {
+
+ float wi, wi_dx;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (pi->time_bin >= time_bin_inhibited)
+ error("Inhibited pi in interaction function!");
+ if (pj->time_bin >= time_bin_inhibited)
+ error("Inhibited pj in interaction function!");
+#endif
+
+ /* Get the masses. */
+ const float mj = pj->mass;
+
+ /* Get r. */
+ const float r = sqrtf(r2);
+
+ const float h_inv = 1.f / hi;
+ const float ui = r * h_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ pi->rho += mj * wi;
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+ hydro_runner_iact_nonsym_density_extra_kernel(pi, pj, dx, wi, wi_dx);
+}
+
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j
+ *
+ * This method wraps around hydro_gradients_collect, which can be an empty
+ * method, in which case no gradients are used.
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_gradient(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ float wi, wj, wi_dx, wj_dx;
+
+ /* Get r. */
+ const float r = sqrtf(r2);
+
+ /* Compute kernel of pi. */
+ const float hi_inv = 1.f / hi;
+ const float ui = r * hi_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ /* Compute kernel of pj. */
+ const float hj_inv = 1.f / hj;
+ const float uj = r * hj_inv;
+ 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_visc_difn(pi, pj, dx, wi, wj, wi_dx, wj_dx);
+}
+
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j:
+ * non-symmetric version
+ *
+ * This method wraps around hydro_gradients_nonsym_collect, which can be an
+ * empty method, in which case no gradients are used.
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ float wi, wj, wi_dx, wj_dx;
+
+ /* Get r. */
+ const float r = sqrtf(r2);
+
+ /* Compute kernel of pi. */
+ const float h_inv = 1.f / hi;
+ const float ui = r * h_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ /* Compute kernel of pj. */
+ const float hj_inv = 1.f / hj;
+ const float uj = r * hj_inv;
+ kernel_deval(uj, &wj, &wj_dx);
+
+ hydro_runner_iact_nonsym_gradient_extra_kernel(pi, pj, dx, wi, wj, wi_dx,
+ wj_dx);
+ hydro_runner_iact_nonsym_gradient_extra_visc_difn(pi, pj, dx, wi, wi_dx);
+}
+
+/**
+ * @brief Force interaction between two particles.
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_force(
+ const float r2, const float dx[3], const float hi, const float hj,
+ struct part *restrict pi, struct part *restrict pj, const float a,
+ const float H) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (pi->time_bin >= time_bin_inhibited)
+ error("Inhibited pi in interaction function!");
+ if (pj->time_bin >= time_bin_inhibited)
+ error("Inhibited pj in interaction function!");
+#endif
+
+ /* Get r. */
+ const float r = sqrtf(r2);
+
+ /* Recover some data */
+ const float mi = pi->mass;
+ const float mj = pj->mass;
+ const float rhoi = pi->rho;
+ const float rhoj = pj->rho;
+ const float rhoi_inv = 1.f / rhoi;
+ const float rhoj_inv = 1.f / rhoj;
+ const float pressurei = pi->force.pressure;
+ const float pressurej = pj->force.pressure;
+
+ /* Get the kernel for hi. */
+ const float hi_inv = 1.0f / hi;
+ const float xi = r * hi_inv;
+ float wi, wi_dx;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ /* Get the kernel for hj. */
+ const float hj_inv = 1.0f / hj;
+ const float xj = r * hj_inv;
+ float wj, wj_dx;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ /* 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);
+
+ /* Antisymmetric kernel grad term for conservation of momentum and energy */
+ G_mean[0] = 0.5f * (Gi[0] - Gj[0]);
+ G_mean[1] = 0.5f * (Gi[1] - Gj[1]);
+ G_mean[2] = 0.5f * (Gi[2] - Gj[2]);
+
+ /* 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);
+
+ /* Pressure terms to be used in evolution equations */
+ const float P_i_term = pressurei * rhoi_inv * rhoj_inv;
+ const float P_j_term = pressurej * rhoi_inv * rhoj_inv;
+ const float Q_i_term = Qi * rhoi_inv * rhoj_inv;
+ const float Q_j_term = Qj * rhoi_inv * rhoj_inv;
+
+ /* Use the force Luke! */
+ pi->a_hydro[0] -=
+ mj * (P_i_term + P_j_term + Q_i_term + Q_j_term) * G_mean[0];
+ pi->a_hydro[1] -=
+ mj * (P_i_term + P_j_term + Q_i_term + Q_j_term) * G_mean[1];
+ pi->a_hydro[2] -=
+ mj * (P_i_term + P_j_term + Q_i_term + Q_j_term) * G_mean[2];
+
+ pj->a_hydro[0] +=
+ mi * (P_i_term + P_j_term + Q_i_term + Q_j_term) * G_mean[0];
+ pj->a_hydro[1] +=
+ mi * (P_i_term + P_j_term + Q_i_term + Q_j_term) * G_mean[1];
+ pj->a_hydro[2] +=
+ mi * (P_i_term + P_j_term + Q_i_term + Q_j_term) * G_mean[2];
+
+ /* 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];
+
+ /* Internal energy time derivative */
+ 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 * (rhoi * rhoj_inv) * dvdotG;
+ pj->drho_dt += mi * (rhoj * rhoi_inv) * dvdotG;
+
+ /* Get the time derivative for h. */
+ pi->force.h_dt -= mj * dvdotG * rhoj_inv;
+ pj->force.h_dt -= mi * dvdotG * rhoi_inv;
+
+ const float v_sig = visc_signal_velocity;
+
+ /* Update the signal velocity. */
+ pi->force.v_sig = max(pi->force.v_sig, v_sig);
+ pj->force.v_sig = max(pj->force.v_sig, v_sig);
+
+ if ((pi->is_h_max) || (pj->is_h_max)) {
+ /* Do not add diffusion or normalising terms if either particle has h=h_max
+ */
+ return;
+ }
+
+ /* Calculate some quantities */
+ const float mean_rho = 0.5f * (rhoi + rhoj);
+ const float mean_rho_inv = 1.f / mean_rho;
+ const float mean_balsara = 0.5f * (pi->force.balsara + pj->force.balsara);
+ const float mod_G = sqrtf(G_mean[0] * G_mean[0] + G_mean[1] * G_mean[1] +
+ G_mean[2] * G_mean[2]);
+ const float v_sig_norm = sqrtf((pi->v[0] - pj->v[0]) * (pi->v[0] - pj->v[0]) +
+ (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 * rhoi - rhoi) * mod_G * mean_rho_inv;
+ float drho_dt_norm_and_difn_j = alpha_norm * mi * v_sig_norm *
+ pj->force.vac_switch *
+ (pj->m0 * rhoj - rhoj) * mod_G * mean_rho_inv;
+
+ /* Only include diffusion for same-material particle pair */
+ if (pi->mat_id == pj->mat_id) {
+ /* 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);
+ const float v_sig_difn = difn_signal_velocity;
+
+ /* Calculate artificial diffusion of internal energy (Sandnes+2025 Eqn. 42)
+ */
+ const 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_inv;
+ const float du_dt_difn_j = -(a_difn_u + b_difn_u * mean_balsara) * mi *
+ v_sig_difn * (utilde_j - utilde_i) * mod_G *
+ mean_rho_inv;
+
+ /* 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 *
+ (rhoi * rhoj_inv) * v_sig_difn *
+ (rhotilde_i - rhotilde_j) * mod_G * mean_rho_inv;
+ drho_dt_norm_and_difn_j += -(a_difn_rho + b_difn_rho * mean_balsara) * mi *
+ (rhoj * rhoi_inv) * v_sig_difn *
+ (rhotilde_j - rhotilde_i) * mod_G * mean_rho_inv;
+ }
+
+ /* 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;
+}
+
+/**
+ * @brief Force interaction between two particles (non-symmetric).
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
+ const float r2, const float dx[3], const float hi, const float hj,
+ struct part *restrict pi, const struct part *restrict pj, const float a,
+ const float H) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (pi->time_bin >= time_bin_inhibited)
+ error("Inhibited pi in interaction function!");
+ if (pj->time_bin >= time_bin_inhibited)
+ error("Inhibited pj in interaction function!");
+#endif
+
+ /* Get r. */
+ const float r = sqrtf(r2);
+
+ /* Recover some data */
+ const float mj = pj->mass;
+ const float rhoi = pi->rho;
+ const float rhoj = pj->rho;
+ const float rhoi_inv = 1.f / rhoi;
+ const float rhoj_inv = 1.f / rhoj;
+ const float pressurei = pi->force.pressure;
+ const float pressurej = pj->force.pressure;
+
+ /* Get the kernel for hi. */
+ const float hi_inv = 1.0f / hi;
+ const float xi = r * hi_inv;
+ float wi, wi_dx;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ /* Get the kernel for hj. */
+ const float hj_inv = 1.0f / hj;
+ const float xj = r * hj_inv;
+ float wj, wj_dx;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ /* 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);
+
+ /* Antisymmetric kernel grad term for conservation of momentum and energy */
+ G_mean[0] = 0.5f * (Gi[0] - Gj[0]);
+ G_mean[1] = 0.5f * (Gi[1] - Gj[1]);
+ G_mean[2] = 0.5f * (Gi[2] - Gj[2]);
+
+ /* 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);
+
+ /* Pressure terms to be used in evolution equations */
+ const float P_i_term = pressurei * rhoi_inv * rhoj_inv;
+ const float P_j_term = pressurej * rhoi_inv * rhoj_inv;
+ const float Q_i_term = Qi * rhoi_inv * rhoj_inv;
+ const float Q_j_term = Qj * rhoi_inv * rhoj_inv;
+
+ /* Use the force Luke! */
+ pi->a_hydro[0] -=
+ mj * (P_i_term + P_j_term + Q_i_term + Q_j_term) * G_mean[0];
+ pi->a_hydro[1] -=
+ mj * (P_i_term + P_j_term + Q_i_term + Q_j_term) * G_mean[1];
+ pi->a_hydro[2] -=
+ mj * (P_i_term + P_j_term + Q_i_term + Q_j_term) * G_mean[2];
+
+ /* 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];
+
+ /* Internal energy time derivative */
+ pi->u_dt += mj * (P_i_term + Q_i_term) * dvdotG;
+
+ /* Density time derivative */
+ pi->drho_dt += mj * (rhoi * rhoj_inv) * dvdotG;
+
+ /* Get the time derivative for h. */
+ pi->force.h_dt -= mj * dvdotG * rhoj_inv;
+
+ const float v_sig = visc_signal_velocity;
+
+ /* Update the signal velocity. */
+ pi->force.v_sig = max(pi->force.v_sig, v_sig);
+
+ if ((pi->is_h_max) || (pj->is_h_max)) {
+ /* Do not add diffusion or normalising terms if either particle has h=h_max
+ */
+ return;
+ }
+
+ /* Calculate some quantities */
+ const float mean_rho = 0.5f * (rhoi + rhoj);
+ const float mean_rho_inv = 1.f / mean_rho;
+ const float mean_balsara = 0.5f * (pi->force.balsara + pj->force.balsara);
+ const float mod_G = sqrtf(G_mean[0] * G_mean[0] + G_mean[1] * G_mean[1] +
+ G_mean[2] * G_mean[2]);
+
+ const float v_sig_norm = sqrtf((pi->v[0] - pj->v[0]) * (pi->v[0] - pj->v[0]) +
+ (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 * rhoi - rhoi) * mod_G * mean_rho_inv;
+
+ /* Only include diffusion for same-material particle pair */
+ if (pi->mat_id == pj->mat_id) {
+ /* 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);
+ const float v_sig_difn = difn_signal_velocity;
+
+ /* Calculate artificial diffusion of internal energy (Sandnes+2025 Eqn. 42)
+ */
+ const 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_inv;
+
+ /* 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 *
+ (rhoi * rhoj_inv) * v_sig_difn *
+ (rhotilde_i - rhotilde_j) * mod_G * mean_rho_inv;
+ }
+
+ /* Add normalising term and artificial diffusion to evolution of density */
+ pi->drho_dt += drho_dt_norm_and_difn_i;
+}
+
+#endif /* SWIFT_REMIX_HYDRO_IACT_H */
diff --git a/src/hydro/REMIX/hydro_io.h b/src/hydro/REMIX/hydro_io.h
new file mode 100644
index 0000000000..1972ef029d
--- /dev/null
+++ b/src/hydro/REMIX/hydro_io.h
@@ -0,0 +1,251 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2024 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+ * 2024 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ * 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_REMIX_HYDRO_IO_H
+#define SWIFT_REMIX_HYDRO_IO_H
+
+/**
+ * @file REMIX/hydro_io.h
+ * @brief REMIX implementation of SPH (Sandnes et al. 2025) i/o routines
+ */
+
+#include "adiabatic_index.h"
+#include "hydro.h"
+#include "hydro_parameters.h"
+#include "io_properties.h"
+#include "kernel_hydro.h"
+
+/**
+ * @brief Specifies which particle fields to read from a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+INLINE static void hydro_read_particles(struct part* parts,
+ struct io_props* list,
+ int* num_fields) {
+
+#ifdef PLANETARY_FIXED_ENTROPY
+ *num_fields = 10;
+#else
+ *num_fields = 9;
+#endif
+
+ /* List what we want to read */
+ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+ UNIT_CONV_LENGTH, parts, x);
+ list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+ UNIT_CONV_SPEED, parts, v);
+ list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+ parts, mass);
+ list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_LENGTH, parts, h);
+ list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
+ list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+ UNIT_CONV_NO_UNITS, parts, id);
+ 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, COMPULSORY,
+ UNIT_CONV_DENSITY, parts, rho);
+ list[8] = io_make_input_field("MaterialIDs", INT, 1, COMPULSORY,
+ UNIT_CONV_NO_UNITS, parts, mat_id);
+#ifdef PLANETARY_FIXED_ENTROPY
+ list[9] = io_make_input_field("Entropies", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS, parts,
+ s_fixed);
+#endif
+}
+
+INLINE static void convert_S(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ ret[0] = hydro_get_comoving_entropy(p, xp);
+}
+
+INLINE static void convert_P(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ ret[0] = hydro_get_comoving_pressure(p);
+}
+
+INLINE static void convert_part_pos(const struct engine* e,
+ const struct part* p,
+ const struct xpart* xp, double* ret) {
+ const struct space* s = e->s;
+ if (s->periodic) {
+ ret[0] = box_wrap(p->x[0], 0.0, s->dim[0]);
+ ret[1] = box_wrap(p->x[1], 0.0, s->dim[1]);
+ ret[2] = box_wrap(p->x[2], 0.0, s->dim[2]);
+ } else {
+ ret[0] = p->x[0];
+ ret[1] = p->x[1];
+ ret[2] = p->x[2];
+ }
+ if (e->snapshot_use_delta_from_edge) {
+ ret[0] = min(ret[0], s->dim[0] - e->snapshot_delta_from_edge);
+ ret[1] = min(ret[1], s->dim[1] - e->snapshot_delta_from_edge);
+ ret[2] = min(ret[2], s->dim[2] - e->snapshot_delta_from_edge);
+ }
+}
+
+INLINE static void convert_part_vel(const struct engine* e,
+ const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+ const float dt_kick_grav_mesh = e->dt_kick_grav_mesh_for_io;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time (hydro term)*/
+ ret[0] = xp->v_full[0] + p->a_hydro[0] * dt_kick_hydro;
+ ret[1] = xp->v_full[1] + p->a_hydro[1] * dt_kick_hydro;
+ ret[2] = xp->v_full[2] + p->a_hydro[2] * dt_kick_hydro;
+
+ /* Add the gravity term */
+ if (p->gpart != NULL) {
+ ret[0] += p->gpart->a_grav[0] * dt_kick_grav;
+ ret[1] += p->gpart->a_grav[1] * dt_kick_grav;
+ ret[2] += p->gpart->a_grav[2] * dt_kick_grav;
+ }
+
+ /* And the mesh gravity term */
+ if (p->gpart != NULL) {
+ ret[0] += p->gpart->a_grav_mesh[0] * dt_kick_grav_mesh;
+ ret[1] += p->gpart->a_grav_mesh[1] * dt_kick_grav_mesh;
+ ret[2] += p->gpart->a_grav_mesh[2] * dt_kick_grav_mesh;
+ }
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a_inv;
+ ret[1] *= cosmo->a_inv;
+ ret[2] *= cosmo->a_inv;
+}
+
+INLINE static void convert_part_potential(const struct engine* e,
+ const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ if (p->gpart != NULL)
+ ret[0] = gravity_get_comoving_potential(p->gpart);
+ else
+ ret[0] = 0.f;
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param xparts The extended particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+INLINE static void hydro_write_particles(const struct part* parts,
+ const struct xpart* xparts,
+ struct io_props* list,
+ int* num_fields) {
+
+ *num_fields = 11;
+
+ /* List what we want to write */
+ list[0] = io_make_output_field_convert_part(
+ "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, parts, xparts,
+ convert_part_pos, "Positions of the particles");
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, parts, xparts,
+ convert_part_vel, "Velocities of the particles");
+ list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f, parts,
+ mass, "Masses of the particles");
+ list[3] = io_make_output_field(
+ "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, parts, h,
+ "Smoothing lengths (FWHM of the kernel) of the particles");
+ list[4] = io_make_output_field(
+ "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS,
+ -3.f * hydro_gamma_minus_one, parts, u,
+ "Thermal energies per unit mass of the particles");
+ list[5] =
+ io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
+ parts, id, "Unique IDs of the particles");
+ list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
+ parts, rho, "Densities of the particles");
+ list[7] = io_make_output_field_convert_part(
+ "Entropies", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, 0.f, parts,
+ xparts, convert_S, "Entropies per unit mass of the particles");
+ list[8] =
+ io_make_output_field("MaterialIDs", INT, 1, UNIT_CONV_NO_UNITS, 0.f,
+ parts, mat_id, "Material IDs of the particles");
+ list[9] = io_make_output_field_convert_part(
+ "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, -3.f * hydro_gamma, parts,
+ xparts, convert_P, "Pressures of the particles");
+ list[10] = io_make_output_field_convert_part(
+ "Potentials", FLOAT, 1, UNIT_CONV_POTENTIAL, 0.f, parts, xparts,
+ convert_part_potential, "Gravitational potentials of the particles");
+}
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grpsph The HDF5 group in which to write
+ */
+INLINE static void hydro_write_flavour(hid_t h_grpsph) {
+
+ io_write_attribute_s(h_grpsph, "Density estimate", "Evolved in time");
+ io_write_attribute_s(h_grpsph, "EoM free functions", "Evolved densities");
+ io_write_attribute_s(h_grpsph, "Kernel gradients",
+ "Linear-order reproducing kernels with grad-h terms");
+ io_write_attribute_s(h_grpsph, "Vacuum boundary treatment",
+ "As presented in Sandnes et al. (2025)");
+ io_write_attribute_s(
+ h_grpsph, "Artificial viscosity",
+ "With linear reconstruction, van Leer slope limiter, Balsara switch");
+ io_write_attribute_s(
+ h_grpsph, "Artificial diffusion",
+ "With linear reconstruction, van Leer slope limiter, Balsara switch");
+ io_write_attribute_s(h_grpsph, "Normalising term",
+ "As presented in Sandnes et al. (2025)");
+}
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+INLINE static int writeEntropyFlag(void) { return 0; }
+
+#endif /* SWIFT_REMIX_HYDRO_IO_H */
diff --git a/src/hydro/REMIX/hydro_kernels.h b/src/hydro/REMIX/hydro_kernels.h
new file mode 100644
index 0000000000..8d4edcb69f
--- /dev/null
+++ b/src/hydro/REMIX/hydro_kernels.h
@@ -0,0 +1,644 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2024 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+ * 2024 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_REMIX_HYDRO_KERNELS_H
+#define SWIFT_REMIX_HYDRO_KERNELS_H
+
+/**
+ * @file REMIX/hydro_kernels.h
+ * @brief Utilities for REMIX hydro kernels.
+ */
+
+#include "const.h"
+#include "hydro_parameters.h"
+#include "math.h"
+
+/**
+ * @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;
+}
+
+/**
+ * @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,
+ 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;
+ const float volume_i = pi->mass / pi->rho_evol;
+ const float volume_j = pj->mass / pj->rho_evol;
+
+ /* 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 += wi * volume_j;
+ pj->m0 += wj * volume_i;
+
+ pi->grad_m0[0] += dx[0] * wi_dx * r_inv * volume_j;
+ pi->grad_m0[1] += dx[1] * wi_dx * r_inv * volume_j;
+ pi->grad_m0[2] += dx[2] * wi_dx * r_inv * volume_j;
+
+ pj->grad_m0[0] += -dx[0] * wj_dx * r_inv * volume_i;
+ pj->grad_m0[1] += -dx[1] * wj_dx * r_inv * volume_i;
+ pj->grad_m0[2] += -dx[2] * wj_dx * r_inv * volume_i;
+}
+
+/**
+ * @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,
+ 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;
+ const float volume_j = pj->mass / pj->rho_evol;
+
+ /* 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 += wi * volume_j;
+
+ pi->grad_m0[0] += dx[0] * wi_dx * r_inv * volume_j;
+ pi->grad_m0[1] += dx[1] * wi_dx * r_inv * volume_j;
+ pi->grad_m0[2] += dx[2] * wi_dx * r_inv * volume_j;
+}
+
+/**
+ * @brief Finishes extra kernel parts of the density calculation.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_end_density_extra_kernel(struct part *restrict p) {
+
+ 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 */
+ const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+
+ /* 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;
+ p->grad_m0[0] *= h_inv_dim_plus_one;
+ p->grad_m0[1] *= h_inv_dim_plus_one;
+ p->grad_m0[2] *= h_inv_dim_plus_one;
+}
+
+/**
+ * @brief Prepares extra kernel 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_kernel(struct part *restrict p) {
+
+ /* Initialise geometric moment matrices (ij-mean, for linear-order repr.
+ * kernel) */
+ zero_sym_matrix(&p->gradient.m2_bar);
+ zero_sym_matrix(&p->gradient.grad_m2_bar[0]);
+ zero_sym_matrix(&p->gradient.grad_m2_bar[1]);
+ zero_sym_matrix(&p->gradient.grad_m2_bar[2]);
+ zero_sym_matrix(&p->gradient.grad_m2_bar_gradhterm);
+
+ /* 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;
+ memset(p->gradient.m1_bar, 0.f, 3 * sizeof(float));
+ memset(p->gradient.grad_m0_bar, 0.f, 3 * sizeof(float));
+ memset(p->gradient.grad_m1_bar_gradhterm, 0.f, 3 * sizeof(float));
+ memset(p->gradient.grad_m1_bar, 0.f, 3 * 3 * sizeof(float));
+}
+
+/**
+ * @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,
+ 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;
+
+ const float hi = pi->h;
+ const float hi_inv = 1.0f / hi; /* 1/h */
+ const float hi_inv_dim = pow_dimension(hi_inv); /* 1/h^d */
+ const float hi_inv_dim_plus_one = hi_inv_dim * hi_inv; /* 1/h^(d+1) */
+
+ const float hj = pj->h;
+ const float hj_inv = 1.0f / hj; /* 1/h */
+ 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 */
+ const float volume_i = pi->mass / pi->rho_evol;
+ const float volume_j = pj->mass / pj->rho_evol;
+
+ /* Mean ij kernels and gradients */
+ const float wi_term = 0.5f * (wi * hi_inv_dim + wj * hj_inv_dim);
+ const float wj_term = wi_term;
+ float wi_dx_term[3], wj_dx_term[3];
+ const float mean_dw_dr =
+ 0.5f * (wi_dx * hi_inv_dim_plus_one + wj_dx * hj_inv_dim_plus_one);
+ wi_dx_term[0] = dx[0] * r_inv * mean_dw_dr;
+ wi_dx_term[1] = dx[1] * r_inv * mean_dw_dr;
+ wi_dx_term[2] = dx[2] * r_inv * mean_dw_dr;
+ wj_dx_term[0] = -wi_dx_term[0];
+ wj_dx_term[1] = -wi_dx_term[1];
+ wj_dx_term[2] = -wi_dx_term[2];
+
+ /* Grad-h term, dW/dh */
+ const float wi_dx_gradhterm = -0.5f *
+ (hydro_dimension * wi + (r * hi_inv) * wi_dx) *
+ hi_inv_dim_plus_one;
+ const float wj_dx_gradhterm = -0.5f *
+ (hydro_dimension * wj + (r * hj_inv) * wj_dx) *
+ hj_inv_dim_plus_one;
+
+ /* 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;
+
+ pi->gradient.grad_m0_bar_gradhterm += wi_dx_gradhterm * volume_j;
+ pj->gradient.grad_m0_bar_gradhterm += wj_dx_gradhterm * volume_i;
+ for (int i = 0; i < 3; i++) {
+ pi->gradient.m1_bar[i] += dx[i] * wi_term * volume_j;
+ pj->gradient.m1_bar[i] += -dx[i] * wj_term * volume_i;
+
+ pi->gradient.grad_m0_bar[i] += wi_dx_term[i] * volume_j;
+ pj->gradient.grad_m0_bar[i] += wj_dx_term[i] * volume_i;
+
+ pi->gradient.grad_m1_bar_gradhterm[i] += dx[i] * wi_dx_gradhterm * volume_j;
+ pj->gradient.grad_m1_bar_gradhterm[i] +=
+ -dx[i] * wj_dx_gradhterm * volume_i;
+
+ for (int j = 0; j < 3; j++) {
+ pi->gradient.grad_m1_bar[i][j] += dx[j] * wi_dx_term[i] * volume_j;
+ pj->gradient.grad_m1_bar[i][j] += -dx[j] * wj_dx_term[i] * volume_i;
+ }
+ }
+
+ for (int j = 0; j < 3; j++) {
+ for (int k = j; k < 3; k++) {
+ int i = (j == k) ? j : 2 + j + k;
+ pi->gradient.m2_bar.elements[i] += dx[j] * dx[k] * wi_term * volume_j;
+ pj->gradient.m2_bar.elements[i] += dx[j] * dx[k] * wj_term * volume_i;
+
+ pi->gradient.grad_m2_bar_gradhterm.elements[i] +=
+ dx[j] * dx[k] * wi_dx_gradhterm * volume_j;
+ pj->gradient.grad_m2_bar_gradhterm.elements[i] +=
+ dx[j] * dx[k] * wj_dx_gradhterm * volume_i;
+
+ for (int l = 0; l < 3; l++) {
+ pi->gradient.grad_m2_bar[l].elements[i] +=
+ dx[j] * dx[k] * wi_dx_term[l] * volume_j;
+ pj->gradient.grad_m2_bar[l].elements[i] +=
+ dx[j] * dx[k] * wj_dx_term[l] * volume_i;
+ }
+ }
+ }
+}
+
+/**
+ * @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;
+
+ const float hi = pi->h;
+ const float hi_inv = 1.0f / hi; /* 1/h */
+ const float hi_inv_dim = pow_dimension(hi_inv); /* 1/h^d */
+ const float hi_inv_dim_plus_one = hi_inv_dim * hi_inv; /* 1/h^(d+1) */
+
+ const float hj = pj->h;
+ const float hj_inv = 1.0f / hj; /* 1/h */
+ 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 */
+ const float volume_j = pj->mass / pj->rho_evol;
+
+ /* Mean ij kernel and gradients */
+ const float wi_term = 0.5f * (wi * hi_inv_dim + wj * hj_inv_dim);
+ float wi_dx_term[3];
+ const float mean_dw_dr =
+ 0.5f * (wi_dx * hi_inv_dim_plus_one + wj_dx * hj_inv_dim_plus_one);
+ wi_dx_term[0] = dx[0] * r_inv * mean_dw_dr;
+ wi_dx_term[1] = dx[1] * r_inv * mean_dw_dr;
+ wi_dx_term[2] = dx[2] * r_inv * mean_dw_dr;
+
+ /* Grad-h term, dW/dh */
+ const float wi_dx_gradhterm = -0.5f *
+ (hydro_dimension * wi + (r * hi_inv) * wi_dx) *
+ hi_inv_dim_plus_one;
+
+ /* 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;
+ for (int i = 0; i < 3; i++) {
+ pi->gradient.m1_bar[i] += dx[i] * wi_term * volume_j;
+
+ pi->gradient.grad_m0_bar[i] += wi_dx_term[i] * volume_j;
+
+ pi->gradient.grad_m1_bar_gradhterm[i] += dx[i] * wi_dx_gradhterm * volume_j;
+
+ for (int j = 0; j < 3; j++) {
+ pi->gradient.grad_m1_bar[i][j] += dx[j] * wi_dx_term[i] * volume_j;
+ }
+ }
+
+ for (int j = 0; j < 3; j++) {
+ for (int k = j; k < 3; k++) {
+ int i = (j == k) ? j : 2 + j + k;
+ pi->gradient.m2_bar.elements[i] += dx[j] * dx[k] * wi_term * volume_j;
+
+ pi->gradient.grad_m2_bar_gradhterm.elements[i] +=
+ dx[j] * dx[k] * wi_dx_gradhterm * volume_j;
+
+ for (int l = 0; l < 3; l++) {
+ pi->gradient.grad_m2_bar[l].elements[i] +=
+ dx[j] * dx[k] * wi_dx_term[l] * volume_j;
+ }
+ }
+ }
+}
+
+/**
+ * @brief Finishes extra kernel parts of the gradient calculation.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_end_gradient_extra_kernel(struct part *restrict p) {
+
+ 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 */
+ const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+
+ /* Volume elements */
+ const float volume = p->mass / p->rho_evol;
+
+ /* 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 (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];
+ for (int j = 0; j < 3; j++) {
+ p->gradient.grad_m1_bar[i][j] +=
+ p->gradient.grad_m1_bar_gradhterm[j] * p->dh_norm_kernel[i];
+ }
+ for (int k = 0; k < 6; k++) {
+ p->gradient.grad_m2_bar[i].elements[k] +=
+ p->gradient.grad_m2_bar_gradhterm.elements[k] * p->dh_norm_kernel[i];
+ }
+ }
+}
+
+/**
+ * @brief Prepare a particle for the force calculation.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_prepare_force_extra_kernel(struct part *restrict p) {
+
+ if (p->is_h_max) {
+ /* Use standard kernal gradients if h=h_max */
+ /* Linear-order reproducing kernel parameters are set to "defaults" */
+ p->force.A = 1.f;
+ p->force.vac_switch = 1.f;
+ memset(p->force.B, 0.f, 3 * sizeof(float));
+ memset(p->force.grad_A, 0.f, 3 * sizeof(float));
+ memset(p->force.grad_B, 0.f, 3 * 3 * sizeof(float));
+
+ return;
+ }
+
+ /* 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;
+ }
+
+ p->gradient.grad_m2_bar[0].xx += 2.f * p->gradient.m1_bar[0];
+ p->gradient.grad_m2_bar[0].xy += p->gradient.m1_bar[1];
+ p->gradient.grad_m2_bar[0].xz += p->gradient.m1_bar[2];
+
+ p->gradient.grad_m2_bar[1].yy += 2.f * p->gradient.m1_bar[1];
+ p->gradient.grad_m2_bar[1].xy += p->gradient.m1_bar[0];
+ p->gradient.grad_m2_bar[1].yz += p->gradient.m1_bar[2];
+
+ p->gradient.grad_m2_bar[2].zz += 2.f * p->gradient.m1_bar[2];
+ 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) */
+ struct sym_matrix m2_bar_inv;
+ /* 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 (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];
+ float ABA_mult_m1_bar[3][3];
+
+ sym_matrix_multiply_by_vector(m2_bar_inv_mult_m1_bar, &m2_bar_inv,
+ p->gradient.m1_bar);
+ for (int i = 0; i < 3; i++) {
+ sym_matrix_multiply_by_vector(m2_bar_inv_mult_grad_m1_bar[i], &m2_bar_inv,
+ p->gradient.grad_m1_bar[i]);
+ sym_matrix_multiplication_ABA(
+ &m2_bar_inv_mult_grad_m2_bar_mult_m2_bar_inv[i], &m2_bar_inv,
+ &p->gradient.grad_m2_bar[i]);
+ sym_matrix_multiply_by_vector(
+ ABA_mult_m1_bar[i], &m2_bar_inv_mult_grad_m2_bar_mult_m2_bar_inv[i],
+ p->gradient.m1_bar);
+ }
+
+ /* 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 (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 (Sandnes+2025 Eqn. 23) */
+ for (int i = 0; i < 3; i++) {
+ B[i] = -m2_bar_inv_mult_m1_bar[i];
+ }
+
+ /* Calculate grad A (Sandnes+2025 Eqn. B.8) */
+ for (int i = 0; i < 3; i++) {
+ grad_A[i] = p->gradient.grad_m0_bar[i];
+
+ for (int j = 0; j < 3; j++) {
+ grad_A[i] +=
+ -2 * m2_bar_inv_mult_m1_bar[j] * p->gradient.grad_m1_bar[i][j] +
+ ABA_mult_m1_bar[i][j] * p->gradient.m1_bar[j];
+ }
+
+ grad_A[i] *= -A * A;
+ }
+
+ /* 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] = -m2_bar_inv_mult_grad_m1_bar[j][i] + ABA_mult_m1_bar[j][i];
+ }
+ }
+
+ /* Store final values */
+ p->force.A = A;
+ memcpy(p->force.B, B, 3 * sizeof(float));
+ memcpy(p->force.grad_A, grad_A, 3 * sizeof(float));
+ memcpy(p->force.grad_B, grad_B, 3 * 3 * sizeof(float));
+
+ /* Vacuum-boundary proximity switch (Sandnes+2025 Eqn. 33) */
+ p->force.vac_switch = 1.f;
+ const float hB = p->h * sqrtf(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
+ const float offset = 0.8f;
+
+ if (hB > offset) {
+ const float sigma = 0.2f;
+ p->force.vac_switch =
+ expf(-(hB - offset) * (hB - offset) / (2.f * sigma * sigma));
+ }
+}
+
+/**
+ * @brief Set gradient terms for linear-order reproducing kernel.
+ *
+ * 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;
+
+ const float hi = pi->h;
+ const float hi_inv = 1.0f / hi; /* 1/h */
+ const float hi_inv_dim = pow_dimension(hi_inv); /* 1/h^d */
+ const float hi_inv_dim_plus_one = hi_inv_dim * hi_inv; /* 1/h^(d+1) */
+
+ const float hj = pj->h;
+ const float hj_inv = 1.0f / hj; /* 1/h */
+ 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) */
+
+ const float wi_dr = hi_inv_dim_plus_one * wi_dx;
+ const float wj_dr = hj_inv_dim_plus_one * wj_dx;
+
+ if ((pi->is_h_max) || (pj->is_h_max)) {
+ /* If one or both particles have h=h_max, revert to standard kernel grads,
+ * without grad-h terms */
+ Gi[0] = wi_dr * dx[0] * r_inv;
+ Gi[1] = wi_dr * dx[1] * r_inv;
+ Gi[2] = wi_dr * dx[2] * r_inv;
+
+ Gj[0] = -wj_dr * dx[0] * r_inv;
+ Gj[1] = -wj_dr * dx[1] * r_inv;
+ Gj[2] = -wj_dr * dx[2] * r_inv;
+
+ return;
+ }
+
+ /* Mean ij kernels and gradients with grad-h terms */
+ const float wi_term = 0.5f * (wi * hi_inv_dim + wj * hj_inv_dim);
+ const 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) */
+ const float Ai = pi->force.A;
+ const float Aj = pj->force.A;
+ float grad_Ai[3], grad_Aj[3], Bi[3], Bj[3], grad_Bi[3][3], grad_Bj[3][3];
+ memcpy(grad_Ai, pi->force.grad_A, 3 * sizeof(float));
+ memcpy(grad_Aj, pj->force.grad_A, 3 * sizeof(float));
+ memcpy(Bi, pi->force.B, 3 * sizeof(float));
+ memcpy(Bj, pj->force.B, 3 * sizeof(float));
+ memcpy(grad_Bi, pi->force.grad_B, 3 * 3 * sizeof(float));
+ memcpy(grad_Bj, pj->force.grad_B, 3 * 3 * sizeof(float));
+
+ for (int i = 0; i < 3; i++) {
+
+ /* Assemble Sandnes+2025 Eqn. 29 */
+ const float mean_dw_dr =
+ 0.5f * (wi_dx * hi_inv_dim_plus_one + wj_dx * hj_inv_dim_plus_one);
+ wi_dx_term[i] = dx[i] * r_inv * mean_dw_dr;
+ wj_dx_term[i] = -wi_dx_term[i];
+
+ wi_dx_term[i] += -0.5f * (hydro_dimension * wi + (r * hi_inv) * wi_dx) *
+ hi_inv_dim_plus_one * pi->dh_norm_kernel[i];
+ wj_dx_term[i] += -0.5f * (hydro_dimension * wj + (r * hj_inv) * 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;
+
+ for (int j = 0; j < 3; j++) {
+ Gi[i] += Ai * Bi[j] * dx[j] * wi_dx_term[i] +
+ grad_Ai[i] * Bi[j] * dx[j] * wi_term +
+ Ai * grad_Bi[i][j] * dx[j] * wi_term;
+
+ Gj[i] += -(Aj * Bj[j] * dx[j] * wj_dx_term[i] +
+ grad_Aj[i] * Bj[j] * dx[j] * wj_term +
+ Aj * grad_Bj[i][j] * dx[j] * wj_term);
+ }
+ }
+
+ /* 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;
+ wj_dx_term_vac[i] = -dx[i] * r_inv * wj_dx * hj_inv_dim_plus_one;
+
+ wi_dx_term_vac[i] += -(hydro_dimension * wi + (r * hi_inv) * wi_dx) *
+ hi_inv_dim_plus_one * pi->dh_norm_kernel[i];
+ wj_dx_term_vac[i] += -(hydro_dimension * wj + (r * hj_inv) * wj_dx) *
+ hj_inv_dim_plus_one * pj->dh_norm_kernel[i];
+ }
+
+ /* 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];
+ Gj[i] = pj->force.vac_switch * Gj[i] +
+ (1.f - pj->force.vac_switch) * wj_dx_term_vac[i];
+ }
+}
+
+#endif /* SWIFT_REMIX_HYDRO_KERNELS_H */
diff --git a/src/hydro/REMIX/hydro_parameters.h b/src/hydro/REMIX/hydro_parameters.h
new file mode 100644
index 0000000000..e33568a31b
--- /dev/null
+++ b/src/hydro/REMIX/hydro_parameters.h
@@ -0,0 +1,186 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2024 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+ * 2024 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ * 2019 Josh Borrow (joshua.borrow@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_REMIX_HYDRO_PARAMETERS_H
+#define SWIFT_REMIX_HYDRO_PARAMETERS_H
+
+/* Configuration file */
+#include "config.h"
+
+/* Global headers */
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
+/* Local headers */
+#include "common_io.h"
+#include "error.h"
+#include "inline.h"
+
+/**
+ * @file REMIX/hydro_parameters.h
+ * @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
+ * as well as a number of compile-time parameters.
+ */
+
+/* Viscosity paramaters -- Defaults; can be changed at run-time */
+
+/*! Default REMIX artificial viscosity parameters */
+#define const_remix_visc_alpha 1.5f
+#define const_remix_visc_beta 3.f
+#define const_remix_visc_epsilon 0.1f
+#define const_remix_visc_a 2.0f / 3.0f
+#define const_remix_visc_b 1.0f / 3.0f
+#define const_remix_difn_a_u 0.05f
+#define const_remix_difn_b_u 0.95f
+#define const_remix_difn_a_rho 0.05f
+#define const_remix_difn_b_rho 0.95f
+#define const_remix_norm_alpha 1.0f
+#define const_remix_slope_limiter_exp_denom 0.04f
+
+/* The viscosity that the particles are reset to after being hit by a
+ * feedback event. This should be set to the same value as the
+ * const_viscosity_alpha in fixed schemes, and likely
+ * to const_viscosity_alpha_max in variable schemes. */
+#define hydro_props_default_viscosity_alpha_feedback_reset 1.5f
+
+/* Structs that store the relevant variables */
+
+/*! Artificial viscosity parameters */
+struct viscosity_global_data {};
+
+/*! Artificial diffusion parameters */
+struct diffusion_global_data {};
+
+/* Functions for reading from parameter file */
+
+/* Forward declartions */
+struct swift_params;
+struct phys_const;
+struct unit_system;
+
+/* Viscosity */
+
+/**
+ * @brief Initialises the viscosity parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct viscosity_global_data* viscosity) {}
+
+/**
+ * @brief Initialises a viscosity struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct to be filled.
+ **/
+static INLINE void viscosity_init_no_hydro(
+ struct viscosity_global_data* viscosity) {}
+
+/**
+ * @brief Prints out the viscosity parameters at the start of a run.
+ *
+ * @param viscosity: pointer to the viscosity_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void viscosity_print(
+ const struct viscosity_global_data* viscosity) {}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Prints the viscosity information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param viscosity: pointer to the viscosity_global_data struct.
+ **/
+static INLINE void viscosity_print_snapshot(
+ hid_t h_grpsph, const struct viscosity_global_data* viscosity) {
+
+ io_write_attribute_f(h_grpsph, "Viscosity alpha", const_remix_visc_alpha);
+ io_write_attribute_f(h_grpsph, "Viscosity beta", const_remix_visc_beta);
+ io_write_attribute_f(h_grpsph, "Viscosity epsilon", const_remix_visc_epsilon);
+ io_write_attribute_f(h_grpsph, "Viscosity a_visc", const_remix_visc_a);
+ io_write_attribute_f(h_grpsph, "Viscosity b_visc", const_remix_visc_b);
+}
+#endif
+
+/* Diffusion */
+
+/**
+ * @brief Initialises the diffusion parameters in the struct from
+ * the parameter file, or sets them to defaults.
+ *
+ * @param params: the pointer to the swift_params file
+ * @param unit_system: pointer to the unit system
+ * @param phys_const: pointer to the physical constants system
+ * @param diffusion_global_data: pointer to the diffusion struct to be filled.
+ **/
+static INLINE void diffusion_init(struct swift_params* params,
+ const struct unit_system* us,
+ const struct phys_const* phys_const,
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Initialises a diffusion struct to sensible numbers for mocking
+ * purposes.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct to be filled.
+ **/
+static INLINE void diffusion_init_no_hydro(
+ struct diffusion_global_data* diffusion) {}
+
+/**
+ * @brief Prints out the diffusion parameters at the start of a run.
+ *
+ * @param diffusion: pointer to the diffusion_global_data struct found in
+ * hydro_properties
+ **/
+static INLINE void diffusion_print(
+ const struct diffusion_global_data* diffusion) {}
+
+#ifdef HAVE_HDF5
+/**
+ * @brief Prints the diffusion information to the snapshot when writing.
+ *
+ * @param h_grpsph: the SPH group in the ICs to write attributes to.
+ * @param diffusion: pointer to the diffusion_global_data struct.
+ **/
+static INLINE void diffusion_print_snapshot(
+ hid_t h_grpsph, const struct diffusion_global_data* diffusion) {
+ io_write_attribute_f(h_grpsph, "Diffusion a_difn_u", const_remix_difn_a_u);
+ io_write_attribute_f(h_grpsph, "Diffusion b_difn_u", const_remix_difn_b_u);
+ io_write_attribute_f(h_grpsph, "Diffusion a_difn_rho",
+ const_remix_difn_a_rho);
+ io_write_attribute_f(h_grpsph, "Diffusion b_difn_rho",
+ const_remix_difn_b_rho);
+}
+#endif
+
+#endif /* SWIFT_REMIX_HYDRO_PARAMETERS_H */
diff --git a/src/hydro/REMIX/hydro_part.h b/src/hydro/REMIX/hydro_part.h
new file mode 100644
index 0000000000..e0a642ad0d
--- /dev/null
+++ b/src/hydro/REMIX/hydro_part.h
@@ -0,0 +1,313 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2024 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+ * 2024 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ * 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_REMIX_HYDRO_PART_H
+#define SWIFT_REMIX_HYDRO_PART_H
+
+/**
+ * @file REMIX/hydro_part.h
+ * @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 "feedback_struct.h"
+#include "mhd_struct.h"
+#include "particle_splitting_struct.h"
+#include "rt_struct.h"
+#include "sink_struct.h"
+#include "star_formation_struct.h"
+#include "symmetric_matrix.h"
+#include "timestep_limiter_struct.h"
+#include "tracers_struct.h"
+
+/**
+ * @brief Particle fields not needed during the SPH loops over neighbours.
+ *
+ * This structure contains the particle fields that are not used in the
+ * density or force loops. Quantities should be used in the kick, drift and
+ * potentially ghost tasks only.
+ */
+struct xpart {
+
+ /*! Offset between current position and position at last tree rebuild. */
+ float x_diff[3];
+
+ /*! Offset between the current position and position at the last sort. */
+ float x_diff_sort[3];
+
+ /*! Velocity at the last full step. */
+ float v_full[3];
+
+ /*! Gravitational acceleration at the end of the last step */
+ float a_grav[3];
+
+ /*! Internal energy at the last full step. */
+ float u_full;
+
+ /*! Evolved density at the last full step. */
+ float rho_evol_full;
+
+ /*! Additional data used to record particle splits */
+ struct particle_splitting_data split_data;
+
+ /*! Additional data used to record cooling information */
+ struct cooling_xpart_data cooling_data;
+
+ /*! Additional data used by the tracers */
+ struct tracers_xpart_data tracers_data;
+
+ /*! Additional data used by the star formation */
+ struct star_formation_xpart_data sf_data;
+
+ /*! Additional data used by the feedback */
+ struct feedback_part_data feedback_data;
+
+ /*! Additional data used by the MHD scheme */
+ struct mhd_xpart_data mhd_data;
+
+} SWIFT_STRUCT_ALIGN;
+
+/**
+ * @brief Particle fields for the SPH particles
+ *
+ * The density and force substructures are used to contain variables only used
+ * within the density and force loops over neighbours. All more permanent
+ * variables should be declared in the main part of the part structure,
+ */
+struct part {
+
+ /*! Particle unique ID. */
+ long long id;
+
+ /*! Pointer to corresponding gravity part. */
+ struct gpart* gpart;
+
+ /*! Particle position. */
+ double x[3];
+
+ /*! Particle predicted velocity. */
+ float v[3];
+
+ /*! Particle acceleration. */
+ float a_hydro[3];
+
+ /*! Particle mass. */
+ float mass;
+
+ /*! Particle smoothing length. */
+ float h;
+
+ /*! Particle internal energy. */
+ float u;
+
+ /*! Time derivative of the internal energy. */
+ float u_dt;
+
+ /*! Particle density (standard SPH estimate). */
+ float rho;
+
+ /*! Particle evolved density (primary density for REMIX equations). */
+ float rho_evol;
+
+ /*! Time derivative of the evolved density. */
+ float drho_dt;
+
+ /*! Gradient of velocity, calculated using normalised kernel. */
+ float dv_norm_kernel[3][3];
+
+ /*! Gradient of internal energy, calculated using normalised kernel. */
+ float du_norm_kernel[3];
+
+ /*! Gradient of density, calculated using normalised kernel. */
+ float drho_norm_kernel[3];
+
+ /*! Gradient of smoothing length, calculated using normalised kernel. */
+ float dh_norm_kernel[3];
+
+ /*! Geometric moment m_0. */
+ float m0;
+
+ /*! Gradient of geometric moment m_0. */
+ float grad_m0[3];
+
+ /* Store density/force specific stuff. */
+ union {
+
+ /**
+ * @brief Structure for the variables only used in the density loop over
+ * neighbours.
+ *
+ * Quantities in this sub-structure should only be accessed in the density
+ * loop over neighbours and the ghost task.
+ */
+ struct {
+
+ /*! Neighbour number count. */
+ float wcount;
+
+ /*! Derivative of the neighbour number with respect to h. */
+ float wcount_dh;
+
+ /*! Derivative of density with respect to h */
+ float rho_dh;
+
+ } density;
+
+ /**
+ * @brief Structure for the variables only used in the gradient loop over
+ * neighbours.
+ *
+ * Quantities should only be accessed in the gradient loop over neighbours
+ * and the extra ghost task.
+ */
+ struct {
+
+ /*! Symmetrised kernel geometric moment m0. */
+ float m0_bar;
+
+ /*! Gradient of symmetrised kernel geometric moment m0. */
+ float grad_m0_bar[3];
+
+ /*! Contributing term for grad-h component of grad_m0_bar. */
+ float grad_m0_bar_gradhterm;
+
+ /*! Symmetrised kernel geometric moment m1. */
+ float m1_bar[3];
+
+ /*! Gradient of symmetrised kernel geometric moment m1. */
+ float grad_m1_bar[3][3];
+
+ /*! Contributing term for grad-h component of grad_m1_bar. */
+ float grad_m1_bar_gradhterm[3];
+
+ /*! Symmetrised kernel geometric moment m2. */
+ struct sym_matrix m2_bar;
+
+ /*! Gradient of symmetrised kernel geometric moment m2 (x, y, z
+ * components). */
+ struct sym_matrix grad_m2_bar[3];
+
+ /*! Contributing term for grad-h component of grad_m2_bar. */
+ struct sym_matrix grad_m2_bar_gradhterm;
+
+ } gradient;
+
+ /**
+ * @brief Structure for the variables only used in the force loop over
+ * neighbours.
+ *
+ * Quantities in this sub-structure should only be accessed in the force
+ * loop over neighbours and the ghost, drift and kick tasks.
+ */
+ struct {
+
+ /*! Particle pressure. */
+ float pressure;
+
+ /*! Particle soundspeed. */
+ float soundspeed;
+
+ /*! Particle signal velocity */
+ float v_sig;
+
+ /*! Time derivative of smoothing length */
+ float h_dt;
+
+ /*! Balsara switch */
+ float balsara;
+
+ /*! Linear-order reproducing kernel coefficient */
+ float A;
+
+ /*! Linear-order reproducing kernel coefficient */
+ float B[3];
+
+ /*! Gradient of linear-order reproducing kernel coefficient */
+ float grad_A[3];
+
+ /*! Gradient of linear-order reproducing kernel coefficient */
+ float grad_B[3][3];
+
+ /*! Variable switch to identify proximity to vacuum boundaries. */
+ float vac_switch;
+
+ /*! Constant needed for slope limiter: 1 / eta_neighbours. */
+ float eta_crit;
+
+ } force;
+ };
+
+ /*! Additional data used by the MHD scheme */
+ struct mhd_part_data mhd_data;
+
+ /*! Chemistry information */
+ struct chemistry_part_data chemistry_data;
+
+ /*! Cooling information */
+ struct cooling_part_data cooling_data;
+
+ /*! Black holes information (e.g. swallowing ID) */
+ struct black_holes_part_data black_holes_data;
+
+ /*! Sink information (e.g. swallowing ID) */
+ struct sink_part_data sink_data;
+
+ /*! Material identifier flag */
+ enum eos_planetary_material_id mat_id;
+
+ /*! Additional Radiative Transfer Data */
+ struct rt_part_data rt_data;
+
+ /*! RT sub-cycling time stepping data */
+ struct rt_timestepping_data rt_time_data;
+
+ /*! Time-step length */
+ timebin_t time_bin;
+
+ /*! Tree-depth at which size / 2 <= h * gamma < size */
+ char depth_h;
+
+ /* Whether or not the particle has h=h_max ('1' or '0') */
+ char is_h_max;
+
+ /*! Time-step limiter information */
+ struct timestep_limiter_data limiter_data;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Time of the last drift */
+ integertime_t ti_drift;
+
+ /* Time of the last kick */
+ integertime_t ti_kick;
+
+#endif
+
+#ifdef PLANETARY_FIXED_ENTROPY
+ /* Fixed specific entropy */
+ float s_fixed;
+#endif
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_REMIX_HYDRO_PART_H */
diff --git a/src/hydro/REMIX/hydro_visc_difn.h b/src/hydro/REMIX/hydro_visc_difn.h
new file mode 100644
index 0000000000..8f387c3c4d
--- /dev/null
+++ b/src/hydro/REMIX/hydro_visc_difn.h
@@ -0,0 +1,452 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2024 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+ * 2024 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_REMIX_HYDRO_VISC_DIFN_H
+#define SWIFT_REMIX_HYDRO_VISC_DIFN_H
+
+/**
+ * @file REMIX/hydro_visc_difn.h
+ * @brief Utilities for REMIX artificial viscosity and diffusion calculations.
+ */
+
+#include "const.h"
+#include "hydro_parameters.h"
+#include "math.h"
+
+/**
+ * @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_visc_difn(struct part *restrict p) {
+
+ memset(p->du_norm_kernel, 0.f, 3 * sizeof(float));
+ memset(p->drho_norm_kernel, 0.f, 3 * sizeof(float));
+ memset(p->dh_norm_kernel, 0.f, 3 * sizeof(float));
+ memset(p->dv_norm_kernel, 0.f, 3 * 3 * sizeof(float));
+}
+
+/**
+ * @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_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;
+
+ const float hi = pi->h;
+ const float hi_inv = 1.0f / hi; /* 1/h */
+ const float hi_inv_dim = pow_dimension(hi_inv); /* 1/h^d */
+ const float hi_inv_dim_plus_one = hi_inv_dim * hi_inv; /* 1/h^(d+1) */
+
+ const float hj = pj->h;
+ const float hj_inv = 1.0f / hj; /* 1/h */
+ 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) */
+
+ const float volume_i = pi->mass / pi->rho_evol;
+ const float volume_j = pj->mass / pj->rho_evol;
+
+ /* 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;
+ wj_dx_term[i] = -dx[i] * r_inv * wj_dx * hj_inv_dim_plus_one;
+
+ wi_dx_term[i] *= (1.f / pi->m0);
+ wj_dx_term[i] *= (1.f / pj->m0);
+
+ wi_dx_term[i] += -wi * hi_inv_dim * pi->grad_m0[i] / (pi->m0 * pi->m0);
+ wj_dx_term[i] += -wj * hj_inv_dim * pj->grad_m0[i] / (pj->m0 * pj->m0);
+ }
+
+ /* 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
+ */
+ 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;
+
+ pi->drho_norm_kernel[i] +=
+ (pj->rho_evol - pi->rho_evol) * wi_dx_term[i] * volume_j;
+ pj->drho_norm_kernel[i] +=
+ (pi->rho_evol - pj->rho_evol) * wj_dx_term[i] * volume_i;
+ }
+
+ for (int j = 0; j < 3; j++) {
+ pi->dv_norm_kernel[i][j] +=
+ (pj->v[j] - pi->v[j]) * wi_dx_term[i] * volume_j;
+ pj->dv_norm_kernel[i][j] +=
+ (pi->v[j] - pj->v[j]) * wj_dx_term[i] * volume_i;
+ }
+ }
+}
+
+/**
+ * @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_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;
+
+ const float hi = pi->h;
+ const float hi_inv = 1.0f / hi; /* 1/h */
+ const float hi_inv_dim = pow_dimension(hi_inv); /* 1/h^d */
+ const float hi_inv_dim_plus_one = hi_inv_dim * hi_inv; /* 1/h^(d+1) */
+
+ const float volume_j = pj->mass / pj->rho_evol;
+
+ /* 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;
+ wi_dx_term[i] *= (1.f / pi->m0);
+ wi_dx_term[i] += -wi * hi_inv_dim * pi->grad_m0[i] / (pi->m0 * pi->m0);
+ }
+
+ /* 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
+ */
+ 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] +=
+ (pj->rho_evol - pi->rho_evol) * wi_dx_term[i] * volume_j;
+ }
+
+ for (int j = 0; j < 3; j++) {
+ pi->dv_norm_kernel[i][j] +=
+ (pj->v[j] - pi->v[j]) * wi_dx_term[i] * volume_j;
+ }
+ }
+}
+
+/**
+ * @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_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 r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+ const float r = sqrtf(r2);
+
+ const float hi_inv = 1.0f / pi->h;
+ const float hj_inv = 1.0f / pj->h;
+
+ /* Reconstructed velocities at the halfway point between particles
+ * (Sandnes+2025 Eqn. 36) */
+ float vtilde_i[3], vtilde_j[3];
+
+ /* 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;
+ const float a_visc = const_remix_visc_a;
+ const float b_visc = const_remix_visc_b;
+ const float eta_crit = 0.5f * (pi->force.eta_crit + pj->force.eta_crit);
+ const float slope_limiter_exp_denom = const_remix_slope_limiter_exp_denom;
+
+ if ((pi->is_h_max) || (pj->is_h_max)) {
+ /* Don't reconstruct velocity if either particle has h=h_max */
+ vtilde_i[0] = 0.f;
+ vtilde_i[1] = 0.f;
+ vtilde_i[2] = 0.f;
+
+ vtilde_j[0] = 0.f;
+ vtilde_j[1] = 0.f;
+ vtilde_j[2] = 0.f;
+
+ } else {
+ /* A numerators and denominators (Sandnes+2025 Eqn. 38) */
+ float A_i_v = 0.f;
+ float A_j_v = 0.f;
+
+ /* 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 (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];
+
+ 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;
+ phi_j_v = 1.f;
+
+ } else if ((A_i_v == 0.f && A_j_v != 0.f) ||
+ (A_j_v == 0.f && A_i_v != 0.f) || (A_i_v == -A_j_v)) {
+ phi_i_v = 0.f;
+ phi_j_v = 0.f;
+ } else {
+ /* 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);
+
+ phi_j_v = min(1.f, 4.f * A_j_v / A_i_v / (1.f + A_j_v / A_i_v) /
+ (1.f + A_j_v / A_i_v));
+ phi_j_v = max(0.f, phi_j_v);
+ }
+
+ /* exp in slope limiter (middle case in Sandnes+2025 Eqn. 37) */
+ const 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) /
+ slope_limiter_exp_denom);
+ phi_j_v *= expf(-(eta_ab - eta_crit) * (eta_ab - eta_crit) /
+ slope_limiter_exp_denom);
+ }
+
+ for (int i = 0; i < 3; i++) {
+ /* 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] =
+ pj->v[i] + (1.f - pj->force.balsara) * phi_j_v * v_reconst_j[i];
+ }
+ }
+
+ /* Assemble Sandnes+2025 Eqn. 40 */
+ const float mu_i =
+ min(0.f, ((vtilde_i[0] - vtilde_j[0]) * dx[0] +
+ (vtilde_i[1] - vtilde_j[1]) * dx[1] +
+ (vtilde_i[2] - vtilde_j[2]) * dx[2]) *
+ hi_inv / (r2 * hi_inv * hi_inv + epsilon * epsilon));
+ const float mu_j =
+ min(0.f, ((vtilde_i[0] - vtilde_j[0]) * dx[0] +
+ (vtilde_i[1] - vtilde_j[1]) * dx[1] +
+ (vtilde_i[2] - vtilde_j[2]) * dx[2]) *
+ hj_inv / (r2 * hj_inv * hj_inv + epsilon * epsilon));
+
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
+
+ /* 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 */
+ 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 r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+ const float r = sqrtf(r2);
+
+ const float hi_inv = 1.0f / pi->h;
+ const float hj_inv = 1.0f / pj->h;
+
+ const float eta_crit = 0.5f * (pi->force.eta_crit + pj->force.eta_crit);
+ const float slope_limiter_exp_denom = const_remix_slope_limiter_exp_denom;
+
+ if ((pi->is_h_max) || (pj->is_h_max)) {
+ /* Don't reconstruct internal energy of density if either particle has
+ * h=h_max */
+ *utilde_i = 0.f;
+ *utilde_j = 0.f;
+ *rhotilde_i = 0.f;
+ *rhotilde_j = 0.f;
+
+ return;
+ }
+
+ /* 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;
+
+ /* 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 (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];
+
+ 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];
+ rho_reconst_j += 0.5 * pj->drho_norm_kernel[i] * dx[i];
+ }
+
+ 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;
+
+ } else if ((A_i_u == 0.f && A_j_u != 0.f) || (A_j_u == 0.f && A_i_u != 0.f) ||
+ (A_i_u == -A_j_u)) {
+ phi_i_u = 0.f;
+ phi_j_u = 0.f;
+ } else {
+ /* 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);
+
+ phi_j_u = min(1.f, 4.f * A_j_u / A_i_u / (1.f + A_j_u / A_i_u) /
+ (1.f + A_j_u / A_i_u));
+ 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;
+
+ } else if ((A_i_rho == 0.f && A_j_rho != 0.f) ||
+ (A_j_rho == 0.f && A_i_rho != 0.f) || (A_i_rho == -A_j_rho)) {
+ phi_i_rho = 0.f;
+ phi_j_rho = 0.f;
+ } else {
+ /* 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);
+
+ phi_j_rho = min(1.f, 4.f * A_j_rho / A_i_rho / (1.f + A_j_rho / A_i_rho) /
+ (1.f + A_j_rho / A_i_rho));
+ phi_j_rho = max(0.f, phi_j_rho);
+ }
+
+ /* exp in slope limiter (middle case in Sandnes+2025 Eqn. 37) */
+ const 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);
+ phi_j_u *= expf(-(eta_ab - eta_crit) * (eta_ab - eta_crit) /
+ slope_limiter_exp_denom);
+ phi_i_rho *= expf(-(eta_ab - eta_crit) * (eta_ab - eta_crit) /
+ slope_limiter_exp_denom);
+ phi_j_rho *= expf(-(eta_ab - eta_crit) * (eta_ab - eta_crit) /
+ slope_limiter_exp_denom);
+ }
+
+ /* 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;
+}
+
+#endif /* SWIFT_REMIX_HYDRO_VISC_DIFN_H */
diff --git a/src/hydro_csds.h b/src/hydro_csds.h
index 96cb374423..a0593ca7b3 100644
--- a/src/hydro_csds.h
+++ b/src/hydro_csds.h
@@ -47,6 +47,8 @@
#error TODO
#elif defined(PLANETARY_SPH)
#error TODO
+#elif defined(REMIX_SPH)
+#error TODO
#elif defined(SPHENIX_SPH)
#include "./hydro/SPHENIX/hydro_csds.h"
#elif defined(GASOLINE_SPH)
diff --git a/src/hydro_io.h b/src/hydro_io.h
index ed78acf691..5a64a284cc 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -43,6 +43,8 @@
#include "./hydro/Shadowswift/hydro_io.h"
#elif defined(PLANETARY_SPH)
#include "./hydro/Planetary/hydro_io.h"
+#elif defined(REMIX_SPH)
+#include "./hydro/REMIX/hydro_io.h"
#elif defined(SPHENIX_SPH)
#include "./hydro/SPHENIX/hydro_io.h"
#elif defined(GASOLINE_SPH)
diff --git a/src/hydro_parameters.h b/src/hydro_parameters.h
index 645a5fe5af..46d93ad43a 100644
--- a/src/hydro_parameters.h
+++ b/src/hydro_parameters.h
@@ -52,6 +52,8 @@
#include "./hydro/Shadowswift/hydro_parameters.h"
#elif defined(PLANETARY_SPH)
#include "./hydro/Planetary/hydro_parameters.h"
+#elif defined(REMIX_SPH)
+#include "./hydro/REMIX/hydro_parameters.h"
#elif defined(SPHENIX_SPH)
#include "./hydro/SPHENIX/hydro_parameters.h"
#elif defined(GASOLINE_SPH)
diff --git a/src/part.h b/src/part.h
index bb574651b5..c68c91cfa8 100644
--- a/src/part.h
+++ b/src/part.h
@@ -80,6 +80,11 @@ struct threadpool;
#elif defined(PLANETARY_SPH)
#include "./hydro/Planetary/hydro_part.h"
#define hydro_need_extra_init_loop 0
+#elif defined(REMIX_SPH)
+#include "./hydro/REMIX/hydro_part.h"
+#define hydro_need_extra_init_loop 0
+#define EXTRA_HYDRO_LOOP
+#define EXTRA_HYDRO_LOOP_TYPE2
#elif defined(SPHENIX_SPH)
#include "./hydro/SPHENIX/hydro_part.h"
#define hydro_need_extra_init_loop 0
diff --git a/src/runner_main.c b/src/runner_main.c
index 8b2d268144..e6b3b1ef05 100644
--- a/src/runner_main.c
+++ b/src/runner_main.c
@@ -214,8 +214,13 @@ void *runner_main(void *data) {
/*limit_h_max=*/0);
#ifdef EXTRA_HYDRO_LOOP
else if (t->subtype == task_subtype_gradient)
+#ifdef EXTRA_HYDRO_LOOP_TYPE2
+ runner_doself2_branch_gradient(r, ci, /*limit_h_min=*/0,
+ /*limit_h_max=*/0);
+#else
runner_doself1_branch_gradient(r, ci, /*limit_h_min=*/0,
/*limit_h_max=*/0);
+#endif
#endif
else if (t->subtype == task_subtype_force)
runner_doself2_branch_force(r, ci, /*limit_h_min=*/0,
@@ -274,8 +279,13 @@ void *runner_main(void *data) {
/*limit_h_max=*/0);
#ifdef EXTRA_HYDRO_LOOP
else if (t->subtype == task_subtype_gradient)
+#ifdef EXTRA_HYDRO_LOOP_TYPE2
+ runner_dopair2_branch_gradient(r, ci, cj, /*limit_h_min=*/0,
+ /*limit_h_max=*/0);
+#else
runner_dopair1_branch_gradient(r, ci, cj, /*limit_h_min=*/0,
/*limit_h_max=*/0);
+#endif
#endif
else if (t->subtype == task_subtype_force)
runner_dopair2_branch_force(r, ci, cj, /*limit_h_min=*/0,
@@ -331,7 +341,11 @@ void *runner_main(void *data) {
runner_dosub_self1_density(r, ci, /*below_h_max=*/0, 1);
#ifdef EXTRA_HYDRO_LOOP
else if (t->subtype == task_subtype_gradient)
+#ifdef EXTRA_HYDRO_LOOP_TYPE2
+ runner_dosub_self2_gradient(r, ci, /*below_h_max=*/0, 1);
+#else
runner_dosub_self1_gradient(r, ci, /*below_h_max=*/0, 1);
+#endif
#endif
else if (t->subtype == task_subtype_force)
runner_dosub_self2_force(r, ci, /*below_h_max=*/0, 1);
@@ -379,7 +393,11 @@ void *runner_main(void *data) {
runner_dosub_pair1_density(r, ci, cj, /*below_h_max=*/0, 1);
#ifdef EXTRA_HYDRO_LOOP
else if (t->subtype == task_subtype_gradient)
+#ifdef EXTRA_HYDRO_LOOP_TYPE2
+ runner_dosub_pair2_gradient(r, ci, cj, /*below_h_max=*/0, 1);
+#else
runner_dosub_pair1_gradient(r, ci, cj, /*below_h_max=*/0, 1);
+#endif
#endif
else if (t->subtype == task_subtype_force)
runner_dosub_pair2_force(r, ci, cj, /*below_h_max=*/0, 1);
diff --git a/src/symmetric_matrix.h b/src/symmetric_matrix.h
new file mode 100644
index 0000000000..88c82df542
--- /dev/null
+++ b/src/symmetric_matrix.h
@@ -0,0 +1,167 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2023 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_SYMMETRIC_MATRIX_H
+#define SWIFT_SYMMETRIC_MATRIX_H
+
+/* Local includes */
+#include "dimension.h"
+#include "error.h"
+
+/**
+ * @brief Symmetric matrix definition in 3D.
+ *
+ * The matrix elements can be accessed as an array or via their "coordinates".
+ * Replicated elements are not stored.
+ */
+struct sym_matrix {
+
+ union {
+ struct {
+ float elements[6];
+ };
+ struct {
+ float xx;
+ float yy;
+ float zz;
+ float xy;
+ float xz;
+ float yz;
+ };
+ };
+};
+
+/**
+ * @brief Zero the matrix
+ */
+__attribute__((always_inline)) INLINE static void zero_sym_matrix(
+ struct sym_matrix *M) {
+ for (int i = 0; i < 6; ++i) M->elements[i] = 0.f;
+}
+
+/**
+ * @brief Construct a 3x3 array from a symmetric matrix.
+ */
+__attribute__((always_inline)) INLINE static void get_matrix_from_sym_matrix(
+ float out[3][3], const struct sym_matrix *in) {
+
+ out[0][0] = in->xx;
+ out[0][1] = in->xy;
+ out[0][2] = in->xz;
+ out[1][0] = in->xy;
+ out[1][1] = in->yy;
+ out[1][2] = in->yz;
+ out[2][0] = in->xz;
+ out[2][1] = in->yz;
+ out[2][2] = in->zz;
+}
+
+/**
+ * @brief Construct a symmetric matrix from a 3x3 array.
+ *
+ * No check is performed to verify the input 3x3 array is indeed symmetric.
+ */
+__attribute__((always_inline)) INLINE static void get_sym_matrix_from_matrix(
+ struct sym_matrix *out, const float in[3][3]) {
+ out->xx = in[0][0];
+ out->yy = in[1][1];
+ out->zz = in[2][2];
+ out->xy = in[0][1];
+ out->xz = in[0][2];
+ out->yz = in[1][2];
+}
+
+/**
+ * @brief Compute the product of a symmetric matrix and a vector.
+ */
+__attribute__((always_inline)) INLINE static void sym_matrix_multiply_by_vector(
+ float out[3], const struct sym_matrix *M, const float v[3]) {
+
+ out[0] = M->xx * v[0] + M->xy * v[1] + M->xz * v[2];
+ out[1] = M->xy * v[0] + M->yy * v[1] + M->yz * v[2];
+ out[2] = M->xz * v[0] + M->yz * v[1] + M->zz * v[2];
+}
+
+/**
+ * @brief Multiply two symmetric matrices in two operations, ABA.
+ */
+__attribute__((always_inline)) INLINE static void sym_matrix_multiplication_ABA(
+ struct sym_matrix *M_out, const struct sym_matrix *A,
+ const struct sym_matrix *B) {
+
+ float BA_array[3][3] = {0};
+ float A_array[3][3], B_array[3][3];
+ get_matrix_from_sym_matrix(A_array, A);
+ get_matrix_from_sym_matrix(B_array, B);
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ for (int k = 0; k < 3; k++) {
+ BA_array[i][j] += B_array[i][k] * A_array[k][j];
+ }
+ }
+ }
+
+ M_out->xx = (A->xx * BA_array[0][0] + A->xy * BA_array[1][0] +
+ A->xz * BA_array[2][0]);
+ M_out->yy = (A->xy * BA_array[0][1] + A->yy * BA_array[1][1] +
+ A->yz * BA_array[2][1]);
+ M_out->zz = (A->xz * BA_array[0][2] + A->yz * BA_array[1][2] +
+ A->zz * BA_array[2][2]);
+ M_out->xy = (A->xy * BA_array[0][0] + A->yy * BA_array[1][0] +
+ A->yz * BA_array[2][0]);
+ M_out->xz = (A->xz * BA_array[0][0] + A->yz * BA_array[1][0] +
+ A->zz * BA_array[2][0]);
+ M_out->yz = (A->xz * BA_array[0][1] + A->yz * BA_array[1][1] +
+ A->zz * BA_array[2][1]);
+}
+
+/**
+ * @brief Print a symmetric matrix.
+ */
+__attribute__((always_inline)) INLINE static void sym_matrix_print(
+ const struct sym_matrix *M) {
+ message("|%.4f %.4f %.4f|", M->xx, M->xy, M->xz);
+ message("|%.4f %.4f %.4f|", M->xy, M->yy, M->yz);
+ message("|%.4f %.4f %.4f|", M->xz, M->yz, M->zz);
+}
+
+/**
+ * @brief Compute the inverse of a symmetric matrix and a vector.
+ *
+ * Returned as a symmetric matrix
+ */
+__attribute__((always_inline)) INLINE static void sym_matrix_invert(
+ struct sym_matrix *M_inv, const struct sym_matrix *M) {
+
+ float M_inv_matrix[3][3];
+ get_matrix_from_sym_matrix(M_inv_matrix, M);
+ int res = invert_dimension_by_dimension_matrix(M_inv_matrix);
+ if (res) {
+ sym_matrix_print(M);
+ error("Error inverting matrix");
+ }
+
+ M_inv->xx = M_inv_matrix[0][0];
+ M_inv->yy = M_inv_matrix[1][1];
+ M_inv->zz = M_inv_matrix[2][2];
+ M_inv->xy = M_inv_matrix[0][1];
+ M_inv->xz = M_inv_matrix[0][2];
+ M_inv->yz = M_inv_matrix[1][2];
+}
+
+#endif /* SWIFT_SYMMETRIC_MATRIX_H */
diff --git a/src/tools.c b/src/tools.c
index 65d7d976cc..44c8fbdc39 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -55,6 +55,7 @@
#include "runner.h"
#include "sink_iact.h"
#include "sink_properties.h"
+#include "space_getsid.h"
#include "star_formation_iact.h"
#include "stars.h"
@@ -205,13 +206,18 @@ void pairs_single_density(double *dim, long long int pid,
void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
- float r2, hi, hj, hig2, hjg2, dx[3];
+ float hi, hj, hig2, hjg2;
struct part *pi, *pj;
const double dim[3] = {r->e->s->dim[0], r->e->s->dim[1], r->e->s->dim[2]};
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
const float a = cosmo->a;
const float H = cosmo->H;
+ double shift[3] = {0.0, 0.0, 0.0};
+ space_getsid_and_swap_cells(e->s, &ci, &cj, shift);
+ const double shift_i[3] = {cj->loc[0] + shift[0], cj->loc[1] + shift[1],
+ cj->loc[2] + shift[2]};
+ const double shift_j[3] = {cj->loc[0], cj->loc[1], cj->loc[2]};
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->hydro.count; ++i) {
@@ -223,17 +229,23 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
/* Skip inactive particles. */
if (!part_is_active(pi, e)) continue;
+ const float pix = pi->x[0] - shift_i[0];
+ const float piy = pi->x[1] - shift_i[1];
+ const float piz = pi->x[2] - shift_i[2];
+
for (int j = 0; j < cj->hydro.count; ++j) {
pj = &cj->hydro.parts[j];
+ const float pjx = pj->x[0] - shift_j[0];
+ const float pjy = pj->x[1] - shift_j[1];
+ const float pjz = pj->x[2] - shift_j[2];
+
/* Pairwise distance */
- r2 = 0.0f;
- for (int k = 0; k < 3; k++) {
- dx[k] = ci->hydro.parts[i].x[k] - cj->hydro.parts[j].x[k];
- dx[k] = nearest(dx[k], dim[k]);
- r2 += dx[k] * dx[k];
- }
+ const float dx[3] = {nearest(pix - pjx, dim[0]),
+ nearest(piy - pjy, dim[1]),
+ nearest(piz - pjz, dim[2])};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Hit or miss? */
if (r2 < hig2 && !part_is_inhibited(pj, e)) {
@@ -258,17 +270,23 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
/* Skip inactive particles. */
if (!part_is_active(pj, e)) continue;
+ const float pjx = pj->x[0] - shift_j[0];
+ const float pjy = pj->x[1] - shift_j[1];
+ const float pjz = pj->x[2] - shift_j[2];
+
for (int i = 0; i < ci->hydro.count; ++i) {
pi = &ci->hydro.parts[i];
+ const float pix = pi->x[0] - shift_i[0];
+ const float piy = pi->x[1] - shift_i[1];
+ const float piz = pi->x[2] - shift_i[2];
+
/* Pairwise distance */
- r2 = 0.0f;
- for (int k = 0; k < 3; k++) {
- dx[k] = cj->hydro.parts[j].x[k] - ci->hydro.parts[i].x[k];
- dx[k] = nearest(dx[k], dim[k]);
- r2 += dx[k] * dx[k];
- }
+ const float dx[3] = {nearest(pjx - pix, dim[0]),
+ nearest(pjy - piy, dim[1]),
+ nearest(pjz - piz, dim[2])};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Hit or miss? */
if (r2 < hjg2 && !part_is_inhibited(pi, e)) {
@@ -287,13 +305,18 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
#ifdef EXTRA_HYDRO_LOOP
void pairs_all_gradient(struct runner *r, struct cell *ci, struct cell *cj) {
- float r2, hi, hj, hig2, hjg2, dx[3];
+ float hi, hj, hig2, hjg2;
struct part *pi, *pj;
const double dim[3] = {r->e->s->dim[0], r->e->s->dim[1], r->e->s->dim[2]};
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
const float a = cosmo->a;
const float H = cosmo->H;
+ double shift[3] = {0.0, 0.0, 0.0};
+ space_getsid_and_swap_cells(e->s, &ci, &cj, shift);
+ const double shift_i[3] = {cj->loc[0] + shift[0], cj->loc[1] + shift[1],
+ cj->loc[2] + shift[2]};
+ const double shift_j[3] = {cj->loc[0], cj->loc[1], cj->loc[2]};
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->hydro.count; ++i) {
@@ -305,25 +328,39 @@ void pairs_all_gradient(struct runner *r, struct cell *ci, struct cell *cj) {
/* Skip inactive particles. */
if (!part_is_active(pi, e)) continue;
+ const float pix = pi->x[0] - shift_i[0];
+ const float piy = pi->x[1] - shift_i[1];
+ const float piz = pi->x[2] - shift_i[2];
+
for (int j = 0; j < cj->hydro.count; ++j) {
pj = &cj->hydro.parts[j];
hj = pj->h;
hjg2 = hj * hj * kernel_gamma2;
+ const float pjx = pj->x[0] - shift_j[0];
+ const float pjy = pj->x[1] - shift_j[1];
+ const float pjz = pj->x[2] - shift_j[2];
+
/* Pairwise distance */
- r2 = 0.0f;
- for (int k = 0; k < 3; k++) {
- dx[k] = ci->hydro.parts[i].x[k] - cj->hydro.parts[j].x[k];
- dx[k] = nearest(dx[k], dim[k]);
- r2 += dx[k] * dx[k];
- }
+ const float dx[3] = {nearest(pix - pjx, dim[0]),
+ nearest(piy - pjy, dim[1]),
+ nearest(piz - pjz, dim[2])};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Hit or miss? */
- if (r2 < hig2 && !part_is_inhibited(pj, e)) {
-
- /* Interact */
- runner_iact_nonsym_gradient(r2, dx, hi, hj, pi, pj, a, H);
+ if (!part_is_inhibited(pj, e)) {
+#ifdef EXTRA_HYDRO_LOOP_TYPE2
+ if (r2 < hig2 || r2 < hjg2) {
+ /* Interact */
+ runner_iact_nonsym_gradient(r2, dx, hi, hj, pi, pj, a, H);
+ }
+#else
+ if (r2 < hig2) {
+ /* Interact */
+ runner_iact_nonsym_gradient(r2, dx, hi, hj, pi, pj, a, H);
+ }
+#endif /* EXTRA_HYDRO_LOOP_TYPE2 */
}
}
}
@@ -338,25 +375,39 @@ void pairs_all_gradient(struct runner *r, struct cell *ci, struct cell *cj) {
/* Skip inactive particles. */
if (!part_is_active(pj, e)) continue;
+ const float pjx = pj->x[0] - shift_j[0];
+ const float pjy = pj->x[1] - shift_j[1];
+ const float pjz = pj->x[2] - shift_j[2];
+
for (int i = 0; i < ci->hydro.count; ++i) {
pi = &ci->hydro.parts[i];
hi = pi->h;
hig2 = hi * hi * kernel_gamma2;
+ const float pix = pi->x[0] - shift_i[0];
+ const float piy = pi->x[1] - shift_i[1];
+ const float piz = pi->x[2] - shift_i[2];
+
/* Pairwise distance */
- r2 = 0.0f;
- for (int k = 0; k < 3; k++) {
- dx[k] = cj->hydro.parts[j].x[k] - ci->hydro.parts[i].x[k];
- dx[k] = nearest(dx[k], dim[k]);
- r2 += dx[k] * dx[k];
- }
+ const float dx[3] = {nearest(pjx - pix, dim[0]),
+ nearest(pjy - piy, dim[1]),
+ nearest(pjz - piz, dim[2])};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Hit or miss? */
- if (r2 < hjg2 && !part_is_inhibited(pi, e)) {
-
- /* Interact */
- runner_iact_nonsym_gradient(r2, dx, hj, pi->h, pj, pi, a, H);
+ if (!part_is_inhibited(pi, e)) {
+#ifdef EXTRA_HYDRO_LOOP_TYPE2
+ if (r2 < hig2 || r2 < hjg2) {
+ /* Interact */
+ runner_iact_nonsym_gradient(r2, dx, hj, pi->h, pj, pi, a, H);
+ }
+#else
+ if (r2 < hjg2) {
+ /* Interact */
+ runner_iact_nonsym_gradient(r2, dx, hj, pi->h, pj, pi, a, H);
+ }
+#endif /* EXTRA_HYDRO_LOOP_TYPE2 */
}
}
}
@@ -365,13 +416,18 @@ void pairs_all_gradient(struct runner *r, struct cell *ci, struct cell *cj) {
void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
- float r2, hi, hj, hig2, hjg2, dx[3];
+ float hi, hj, hig2, hjg2;
struct part *pi, *pj;
const double dim[3] = {r->e->s->dim[0], r->e->s->dim[1], r->e->s->dim[2]};
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
const float a = cosmo->a;
const float H = cosmo->H;
+ double shift[3] = {0.0, 0.0, 0.0};
+ space_getsid_and_swap_cells(e->s, &ci, &cj, shift);
+ const double shift_i[3] = {cj->loc[0] + shift[0], cj->loc[1] + shift[1],
+ cj->loc[2] + shift[2]};
+ const double shift_j[3] = {cj->loc[0], cj->loc[1], cj->loc[2]};
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->hydro.count; ++i) {
@@ -383,19 +439,25 @@ void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
/* Skip inactive particles. */
if (!part_is_active(pi, e)) continue;
+ const float pix = pi->x[0] - shift_i[0];
+ const float piy = pi->x[1] - shift_i[1];
+ const float piz = pi->x[2] - shift_i[2];
+
for (int j = 0; j < cj->hydro.count; ++j) {
pj = &cj->hydro.parts[j];
hj = pj->h;
hjg2 = hj * hj * kernel_gamma2;
+ const float pjx = pj->x[0] - shift_j[0];
+ const float pjy = pj->x[1] - shift_j[1];
+ const float pjz = pj->x[2] - shift_j[2];
+
/* Pairwise distance */
- r2 = 0.0f;
- for (int k = 0; k < 3; k++) {
- dx[k] = ci->hydro.parts[i].x[k] - cj->hydro.parts[j].x[k];
- dx[k] = nearest(dx[k], dim[k]);
- r2 += dx[k] * dx[k];
- }
+ const float dx[3] = {nearest(pix - pjx, dim[0]),
+ nearest(piy - pjy, dim[1]),
+ nearest(piz - pjz, dim[2])};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Hit or miss? */
if (r2 < hig2 || r2 < hjg2) {
@@ -416,19 +478,25 @@ void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
/* Skip inactive particles. */
if (!part_is_active(pj, e)) continue;
+ const float pjx = pj->x[0] - shift_j[0];
+ const float pjy = pj->x[1] - shift_j[1];
+ const float pjz = pj->x[2] - shift_j[2];
+
for (int i = 0; i < ci->hydro.count; ++i) {
pi = &ci->hydro.parts[i];
hi = pi->h;
hig2 = hi * hi * kernel_gamma2;
+ const float pix = pi->x[0] - shift_i[0];
+ const float piy = pi->x[1] - shift_i[1];
+ const float piz = pi->x[2] - shift_i[2];
+
/* Pairwise distance */
- r2 = 0.0f;
- for (int k = 0; k < 3; k++) {
- dx[k] = cj->hydro.parts[j].x[k] - ci->hydro.parts[i].x[k];
- dx[k] = nearest(dx[k], dim[k]);
- r2 += dx[k] * dx[k];
- }
+ const float dx[3] = {nearest(pjx - pix, dim[0]),
+ nearest(pjy - piy, dim[1]),
+ nearest(pjz - piz, dim[2])};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Hit or miss? */
if (r2 < hjg2 || r2 < hig2) {
@@ -520,7 +588,7 @@ void pairs_all_stars_density(struct runner *r, struct cell *ci,
}
void self_all_density(struct runner *r, struct cell *ci) {
- float r2, hi, hj, hig2, hjg2, dxi[3]; //, dxj[3];
+ float hi, hj, hig2, hjg2;
struct part *pi, *pj;
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
@@ -534,6 +602,8 @@ void self_all_density(struct runner *r, struct cell *ci) {
hi = pi->h;
hig2 = hi * hi * kernel_gamma2;
+ const double pix[3] = {pi->x[0], pi->x[1], pi->x[2]};
+
for (int j = i + 1; j < ci->hydro.count; ++j) {
pj = &ci->hydro.parts[j];
@@ -542,37 +612,37 @@ void self_all_density(struct runner *r, struct cell *ci) {
if (pi == pj) continue;
+ const double pjx[3] = {pj->x[0], pj->x[1], pj->x[2]};
+
/* Pairwise distance */
- r2 = 0.0f;
- for (int k = 0; k < 3; k++) {
- dxi[k] = ci->hydro.parts[i].x[k] - ci->hydro.parts[j].x[k];
- r2 += dxi[k] * dxi[k];
- }
+ float dx[3] = {(float)(pix[0] - pjx[0]), (float)(pix[1] - pjx[1]),
+ (float)(pix[2] - pjx[2])};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Hit or miss? */
if (r2 < hig2 && part_is_active(pi, e) && !part_is_inhibited(pj, e)) {
/* Interact */
- runner_iact_nonsym_density(r2, dxi, hi, hj, pi, pj, a, H);
- runner_iact_nonsym_chemistry(r2, dxi, hi, hj, pi, pj, a, H);
- runner_iact_nonsym_pressure_floor(r2, dxi, hi, hj, pi, pj, a, H);
- runner_iact_nonsym_star_formation(r2, dxi, hi, hj, pi, pj, a, H);
- runner_iact_nonsym_sink(r2, dxi, hi, hj, pi, pj, a, H);
+ runner_iact_nonsym_density(r2, dx, hi, hj, pi, pj, a, H);
+ runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, a, H);
+ runner_iact_nonsym_pressure_floor(r2, dx, hi, hj, pi, pj, a, H);
+ runner_iact_nonsym_star_formation(r2, dx, hi, hj, pi, pj, a, H);
+ runner_iact_nonsym_sink(r2, dx, hi, hj, pi, pj, a, H);
}
/* Hit or miss? */
if (r2 < hjg2 && part_is_active(pj, e) && !part_is_inhibited(pi, e)) {
- dxi[0] = -dxi[0];
- dxi[1] = -dxi[1];
- dxi[2] = -dxi[2];
+ dx[0] = -dx[0];
+ dx[1] = -dx[1];
+ dx[2] = -dx[2];
/* Interact */
- runner_iact_nonsym_density(r2, dxi, hj, hi, pj, pi, a, H);
- runner_iact_nonsym_chemistry(r2, dxi, hj, hi, pj, pi, a, H);
- runner_iact_nonsym_pressure_floor(r2, dxi, hj, hi, pj, pi, a, H);
- runner_iact_nonsym_star_formation(r2, dxi, hj, hi, pj, pi, a, H);
- runner_iact_nonsym_sink(r2, dxi, hj, hi, pj, pi, a, H);
+ runner_iact_nonsym_density(r2, dx, hj, hi, pj, pi, a, H);
+ runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, a, H);
+ runner_iact_nonsym_pressure_floor(r2, dx, hj, hi, pj, pi, a, H);
+ runner_iact_nonsym_star_formation(r2, dx, hj, hi, pj, pi, a, H);
+ runner_iact_nonsym_sink(r2, dx, hj, hi, pj, pi, a, H);
}
}
}
@@ -580,7 +650,7 @@ void self_all_density(struct runner *r, struct cell *ci) {
#ifdef EXTRA_HYDRO_LOOP
void self_all_gradient(struct runner *r, struct cell *ci) {
- float r2, hi, hj, hig2, hjg2, dxi[3]; //, dxj[3];
+ float hi, hj, hig2, hjg2;
struct part *pi, *pj;
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
@@ -594,6 +664,8 @@ void self_all_gradient(struct runner *r, struct cell *ci) {
hi = pi->h;
hig2 = hi * hi * kernel_gamma2;
+ const double pix[3] = {pi->x[0], pi->x[1], pi->x[2]};
+
for (int j = i + 1; j < ci->hydro.count; ++j) {
pj = &ci->hydro.parts[j];
@@ -602,29 +674,46 @@ void self_all_gradient(struct runner *r, struct cell *ci) {
if (pi == pj) continue;
+ const double pjx[3] = {pj->x[0], pj->x[1], pj->x[2]};
+
/* Pairwise distance */
- r2 = 0.0f;
- for (int k = 0; k < 3; k++) {
- dxi[k] = ci->hydro.parts[i].x[k] - ci->hydro.parts[j].x[k];
- r2 += dxi[k] * dxi[k];
- }
+ float dx[3] = {(float)(pix[0] - pjx[0]), (float)(pix[1] - pjx[1]),
+ (float)(pix[2] - pjx[2])};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Hit or miss? */
- if (r2 < hig2 && part_is_active(pi, e) && !part_is_inhibited(pj, e)) {
-
- /* Interact */
- runner_iact_nonsym_gradient(r2, dxi, hi, hj, pi, pj, a, H);
+ if (part_is_active(pi, e) && !part_is_inhibited(pj, e)) {
+#ifdef EXTRA_HYDRO_LOOP_TYPE2
+ if (r2 < hig2 || r2 < hjg2) {
+ /* Interact */
+ runner_iact_nonsym_gradient(r2, dx, hi, hj, pi, pj, a, H);
+ }
+#else
+ if (r2 < hig2) {
+ /* Interact */
+ runner_iact_nonsym_gradient(r2, dx, hi, hj, pi, pj, a, H);
+ }
+#endif /* EXTRA_HYDRO_LOOP_TYPE2 */
}
/* Hit or miss? */
- if (r2 < hjg2 && part_is_active(pj, e) && !part_is_inhibited(pi, e)) {
+ if (part_is_active(pj, e) && !part_is_inhibited(pi, e)) {
- dxi[0] = -dxi[0];
- dxi[1] = -dxi[1];
- dxi[2] = -dxi[2];
+ dx[0] = -dx[0];
+ dx[1] = -dx[1];
+ dx[2] = -dx[2];
- /* Interact */
- runner_iact_nonsym_gradient(r2, dxi, hj, hi, pj, pi, a, H);
+#ifdef EXTRA_HYDRO_LOOP_TYPE2
+ if (r2 < hig2 || r2 < hjg2) {
+ /* Interact */
+ runner_iact_nonsym_gradient(r2, dx, hj, hi, pj, pi, a, H);
+ }
+#else
+ if (r2 < hjg2) {
+ /* Interact */
+ runner_iact_nonsym_gradient(r2, dx, hj, hi, pj, pi, a, H);
+ }
+#endif /* EXTRA_HYDRO_LOOP_TYPE2 */
}
}
}
@@ -632,7 +721,7 @@ void self_all_gradient(struct runner *r, struct cell *ci) {
#endif /* EXTRA_HYDRO_LOOP */
void self_all_force(struct runner *r, struct cell *ci) {
- float r2, hi, hj, hig2, hjg2, dxi[3]; //, dxj[3];
+ float hi, hj, hig2, hjg2;
struct part *pi, *pj;
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
@@ -646,6 +735,8 @@ void self_all_force(struct runner *r, struct cell *ci) {
hi = pi->h;
hig2 = hi * hi * kernel_gamma2;
+ const double pix[3] = {pi->x[0], pi->x[1], pi->x[2]};
+
for (int j = i + 1; j < ci->hydro.count; ++j) {
pj = &ci->hydro.parts[j];
@@ -654,18 +745,18 @@ void self_all_force(struct runner *r, struct cell *ci) {
if (pi == pj) continue;
+ const double pjx[3] = {pj->x[0], pj->x[1], pj->x[2]};
+
/* Pairwise distance */
- r2 = 0.0f;
- for (int k = 0; k < 3; k++) {
- dxi[k] = ci->hydro.parts[i].x[k] - ci->hydro.parts[j].x[k];
- r2 += dxi[k] * dxi[k];
- }
+ float dx[3] = {(float)(pix[0] - pjx[0]), (float)(pix[1] - pjx[1]),
+ (float)(pix[2] - pjx[2])};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Hit or miss? */
if (r2 < hig2 || r2 < hjg2) {
/* Interact */
- runner_iact_force(r2, dxi, hi, hj, pi, pj, a, H);
+ runner_iact_force(r2, dx, hi, hj, pi, pj, a, H);
}
}
}
diff --git a/tests/makeInput.py b/tests/makeInput.py
index c3544d7a25..43fe3da976 100644
--- a/tests/makeInput.py
+++ b/tests/makeInput.py
@@ -27,18 +27,19 @@ from numpy import *
periodic = 1 # 1 For periodic box
boxSize = 1.0
L = 4 # Number of particles along one axis
-rho = 2.0 # Density
+density = 2.0 # Density
P = 1.0 # Pressure
gamma = 5.0 / 3.0 # Gas adiabatic index
+material = 0 # Ideal gas
fileName = "input.hdf5"
# ---------------------------------------------------
numPart = L ** 3
-mass = boxSize ** 3 * rho / numPart
-internalEnergy = P / ((gamma - 1.0) * rho)
+mass = boxSize ** 3 * density / numPart
+internalEnergy = P / ((gamma - 1.0) * density)
# chemistry data
-he_density = rho * 0.24
+he_density = density * 0.24
# Generate particles
coords = zeros((numPart, 3))
@@ -46,7 +47,9 @@ v = zeros((numPart, 3))
m = zeros((numPart, 1))
h = zeros((numPart, 1))
u = zeros((numPart, 1))
+rho = zeros((numPart, 1))
ids = zeros((numPart, 1), dtype="L")
+mat = zeros((numPart, 1), dtype="i")
# chemistry data
he = zeros((numPart, 1))
@@ -67,7 +70,9 @@ for i in range(L):
m[index] = mass
h[index] = 2.251 * boxSize / L
u[index] = internalEnergy
+ rho[index] = density
ids[index] = index
+ mat[index] = material
# chemistry data
he[index] = he_density
@@ -114,8 +119,12 @@ ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
ds[()] = h
ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
ds[()] = u
+ds = grp.create_dataset("Density", (numPart, 1), "f")
+ds[()] = rho
ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
ds[()] = ids
+ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
+ds[()] = mat
# chemistry
ds = grp.create_dataset("HeDensity", (numPart, 1), "f")
ds[()] = he
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 063d414098..128c36480f 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -115,7 +115,9 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH) || \
defined(SPHENIX_SPH) || defined(PHANTOM_SPH) || defined(GASOLINE_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
-#elif defined(PLANETARY_SPH)
+#elif defined(PLANETARY_SPH) || defined(REMIX_SPH)
+ set_idg_def(&eos.idg_def, 0);
+ part->mat_id = 0;
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH)
part->conserved.energy = pressure / (hydro_gamma_minus_one * density);
@@ -219,6 +221,11 @@ void reset_particles(struct cell *c, struct hydro_space *hs,
hydro_first_init_part(p, &c->hydro.xparts[i]);
p->time_bin = 1;
#endif
+#if defined(REMIX_SPH)
+ p->rho = density;
+ hydro_first_init_part(p, &c->hydro.xparts[i]);
+ p->time_bin = 1;
+#endif
hydro_init_part(p, hs);
adaptive_softening_init_part(p);
@@ -289,6 +296,9 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
#else
part->mass = density * volume / count;
#endif
+#if defined(REMIX_SPH)
+ part->rho = density;
+#endif
set_velocity(part, vel, size);
set_energy_state(part, press, size, density);
@@ -379,7 +389,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
#if defined(MINIMAL_SPH) || defined(PLANETARY_SPH) || \
defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH) || \
defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN) || \
- defined(GASOLINE_SPH)
+ defined(GASOLINE_SPH) || defined(REMIX_SPH)
0.f,
#elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
main_cell->hydro.parts[pid].viscosity.div_v,
@@ -442,11 +452,19 @@ void runner_dopair1_branch_density(struct runner *r, struct cell *ci,
void runner_doself1_branch_density(struct runner *r, struct cell *ci,
int limit_h_min, int limit_h_max);
#ifdef EXTRA_HYDRO_LOOP
+#ifdef EXTRA_HYDRO_LOOP_TYPE2
+void runner_dopair2_branch_gradient(struct runner *r, struct cell *ci,
+ struct cell *cj, int limit_h_min,
+ int limit_h_max);
+void runner_doself2_branch_gradient(struct runner *r, struct cell *ci,
+ int limit_h_min, int limit_h_max);
+#else
void runner_dopair1_branch_gradient(struct runner *r, struct cell *ci,
struct cell *cj, int limit_h_min,
int limit_h_max);
void runner_doself1_branch_gradient(struct runner *r, struct cell *ci,
int limit_h_min, int limit_h_max);
+#endif /* EXTRA_HYDRO_LOOP_TYPE2 */
#endif /* EXTRA_HYDRO LOOP */
void runner_dopair2_branch_force(struct runner *r, struct cell *ci,
struct cell *cj, int limit_h_min,
@@ -750,9 +768,15 @@ int main(int argc, char *argv[]) {
struct cell *cj = cells[iii * 25 + jjj * 5 + kkk];
- if (cj > ci)
+ if (cj > ci) {
+#ifdef EXTRA_HYDRO_LOOP_TYPE2
+ runner_dopair2_branch_gradient(
+ &runner, ci, cj, /*limit_h_min=*/0, /*limit_h_max=*/0);
+#else
runner_dopair1_branch_gradient(
&runner, ci, cj, /*limit_h_min=*/0, /*limit_h_max=*/0);
+#endif
+ }
}
}
}
@@ -761,9 +785,15 @@ int main(int argc, char *argv[]) {
}
/* And now the self-interaction for the central cells */
- for (int j = 0; j < 27; ++j)
+ for (int j = 0; j < 27; ++j) {
+#ifdef EXTRA_HYDRO_LOOP_TYPE2
+ runner_doself2_branch_gradient(&runner, inner_cells[j], /*limit_h_min=*/0,
+ /*limit_h_max=*/0);
+#else
runner_doself1_branch_gradient(&runner, inner_cells[j], /*limit_h_min=*/0,
/*limit_h_max=*/0);
+#endif
+ }
/* Extra ghost to finish everything on the central cells */
for (int j = 0; j < 27; ++j)
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 1a91613f6a..e90d4e1781 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -160,6 +160,9 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
#else
part->mass = density * volume / count;
#endif
+#if defined(REMIX_SPH)
+ part->rho_evol = density;
+#endif
#if defined(HOPKINS_PE_SPH)
part->entropy = 1.f;
diff --git a/tests/testActivePair.c b/tests/testActivePair.c
index d2eb6b29eb..1a9f5e849c 100644
--- a/tests/testActivePair.c
+++ b/tests/testActivePair.c
@@ -128,6 +128,15 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
part->entropy_one_over_gamma = 1.f;
#elif defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH)
part->conserved.energy = 1.f;
+#elif defined(PLANETARY_SPH)
+ set_idg_def(&eos.idg_def, 0);
+ part->mat_id = 0;
+ part->u = 1.f;
+#elif defined(REMIX_SPH)
+ set_idg_def(&eos.idg_def, 0);
+ part->mat_id = 0;
+ part->u = 1.f;
+ part->rho_evol = 1.f;
#endif
#if defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH)
@@ -282,6 +291,44 @@ void zero_particle_fields_force(
}
p->geometry.volume = 1.0f;
#endif
+#ifdef PLANETARY_SPH
+ p->rho = 1.f;
+ p->density.rho_dh = 0.f;
+ p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
+ p->density.wcount_dh = 0.f;
+ p->density.rot_v[0] = 0.f;
+ p->density.rot_v[1] = 0.f;
+ p->density.rot_v[2] = 0.f;
+ p->density.div_v = 0.f;
+#endif /* PLANETARY_SPH */
+#if defined(REMIX_SPH)
+ p->rho = 1.f;
+ p->m0 = 1.f;
+ p->density.rho_dh = 0.f;
+ p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
+ p->density.wcount_dh = 0.f;
+ p->gradient.m0_bar = 1.f;
+ p->gradient.grad_m0_bar_gradhterm = 0.f;
+ memset(p->grad_m0, 0.f, 3 * sizeof(float));
+ memset(p->dv_norm_kernel, 0.f, 3 * 3 * sizeof(float));
+ memset(p->du_norm_kernel, 0.f, 3 * sizeof(float));
+ memset(p->drho_norm_kernel, 0.f, 3 * sizeof(float));
+ memset(p->dh_norm_kernel, 0.f, 3 * sizeof(float));
+ memset(p->gradient.grad_m0_bar, 0.f, 3 * sizeof(float));
+ memset(p->gradient.m1_bar, 0.f, 3 * sizeof(float));
+ memset(p->gradient.grad_m1_bar, 0.f, 3 * 3 * sizeof(float));
+ memset(p->gradient.grad_m1_bar_gradhterm, 0.f, 3 * sizeof(float));
+ p->gradient.m2_bar.xx = 1.f;
+ p->gradient.m2_bar.yy = 1.f;
+ p->gradient.m2_bar.zz = 1.f;
+ p->gradient.m2_bar.xy = 0.f;
+ p->gradient.m2_bar.xz = 0.f;
+ p->gradient.m2_bar.yz = 0.f;
+ zero_sym_matrix(&p->gradient.grad_m2_bar[0]);
+ zero_sym_matrix(&p->gradient.grad_m2_bar[1]);
+ zero_sym_matrix(&p->gradient.grad_m2_bar[2]);
+ zero_sym_matrix(&p->gradient.grad_m2_bar_gradhterm);
+#endif /* REMIX_SPH */
/* And prepare for a round of force tasks. */
hydro_prepare_force(p, xp, cosmo, hydro_props, pressure_floor, 0., 0.);
diff --git a/tests/testEOS.c b/tests/testEOS.c
index db6fa4fbef..514a8dc73e 100644
--- a/tests/testEOS.c
+++ b/tests/testEOS.c
@@ -48,300 +48,15 @@
#endif
/**
- * @brief Write a list of densities, energies, and resulting pressures to file
- * from an equation of state.
+ * @brief Test planetary equations of state.
*
* WORK IN PROGRESS
*
- * So far only does the Hubbard & MacFarlane (1980) equations of state.
- *
- * Usage:
- * $ ./testEOS (mat_id) (do_output)
- *
- * Sys args (optional):
- * mat_id | int | Material ID, see equation_of_state.h for the options.
- * Default: 201 (= id_HM80_ice).
- *
- * do_output | int | Set 1 to write the output file of rho, u, P values,
- * set 0 for no output. Default: 0.
- *
- * Output text file contains:
- * header
- * num_rho num_u mat_id # Header values
- * rho_0 rho_1 rho_2 ... rho_num_rho # Array of densities, rho
- * u_0 u_1 u_2 ... u_num_u # Array of energies, u
- * P_0_0 P_0_1 ... P_0_num_u # Array of pressures, P(rho, u)
- * P_1_0 ... ... P_1_num_u
- * ... ... ... ...
- * P_num_rho_0 ... P_num_rho_num_u
- * c_0_0 c_0_1 ... c_0_num_u # Array of sound speeds, c(rho,
- * u)
- * c_1_0 ... ... c_1_num_u
- * ... ... ... ...
- * c_num_rho_0 ... c_num_rho_num_u
- *
- * Note that the values tested extend beyond the range that most EOS are
- * designed for (e.g. outside table limits), to help test the EOS in case of
- * unexpected particle behaviour.
- *
+ * *
*/
#ifdef EOS_PLANETARY
-int main(int argc, char *argv[]) {
- float rho, u, log_rho, log_u, P, c;
- struct unit_system us;
- struct swift_params *params =
- (struct swift_params *)malloc(sizeof(struct swift_params));
- if (params == NULL) error("Error allocating memory for the parameter file.");
- const struct phys_const *phys_const = 0; // Unused placeholder
- const float J_kg_to_erg_g = 1e4; // Convert J/kg to erg/g
- char filename[64];
- // Output table params
- const int num_rho = 100, num_u = 100;
- float log_rho_min = logf(1e-4f), log_rho_max = logf(1e3f), // Densities (cgs)
- log_u_min = logf(1e4f),
- log_u_max = logf(1e13f), // Sp. int. energies (SI)
- log_rho_step = (log_rho_max - log_rho_min) / (num_rho - 1.f),
- log_u_step = (log_u_max - log_u_min) / (num_u - 1.f);
- float A1_rho[num_rho], A1_u[num_u];
- // Sys args
- int mat_id_in, do_output;
- // Default sys args
- const int mat_id_def = eos_planetary_id_HM80_ice;
- const int do_output_def = 0;
-
- // Check the number of system arguments and use defaults if not provided
- switch (argc) {
- case 1:
- // Default both
- mat_id_in = mat_id_def;
- do_output = do_output_def;
- break;
-
- case 2:
- // Read mat_id, default do_output
- mat_id_in = atoi(argv[1]);
- do_output = do_output_def;
- break;
-
- case 3:
- // Read both
- mat_id_in = atoi(argv[1]);
- do_output = atoi(argv[2]);
- break;
-
- default:
- error("Invalid number of system arguments!\n");
- mat_id_in = mat_id_def; // Ignored, just here to keep the compiler happy
- do_output = do_output_def;
- };
-
- enum eos_planetary_material_id mat_id =
- (enum eos_planetary_material_id)mat_id_in;
-
- /* Greeting message */
- printf("This is %s\n", package_description());
-
- // Check material ID
- const enum eos_planetary_type_id type =
- (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
-
- // Select the material base type
- switch (type) {
- // Tillotson
- case eos_planetary_type_Til:
- switch (mat_id) {
- case eos_planetary_id_Til_iron:
- printf(" Tillotson iron \n");
- break;
-
- case eos_planetary_id_Til_granite:
- printf(" Tillotson granite \n");
- break;
-
- case eos_planetary_id_Til_water:
- printf(" Tillotson water \n");
- break;
-
- default:
- error("Unknown material ID! mat_id = %d \n", mat_id);
- };
- break;
-
- // Hubbard & MacFarlane (1980)
- case eos_planetary_type_HM80:
- switch (mat_id) {
- case eos_planetary_id_HM80_HHe:
- printf(" Hubbard & MacFarlane (1980) hydrogen-helium atmosphere \n");
- break;
-
- case eos_planetary_id_HM80_ice:
- printf(" Hubbard & MacFarlane (1980) ice mix \n");
- break;
-
- case eos_planetary_id_HM80_rock:
- printf(" Hubbard & MacFarlane (1980) rock mix \n");
- break;
-
- default:
- error("Unknown material ID! mat_id = %d \n", mat_id);
- };
- break;
-
- // SESAME
- case eos_planetary_type_SESAME:
- switch (mat_id) {
- case eos_planetary_id_SESAME_iron:
- printf(" SESAME basalt 7530 \n");
- break;
-
- case eos_planetary_id_SESAME_basalt:
- printf(" SESAME basalt 7530 \n");
- break;
-
- case eos_planetary_id_SESAME_water:
- printf(" SESAME water 7154 \n");
- break;
-
- case eos_planetary_id_SS08_water:
- printf(" Senft & Stewart (2008) SESAME-like water \n");
- break;
-
- default:
- error("Unknown material ID! mat_id = %d \n", mat_id);
- };
- break;
-
- default:
- error("Unknown material type! mat_id = %d \n", mat_id);
- }
-
- // Convert to internal units
- // Earth masses and radii
- // units_init(&us, 5.9724e27, 6.3710e8, 1.f, 1.f, 1.f);
- // SI
- units_init(&us, 1000.f, 100.f, 1.f, 1.f, 1.f);
- log_rho_min -= logf(units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY));
- log_rho_max -= logf(units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY));
- log_u_min += logf(J_kg_to_erg_g / units_cgs_conversion_factor(
- &us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
- log_u_max += logf(J_kg_to_erg_g / units_cgs_conversion_factor(
- &us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
-
- // Set the input parameters
- // Which EOS to initialise
- parser_set_param(params, "EoS:planetary_use_Til:1");
- parser_set_param(params, "EoS:planetary_use_HM80:1");
- parser_set_param(params, "EoS:planetary_use_SESAME:1");
- // Table file names
- parser_set_param(params,
- "EoS:planetary_HM80_HHe_table_file:"
- "../examples/planetary_HM80_HHe.txt");
- parser_set_param(params,
- "EoS:planetary_HM80_ice_table_file:"
- "../examples/planetary_HM80_ice.txt");
- parser_set_param(params,
- "EoS:planetary_HM80_rock_table_file:"
- "../examples/planetary_HM80_rock.txt");
- parser_set_param(params,
- "EoS:planetary_SESAME_iron_table_file:"
- "../examples/planetary_SESAME_iron_2140.txt");
- parser_set_param(params,
- "EoS:planetary_SESAME_basalt_table_file:"
- "../examples/planetary_SESAME_basalt_7530.txt");
- parser_set_param(params,
- "EoS:planetary_SESAME_water_table_file:"
- "../examples/planetary_SESAME_water_7154.txt");
- parser_set_param(params,
- "EoS:planetary_SS08_water_table_file:"
- "../examples/planetary_SS08_water.txt");
-
- // Initialise the EOS materials
- eos_init(&eos, phys_const, &us, params);
-
- // Manual debug testing
- if (1) {
- printf("\n ### MANUAL DEBUG TESTING ### \n");
-
- rho = 5960;
- u = 1.7e8;
- P = gas_pressure_from_internal_energy(rho, u, eos_planetary_id_HM80_ice);
- printf("u = %.2e, rho = %.2e, P = %.2e \n", u, rho, P);
-
- return 0;
- }
-
- // Output file
- sprintf(filename, "testEOS_rho_u_P_c_%d.txt", mat_id);
- FILE *f = fopen(filename, "w");
- if (f == NULL) {
- printf("Could not open output file!\n");
- exit(EXIT_FAILURE);
- }
-
- if (do_output == 1) {
- fprintf(f, "Density Sp.Int.Energy mat_id \n");
- fprintf(f, "%d %d %d \n", num_rho, num_u, mat_id);
- }
-
- // Densities
- log_rho = log_rho_min;
- for (int i = 0; i < num_rho; i++) {
- A1_rho[i] = exp(log_rho);
- log_rho += log_rho_step;
-
- if (do_output == 1)
- fprintf(f, "%.6e ",
- A1_rho[i] * units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY));
- }
- if (do_output == 1) fprintf(f, "\n");
-
- // Sp. int. energies
- log_u = log_u_min;
- for (int i = 0; i < num_u; i++) {
- A1_u[i] = exp(log_u);
- log_u += log_u_step;
-
- if (do_output == 1)
- fprintf(f, "%.6e ",
- A1_u[i] * units_cgs_conversion_factor(
- &us, UNIT_CONV_ENERGY_PER_UNIT_MASS));
- }
- if (do_output == 1) fprintf(f, "\n");
-
- // Pressures
- for (int i = 0; i < num_rho; i++) {
- rho = A1_rho[i];
-
- for (int j = 0; j < num_u; j++) {
- P = gas_pressure_from_internal_energy(rho, A1_u[j], mat_id);
-
- if (do_output == 1)
- fprintf(f, "%.6e ",
- P * units_cgs_conversion_factor(&us, UNIT_CONV_PRESSURE));
- }
-
- if (do_output == 1) fprintf(f, "\n");
- }
-
- // Sound speeds
- for (int i = 0; i < num_rho; i++) {
- rho = A1_rho[i];
-
- for (int j = 0; j < num_u; j++) {
- c = gas_soundspeed_from_internal_energy(rho, A1_u[j], mat_id);
-
- if (do_output == 1)
- fprintf(f, "%.6e ",
- c * units_cgs_conversion_factor(&us, UNIT_CONV_SPEED));
- }
-
- if (do_output == 1) fprintf(f, "\n");
- }
- fclose(f);
-
- return 0;
-}
+int main(int argc, char *argv[]) { return 0; }
#else
int main(int argc, char *argv[]) { return 0; }
#endif
diff --git a/tests/testEOS.sh b/tests/testEOS.sh
deleted file mode 100755
index bcd87eabbf..0000000000
--- a/tests/testEOS.sh
+++ /dev/null
@@ -1,27 +0,0 @@
-#!/bin/bash
-
-echo ""
-
-rm -f testEOS_rho_u_P_*.txt
-
-echo "Running testEOS for each planetary material"
-
-A1_mat_id=(
- 100
- 101
- 102
- 200
- 201
- 202
- 300
- 301
- 302
- 303
-)
-
-for mat_id in "${A1_mat_id[@]}"
-do
- ./testEOS "$mat_id" 1
-done
-
-exit $?
diff --git a/tests/testInteractions.c b/tests/testInteractions.c
index 005c4666b7..960c2f9254 100644
--- a/tests/testInteractions.c
+++ b/tests/testInteractions.c
@@ -114,7 +114,8 @@ void prepare_force(struct part *parts, size_t count) {
#if !defined(GIZMO_MFV_SPH) && !defined(MINIMAL_SPH) && \
!defined(PLANETARY_SPH) && !defined(HOPKINS_PU_SPH) && \
!defined(HOPKINS_PU_SPH_MONAGHAN) && !defined(ANARCHY_PU_SPH) && \
- !defined(SPHENIX_SPH) && !defined(PHANTOM_SPH) && !defined(GASOLINE_SPH)
+ !defined(SPHENIX_SPH) && !defined(PHANTOM_SPH) && \
+ !defined(GASOLINE_SPH) && !defined(REMIX_SPH)
struct part *p;
for (size_t i = 0; i < count; ++i) {
p = &parts[i];
diff --git a/tests/testPeriodicBC.c b/tests/testPeriodicBC.c
index 05918b5b5e..be31b66051 100644
--- a/tests/testPeriodicBC.c
+++ b/tests/testPeriodicBC.c
@@ -139,6 +139,9 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
#else
part->mass = density * volume / count;
#endif
+#if defined(REMIX_SPH)
+ part->rho_evol = density;
+#endif
#if defined(HOPKINS_PE_SPH)
part->entropy = 1.f;
diff --git a/tests/testRiemannTRRS.c b/tests/testRiemannTRRS.c
index 16956bec0f..db988c3021 100644
--- a/tests/testRiemannTRRS.c
+++ b/tests/testRiemannTRRS.c
@@ -18,6 +18,8 @@
******************************************************************************/
#include <config.h>
+#if !defined(RIEMANN_SOLVER_NONE)
+
/* Local headers. */
#include <string.h>
@@ -326,3 +328,7 @@ int main(int argc, char* argv[]) {
return 0;
}
+
+#else
+int main(int argc, char* argv[]) { return 0; }
+#endif /* !RIEMANN_SOLVER_NONE */
--
GitLab