Added (so-far unused) Gizmo files to the master to simplify the unification merge request

src/Makefile.am

@@ -54,7 +54,9 @@ nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel.h
 		 hydro/Default/hydro.h hydro/Default/hydro_iact.h hydro/Default/hydro_io.h \
                  hydro/Default/hydro_debug.h hydro/Default/hydro_part.h \
 		 hydro/Gadget2/hydro.h hydro/Gadget2/hydro_iact.h hydro/Gadget2/hydro_io.h \
-                 hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h
+                 hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \
+		 hydro/Gizmo/hydro.h hydro/Gizmo/hydro_iact.h hydro/Gizmo/hydro_io.h \
+                 hydro/Gizmo/hydro_debug.h hydro/Gizmo/hydro_part.h
 
 # Sources and flags for regular library
 libswiftsim_la_SOURCES = $(AM_SOURCES)

src/hydro/Gizmo/hydro_debug.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 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/>.
+ *
+ ******************************************************************************/
+
+__attribute__((always_inline))
+    INLINE static void hydro_debug_particle(struct part* p, struct xpart* xp) {
+  printf(
+      "x=[%.16e,%.16e,%.16e], "
+      "v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], volume=%.3e\n",
+      p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
+      p->a_hydro[1], p->a_hydro[2], p->geometry.volume);
+}

src/hydro/Gizmo/hydro_io.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 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/>.
+ *
+ ******************************************************************************/
+
+/**
+ * @brief Reads the different particles to the HDF5 file
+ *
+ * @param h_grp The HDF5 group in which to read the arrays.
+ * @param N The number of particles on that MPI rank.
+ * @param N_total The total number of particles (only used in MPI mode)
+ * @param offset The offset of the particles for this MPI rank (only used in MPI
+ *mode)
+ * @param parts The particle array
+ *
+ */
+__attribute__((always_inline)) INLINE static void hydro_read_particles(
+    hid_t h_grp, int N, long long N_total, long long offset,
+    struct part* parts) {
+
+  /* Read arrays */
+  readArray(h_grp, "Coordinates", DOUBLE, N, 3, parts, N_total, offset, x,
+            COMPULSORY);
+  readArray(h_grp, "Velocities", FLOAT, N, 3, parts, N_total, offset, v,
+            COMPULSORY);
+  readArray(h_grp, "Masses", FLOAT, N, 1, parts, N_total, offset,
+            conserved.mass, COMPULSORY);
+  readArray(h_grp, "SmoothingLength", FLOAT, N, 1, parts, N_total, offset, h,
+            COMPULSORY);
+  readArray(h_grp, "InternalEnergy", FLOAT, N, 1, parts, N_total, offset,
+            primitives.P, COMPULSORY);
+  readArray(h_grp, "ParticleIDs", ULONGLONG, N, 1, parts, N_total, offset, id,
+            COMPULSORY);
+  readArray(h_grp, "Acceleration", FLOAT, N, 3, parts, N_total, offset, a_hydro,
+            OPTIONAL);
+  readArray(h_grp, "Density", FLOAT, N, 1, parts, N_total, offset,
+            primitives.rho, OPTIONAL);
+}
+
+/**
+ * @brief Writes the different particles to the HDF5 file
+ *
+ * @param h_grp The HDF5 group in which to write the arrays.
+ * @param fileName The name of the file (unsued in MPI mode).
+ * @param xmfFile The XMF file to write to (unused in MPI mode).
+ * @param N The number of particles on that MPI rank.
+ * @param N_total The total number of particles (only used in MPI mode)
+ * @param mpi_rank The MPI rank of this node (only used in MPI mode)
+ * @param offset The offset of the particles for this MPI rank (only used in MPI
+ *mode)
+ * @param parts The particle array
+ * @param us The unit system to use
+ *
+ */
+__attribute__((always_inline)) INLINE static void hydro_write_particles(
+    hid_t h_grp, char* fileName, FILE* xmfFile, int N, long long N_total,
+    int mpi_rank, long long offset, struct part* parts, struct UnitSystem* us) {
+
+  /* Write arrays */
+  writeArray(h_grp, fileName, xmfFile, "Coordinates", DOUBLE, N, 3, parts,
+             N_total, mpi_rank, offset, x, us, UNIT_CONV_LENGTH);
+  writeArray(h_grp, fileName, xmfFile, "Velocities", FLOAT, N, 3, parts,
+             N_total, mpi_rank, offset, v, us, UNIT_CONV_SPEED);
+  writeArray(h_grp, fileName, xmfFile, "Masses", FLOAT, N, 1, parts, N_total,
+             mpi_rank, offset, conserved.mass, us, UNIT_CONV_MASS);
+  writeArray(h_grp, fileName, xmfFile, "SmoothingLength", FLOAT, N, 1, parts,
+             N_total, mpi_rank, offset, h, us, UNIT_CONV_LENGTH);
+  writeArray(h_grp, fileName, xmfFile, "InternalEnergy", FLOAT, N, 1, parts,
+             N_total, mpi_rank, offset, primitives.P, us,
+             UNIT_CONV_ENTROPY_PER_UNIT_MASS);
+  writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N, 1, parts,
+             N_total, mpi_rank, offset, id, us, UNIT_CONV_NO_UNITS);
+  writeArray(h_grp, fileName, xmfFile, "Acceleration", FLOAT, N, 3, parts,
+             N_total, mpi_rank, offset, a_hydro, us, UNIT_CONV_ACCELERATION);
+  writeArray(h_grp, fileName, xmfFile, "Density", FLOAT, N, 1, parts, N_total,
+             mpi_rank, offset, primitives.rho, us, UNIT_CONV_DENSITY);
+}
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grpsph The HDF5 group in which to write
+ */
+void writeSPHflavour(hid_t h_grpsph) {
+
+  /* Kernel function description */
+  writeAttribute_s(h_grpsph, "Kernel", kernel_name);
+  writeAttribute_f(h_grpsph, "Kernel eta", const_eta_kernel);
+  writeAttribute_f(h_grpsph, "Weighted N_ngb", kernel_nwneigh);
+  writeAttribute_f(h_grpsph, "Delta N_ngb", const_delta_nwneigh);
+  writeAttribute_f(h_grpsph, "Hydro gamma", const_hydro_gamma);
+
+  /* Viscosity and thermal conduction */
+  writeAttribute_s(h_grpsph, "Thermal Conductivity Model",
+                   "(No treatment) Legacy Gadget-2 as in Springel (2005)");
+  writeAttribute_s(h_grpsph, "Viscosity Model",
+                   "Legacy Gadget-2 as in Springel (2005)");
+  writeAttribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
+  writeAttribute_f(h_grpsph, "Viscosity beta", 3.f);
+
+  /* Time integration properties */
+  writeAttribute_f(h_grpsph, "CFL parameter", const_cfl);
+  writeAttribute_f(h_grpsph, "Maximal ln(Delta h) change over dt",
+                   const_ln_max_h_change);
+  writeAttribute_f(h_grpsph, "Maximal Delta h change over dt",
+                   exp(const_ln_max_h_change));
+}

src/hydro/Gizmo/hydro_part.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *                    Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/* Some standard headers. */
+#include <stdlib.h>
+
+#define GFLOAT float
+
+/* Extra particle data not needed during the computation. */
+struct xpart {
+
+  /* Old position, at last tree rebuild. */
+  double x_old[3];
+
+  /* Velocity at the last full step. */
+  float v_full[3];
+
+  /* Entropy at the half-step. */
+  float u_hdt;
+
+  /* Old density. */
+  float omega;
+
+} __attribute__((aligned(xpart_align)));
+
+/* Data of a single particle. */
+struct part {
+
+  /* Particle position. */
+  double x[3];
+
+  /* Particle velocity. */
+  float v[3];
+
+  /* Particle acceleration. */
+  float a_hydro[3];
+
+  float mass;  // MATTHIEU
+  float h_dt;
+  float rho;
+  float rho_dh;
+
+  /* Particle cutoff radius. */
+  float h;
+
+  /* Particle time of beginning of time-step. */
+  int ti_begin;
+
+  /* Particle time of end of time-step. */
+  int ti_end;
+
+  /* The primitive hydrodynamical variables */
+  struct {
+
+    /* fluid velocity */
+    GFLOAT v[3];
+
+    /* density */
+    GFLOAT rho;
+
+    /* pressure */
+    GFLOAT P;
+
+    struct {
+
+      GFLOAT rho[3];
+
+      GFLOAT v[3][3];
+
+      GFLOAT P[3];
+
+    } gradients;
+
+    struct {
+
+      /* extreme values among the neighbours */
+      GFLOAT rho[2];
+
+      GFLOAT v[3][2];
+
+      GFLOAT P[2];
+
+      /* maximal distance to all neighbouring faces */
+      float maxr;
+
+    } limiter;
+
+  } primitives;
+
+  /* The conserved hydrodynamical variables */
+  struct {
+
+    /* fluid momentum */
+    GFLOAT momentum[3];
+
+    /* fluid mass */
+    GFLOAT mass;
+
+    /* fluid energy */
+    GFLOAT energy;
+
+  } conserved;
+
+  /* Geometrical quantities used for hydro */
+  struct {
+
+    /* volume of the particle */
+    float volume;
+
+    /* gradient matrix */
+    float matrix_E[3][3];
+
+  } geometry;
+
+  struct {
+
+    float vmax;
+
+  } timestepvars;
+
+  /* Quantities used during the density loop */
+  struct {
+
+    /* Particle velocity divergence. */
+    float div_v;
+
+    /* Derivative of particle number density. */
+    float wcount_dh;
+
+    /* Particle velocity curl. */
+    float curl_v[3];
+
+    /* Particle number density. */
+    float wcount;
+
+  } density;
+
+  /* Particle ID. */
+  unsigned long long id;
+
+  /* Associated gravitas. */
+  struct gpart *gpart;
+
+} __attribute__((aligned(part_align)));

src/riemann.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ *
+ * 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_RIEMANN_H
+#define SWIFT_RIEMANN_H
+
+/* gives us const_hydro_gamma and tells us which floating point type to use */
+#include "const.h"
+#include "math.h"
+#include "stdio.h"
+#include "float.h"
+#include "stdlib.h"
+#include "error.h"
+
+#define HLLC_SOLVER
+
+#ifdef EXACT_SOLVER
+#include "riemann/riemann_exact.h"
+#endif
+
+#ifdef TRRS_SOLVER
+#include "riemann/riemann_trrs.h"
+#endif
+
+#ifdef HLLC_SOLVER
+#include "riemann/riemann_hllc.h"
+#endif
+
+#endif /* SWIFT_RIEMANN_H */

src/riemann/riemann_hllc.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ *
+ * 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_RIEMANN_HLLC_H
+#define SWIFT_RIEMANN_HLLC_H
+
+__attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
+    GFLOAT *WL, GFLOAT *WR, float *n, float *vij, GFLOAT *totflux) {
+
+  GFLOAT uL, uR, aL, aR;
+  GFLOAT rhobar, abar, pPVRS, pstar, qL, qR, SL, SR, Sstar;
+  GFLOAT v2, eL, eR;
+  GFLOAT UstarL[5], UstarR[5];
+
+  /* Handle pure vacuum */
+  if (!WL[0] && !WR[0]) {
+    totflux[0] = 0.;
+    totflux[1] = 0.;
+    totflux[2] = 0.;
+    totflux[3] = 0.;
+    totflux[4] = 0.;
+    return;
+  }
+
+  /* STEP 0: obtain velocity in interface frame */
+  uL = WL[1] * n[0] + WL[2] * n[1] + WL[3] * n[2];
+  uR = WR[1] * n[0] + WR[2] * n[1] + WR[3] * n[2];
+  aL = sqrtf(const_hydro_gamma * WL[4] / WL[0]);
+  aR = sqrtf(const_hydro_gamma * WR[4] / WR[0]);
+
+  /* Handle vacuum: vacuum does not require iteration and is always exact */
+  if (!WL[0] || !WR[0]) {
+    error("Vacuum not yet supported");
+  }
+  if (2. * aL / (const_hydro_gamma - 1.) + 2. * aR / (const_hydro_gamma - 1.) <
+      fabs(uL - uR)) {
+    error("Vacuum not yet supported");
+  }
+
+  /* STEP 1: pressure estimate */
+  rhobar = 0.5 * (WL[0] + WR[0]);
+  abar = 0.5 * (aL + aR);
+  pPVRS = 0.5 * (WL[4] + WR[4]) - 0.5 * (uR - uL) * rhobar * abar;
+  pstar = fmaxf(0., pPVRS);
+
+  /* STEP 2: wave speed estimates
+     all these speeds are along the interface normal, since uL and uR are */
+  qL = 1.;
+  if (pstar > WL[4]) {
+    qL = sqrtf(1. + 0.5 * (const_hydro_gamma + 1.) / const_hydro_gamma *
+                        (pstar / WL[4] - 1.));
+  }
+  qR = 1.;
+  if (pstar > WR[4]) {
+    qR = sqrtf(1. + 0.5 * (const_hydro_gamma + 1.) / const_hydro_gamma *
+                        (pstar / WR[4] - 1.));
+  }
+  SL = uL - aL * qL;
+  SR = uR + aR * qR;
+  Sstar = (WR[4] - WL[4] + WL[0] * uL * (SL - uL) - WR[0] * uR * (SR - uR)) /
+          (WL[0] * (SL - uL) - WR[0] * (SR - uR));
+
+  /* STEP 3: HLLC flux in a frame moving with the interface velocity */
+  if (Sstar >= 0.) {
+    /* flux FL */
+    totflux[0] = WL[0] * uL;
+    /* these are the actual correct fluxes in the boosted lab frame
+       (not rotated to interface frame) */
+    totflux[1] = WL[0] * WL[1] * uL + WL[4] * n[0];
+    totflux[2] = WL[0] * WL[2] * uL + WL[4] * n[1];
+    totflux[3] = WL[0] * WL[2] * uL + WL[4] * n[2];
+    v2 = WL[1] * WL[1] + WL[2] * WL[2] + WL[3] * WL[3];
+    eL = WL[4] / (const_hydro_gamma - 1.) / WL[0] + 0.5 * v2;
+    totflux[4] = WL[0] * eL * uL + WL[4] * uL;
+    if (SL < 0.) {
+      /* add flux FstarL */
+      UstarL[0] = 1.;
+      /* we need UstarL in the lab frame:
+         subtract the velocity in the interface frame from the lab frame
+         velocity and then add Sstar in interface frame */
+      UstarL[1] = WL[1] + (Sstar - uL) * n[0];
+      UstarL[2] = WL[2] + (Sstar - uL) * n[1];
+      UstarL[3] = WL[3] + (Sstar - uL) * n[2];
+      UstarL[4] = eL + (Sstar - uL) * (Sstar + WL[4] / (WL[0] * (SL - uL)));
+      UstarL[0] *= WL[0] * (SL - uL) / (SL - Sstar);
+      UstarL[1] *= WL[0] * (SL - uL) / (SL - Sstar);
+      UstarL[2] *= WL[0] * (SL - uL) / (SL - Sstar);
+      UstarL[3] *= WL[0] * (SL - uL) / (SL - Sstar);
+      UstarL[4] *= WL[0] * (SL - uL) / (SL - Sstar);
+      totflux[0] += SL * (UstarL[0] - WL[0]);
+      totflux[1] += SL * (UstarL[1] - WL[0] * WL[1]);
+      totflux[2] += SL * (UstarL[2] - WL[0] * WL[2]);
+      totflux[3] += SL * (UstarL[3] - WL[0] * WL[3]);
+      totflux[4] += SL * (UstarL[4] - WL[0] * eL);
+    }
+  } else {
+    /* flux FR */
+    totflux[0] = WR[0] * uR;
+    totflux[1] = WR[0] * WR[1] * uR + WR[4] * n[0];
+    totflux[2] = WR[0] * WR[2] * uR + WR[4] * n[1];
+    totflux[3] = WR[0] * WR[3] * uR + WR[4] * n[2];
+    v2 = WR[1] * WR[1] + WR[2] * WR[2] + WR[3] * WR[3];
+    eR = WR[4] / (const_hydro_gamma - 1.) / WR[0] + 0.5 * v2;
+    totflux[4] = WR[0] * eR * uR + WR[4] * uR;
+    if (SR > 0.) {
+      /* add flux FstarR */
+      UstarR[0] = 1.;
+      UstarR[1] = WR[1] + (Sstar - uR) * n[0];
+      UstarR[2] = WR[2] + (Sstar - uR) * n[1];
+      UstarR[3] = WR[3] + (Sstar - uR) * n[2];
+      UstarR[4] = eR + (Sstar - uR) * (Sstar + WR[4] / (WR[0] * (SR - uR)));
+      UstarR[0] *= WR[0] * (SR - uR) / (SR - Sstar);
+      UstarR[1] *= WR[0] * (SR - uR) / (SR - Sstar);
+      UstarR[2] *= WR[0] * (SR - uR) / (SR - Sstar);
+      UstarR[3] *= WR[0] * (SR - uR) / (SR - Sstar);
+      UstarR[4] *= WR[0] * (SR - uR) / (SR - Sstar);
+      totflux[0] += SR * (UstarR[0] - WR[0]);
+      totflux[1] += SR * (UstarR[1] - WR[0] * WR[1]);
+      totflux[2] += SR * (UstarR[2] - WR[0] * WR[2]);
+      totflux[3] += SR * (UstarR[3] - WR[0] * WR[3]);
+      totflux[4] += SR * (UstarR[4] - WR[0] * eR);
+    }
+  }
+
+  /* deboost to lab frame
+     we add the flux contribution due to the movement of the interface
+     the density flux is unchanged
+     we add the extra velocity flux due to the absolute motion of the fluid
+     similarly, we need to add the energy fluxes due to the absolute motion */
+  v2 = vij[0] * vij[0] + vij[1] * vij[1] + vij[2] * vij[2];
+  totflux[1] += vij[0] * totflux[0];
+  totflux[2] += vij[1] * totflux[0];
+  totflux[3] += vij[2] * totflux[0];
+  totflux[4] += vij[0] * totflux[1] + vij[1] * totflux[2] +
+                vij[2] * totflux[3] + 0.5 * v2 * totflux[0];
+}
+
+#endif /* SWIFT_RIEMANN_HLLC_H */

src/riemann/riemann_trrs.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ *
+ * 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_RIEMANN_TRRS_H
+#define SWIFT_RIEMANN_TRRS_H
+
+/* frequently used combinations of const_hydro_gamma */
+#define const_riemann_gp1d2g \
+  (0.5f * (const_hydro_gamma + 1.0f) / const_hydro_gamma)
+#define const_riemann_gm1d2g \
+  (0.5f * (const_hydro_gamma - 1.0f) / const_hydro_gamma)
+#define const_riemann_gm1dgp1 \
+  ((const_hydro_gamma - 1.0f) / (const_hydro_gamma + 1.0f))
+#define const_riemann_tdgp1 (2.0f / (const_hydro_gamma + 1.0f))
+#define const_riemann_tdgm1 (2.0f / (const_hydro_gamma - 1.0f))
+#define const_riemann_gm1d2 (0.5f * (const_hydro_gamma - 1.0f))
+#define const_riemann_tgdgm1 \
+  (2.0f * const_hydro_gamma / (const_hydro_gamma - 1.0f))
+#define const_riemann_ginv (1.0f / const_hydro_gamma)
+
+/**
+ * @brief Solve the Riemann problem using the Two Rarefaction Riemann Solver
+ *
+ * By assuming 2 rarefaction waves, we can analytically solve for the pressure
+ * and velocity in the intermediate region, eliminating the iterative procedure.
+ *
+ * According to Toro: "The two-rarefaction approximation is generally quite
+ * robust; (...) The TRRS is in fact exact when both non-linear waves are
+ * actually rarefaction waves."
+ *
+ * @param WL The left state vector
+ * @param WR The right state vector
+ * @param Whalf Empty state vector in which the result will be stored
+ * @param n_unit Normal vector of the interface
+ */
+__attribute__((always_inline)) INLINE static void riemann_solver_solve(
+    GFLOAT* WL, GFLOAT* WR, GFLOAT* Whalf, float* n_unit) {
+  GFLOAT aL, aR;
+  GFLOAT PLR;
+  GFLOAT vL, vR;
+  GFLOAT ustar, pstar;
+  GFLOAT vhalf;
+  GFLOAT pdpR, SHR, STR;
+  GFLOAT pdpL, SHL, STL;
+
+  /* calculate the velocities along the interface normal */
+  vL = WL[1] * n_unit[0] + WL[2] * n_unit[1] + WL[3] * n_unit[2];
+  vR = WR[1] * n_unit[0] + WR[2] * n_unit[1] + WR[3] * n_unit[2];
+
+  /* calculate the sound speeds */
+  aL = sqrtf(const_hydro_gamma * WL[4] / WL[0]);
+  aR = sqrtf(const_hydro_gamma * WR[4] / WR[0]);
+
+  /* calculate the velocity and pressure in the intermediate state */
+  PLR = pow(WL[4] / WR[4], const_riemann_gm1d2g);
+  ustar = (PLR * vL / aL + vR / aR + const_riemann_tdgm1 * (PLR - 1.0f)) /
+          (PLR / aL + 1.0f / aR);
+  pstar = 0.5f * (WL[4] * pow(1.0f + const_riemann_gm1d2 / aL * (vL - ustar),
+                              const_riemann_tgdgm1) +
+                  WR[4] * pow(1.0f + const_riemann_gm1d2 / aR * (ustar - vR),
+                              const_riemann_tgdgm1));
+
+  /* sample the solution */
+  if (ustar < 0.0f) {
+    /* right state */
+    Whalf[1] = WR[1];
+    Whalf[2] = WR[2];
+    Whalf[3] = WR[3];
+    pdpR = pstar / WR[4];
+    /* always a rarefaction wave, that's the approximation */
+    SHR = vR + aR;
+    if (SHR > 0.0f) {
+      STR = ustar + aR * pow(pdpR, const_riemann_gm1d2g);
+      if (STR <= 0.0f) {
+        Whalf[0] =
+            WR[0] * pow(const_riemann_tdgp1 - const_riemann_gm1dgp1 / aR * vR,
+                        const_riemann_tdgm1);
+        vhalf = const_riemann_tdgp1 * (-aR + const_riemann_gm1d2 * vR) - vR;
+        Whalf[4] =
+            WR[4] * pow(const_riemann_tdgp1 - const_riemann_gm1dgp1 / aR * vR,
+                        const_riemann_tgdgm1);
+      } else {
+        Whalf[0] = WR[0] * pow(pdpR, const_riemann_ginv);
+        vhalf = ustar - vR;
+        Whalf[4] = pstar;
+      }
+    } else {
+      Whalf[0] = WR[0];
+      vhalf = 0.0f;
+      Whalf[4] = WR[4];
+    }
+  } else {
+    /* left state */
+    Whalf[1] = WL[1];
+    Whalf[2] = WL[2];
+    Whalf[3] = WL[3];
+    pdpL = pstar / WL[4];
+    /* rarefaction wave */
+    SHL = vL - aL;
+    if (SHL < 0.0f) {
+      STL = ustar - aL * pow(pdpL, const_riemann_gm1d2g);
+      if (STL > 0.0f) {
+        Whalf[0] =
+            WL[0] * pow(const_riemann_tdgp1 + const_riemann_gm1dgp1 / aL * vL,
+                        const_riemann_tdgm1);
+        vhalf = const_riemann_tdgp1 * (aL + const_riemann_gm1d2 * vL) - vL;
+        Whalf[4] =
+            WL[4] * pow(const_riemann_tdgp1 + const_riemann_gm1dgp1 / aL * vL,
+                        const_riemann_tgdgm1);
+      } else {
+        Whalf[0] = WL[0] * pow(pdpL, const_riemann_ginv);
+        vhalf = ustar - vL;
+        Whalf[4] = pstar;
+      }
+    } else {
+      Whalf[0] = WL[0];
+      vhalf = 0.0f;
+      Whalf[4] = WL[4];
+    }
+  }
+
+  /* add the velocity solution along the interface normal to the velocities */
+  Whalf[1] += vhalf * n_unit[0];
+  Whalf[2] += vhalf * n_unit[1];
+  Whalf[3] += vhalf * n_unit[2];
+}
+
+#endif /* SWIFT_RIEMANN_TRRS_H */