diff --git a/src/const.h b/src/const.h
index 736ba0155e03dc0f595ac06edfd0763b38302e9c..b6182439aa49f36db503a7868ba9e301a29d2593 100644
--- a/src/const.h
+++ b/src/const.h
@@ -51,8 +51,9 @@
/* SPH variant to use */
//#define MINIMAL_SPH
-#define GADGET2_SPH
+//#define GADGET2_SPH
//#define DEFAULT_SPH
+#define GIZMO_SPH
/* Self gravity stuff. */
#define const_gravity_multipole_order 2
diff --git a/src/debug.c b/src/debug.c
index 487fd53e74399ef7bd1802704adbf84bbc3dc0a3..ca77745ec9f905d85403445ff3d8ef2de16034e7 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -43,6 +43,8 @@
#include "./hydro/Gadget2/hydro_debug.h"
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_debug.h"
+#elif defined(GIZMO_SPH)
+#include "./hydro/Gizmo/hydro_debug.h"
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/hydro.h b/src/hydro.h
index b2ae9d57c399ecea818e9f3dc7db238e01487a9a..793b82ea1859fe939484a7b6016a33c404fd5aca 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -38,6 +38,10 @@
#include "./hydro/Default/hydro.h"
#include "./hydro/Default/hydro_iact.h"
#define SPH_IMPLEMENTATION "Default version of SPH"
+#elif defined(GIZMO_SPH)
+#include "./hydro/Gizmo/hydro.h"
+#include "./hydro/Gizmo/hydro_iact.h"
+#define SPH_IMPLEMENTATION "Mesh-free hydrodynamics as in GIZMO (Hopkins 2015)"
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index f69dc3f1798f014e895c4a63760805b1739cec94..95231e2d8fe0d7795f5f7eaa3fdaef18bf3ad292 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -17,6 +17,9 @@
*
******************************************************************************/
+#include <float.h>
+#include "adiabatic_index.h"
+
/**
* @brief Computes the hydro time-step of a given particle
*
@@ -25,9 +28,12 @@
*
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
- struct part* p, struct xpart* xp) {
+ const struct part* restrict p, const struct xpart* restrict xp,
+ const struct hydro_props* restrict hydro_properties) {
- return const_cfl * p->h / fabs(p->timestepvars.vmax);
+ const float CFL_condition = hydro_properties->CFL_condition;
+
+ return CFL_condition * p->h / fabsf(p->timestepvars.vmax);
}
/**
@@ -59,7 +65,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->primitives.v[1] = p->conserved.momentum[1] / p->conserved.mass;
p->primitives.v[2] = p->conserved.momentum[2] / p->conserved.mass;
p->primitives.P =
- (const_hydro_gamma - 1.) *
+ hydro_gamma_minus_one *
(p->conserved.energy -
0.5 * (p->conserved.momentum[0] * p->conserved.momentum[0] +
p->conserved.momentum[1] * p->conserved.momentum[1] +
@@ -158,15 +164,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
float detE, volume;
float E[3][3];
- GFLOAT m, momentum[3], energy;
+ float m, momentum[3], energy;
#ifndef THERMAL_ENERGY
- GFLOAT momentum2;
+ float momentum2;
#endif
#if defined(SPH_GRADIENTS) && defined(SLOPE_LIMITER)
- GFLOAT gradrho[3], gradv[3][3], gradP[3];
- GFLOAT gradtrue, gradmax, gradmin, alpha;
+ float gradrho[3], gradv[3][3], gradP[3];
+ float gradtrue, gradmax, gradmin, alpha;
#endif
/* Final operation on the geometry. */
@@ -374,9 +380,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
p->primitives.v[2] = momentum[2] / m;
#ifndef THERMAL_ENERGY
p->primitives.P =
- (const_hydro_gamma - 1.) * (energy - 0.5 * momentum2 / m) / volume;
+ hydro_gamma_minus_one * (energy - 0.5 * momentum2 / m) / volume;
#else
- p->primitives.P = (const_hydro_gamma - 1.) * energy / volume;
+ p->primitives.P = hydro_gamma_minus_one * energy / volume;
#endif
}
}
@@ -392,8 +398,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
#ifndef SPH_GRADIENTS
float h, ih, ih2, ih3;
#ifdef SLOPE_LIMITER
- GFLOAT gradrho[3], gradv[3][3], gradP[3];
- GFLOAT gradtrue, gradmax, gradmin, alpha;
+ float gradrho[3], gradv[3][3], gradP[3];
+ float gradtrue, gradmax, gradmin, alpha;
#endif
/* add kernel normalization to gradients */
@@ -575,10 +581,10 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part* p) {
float volume;
- GFLOAT m;
- GFLOAT momentum[3];
+ float m;
+ float momentum[3];
#ifndef THERMAL_ENERGY
- GFLOAT momentum2;
+ float momentum2;
#endif
volume = p->geometry.volume;
@@ -587,7 +593,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
p->primitives.v[1] = p->v[1];
p->primitives.v[2] = p->v[2];
/* P actually contains internal energy at this point */
- p->primitives.P *= (const_hydro_gamma - 1.) * p->primitives.rho;
+ p->primitives.P *= hydro_gamma_minus_one * p->primitives.rho;
p->conserved.mass = m = p->primitives.rho * volume;
p->conserved.momentum[0] = momentum[0] = m * p->primitives.v[0];
@@ -597,9 +603,9 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
momentum2 = momentum[0] * momentum[0] + momentum[1] * momentum[1] +
momentum[2] * momentum[2];
p->conserved.energy =
- p->primitives.P / (const_hydro_gamma - 1.) * volume + 0.5 * momentum2 / m;
+ p->primitives.P / hydro_gamma_minus_one * volume + 0.5 * momentum2 / m;
#else
- p->conserved.energy = p->primitives.P / (const_hydro_gamma - 1.) * volume;
+ p->conserved.energy = p->primitives.P / hydro_gamma_minus_one * volume;
#endif
}
@@ -614,7 +620,51 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
struct part* p) {}
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt, float half_dt) {}
+
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
- struct part* p) {
- return 0.f;
+ const struct part* restrict p, float dt) {
+
+ return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+}
+
+/**
+ * @brief Returns the entropy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+ const struct part* restrict p, float dt) {
+
+ return p->primitives.P / pow_gamma(p->primitives.rho);
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
+ const struct part* restrict p, float dt) {
+
+ return sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+}
+
+/**
+ * @brief Returns the pressure of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+ const struct part* restrict p, float dt) {
+
+ return p->primitives.P;
}
diff --git a/src/hydro/Gizmo/hydro_debug.h b/src/hydro/Gizmo/hydro_debug.h
index 365d85a2f651cf98b0713e8d82f11ae70fa9beaa..e42ebb875f257001bf346c82a8c7559b47ef6b86 100644
--- a/src/hydro/Gizmo/hydro_debug.h
+++ b/src/hydro/Gizmo/hydro_debug.h
@@ -18,7 +18,7 @@
******************************************************************************/
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
- struct part* p, struct xpart* xp) {
+ const struct part* p, const struct xpart* xp) {
printf(
"x=[%.16e,%.16e,%.16e], "
"v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], volume=%.3e\n",
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index 30a8d6cbebc851b44a5ee2339950aec9e15057c0..7f5289c8d519ebc32ba3816fe253776e3350d1ed 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -19,6 +19,7 @@
*
******************************************************************************/
+#include "adiabatic_index.h"
#include "riemann.h"
#define USE_GRADIENTS
@@ -26,7 +27,7 @@
/* #define PRINT_ID 0 */
/* this corresponds to task_subtype_hydro_loop1 */
-__attribute__((always_inline)) INLINE static void runner_iact_hydro_loop1(
+__attribute__((always_inline)) INLINE static void runner_iact_density(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
float r = sqrtf(r2);
@@ -194,9 +195,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_hydro_loop1(
}
/* this corresponds to task_subtype_hydro_loop1 */
-__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_hydro_loop1(float r2, float *dx, float hi, float hj,
- struct part *pi, struct part *pj) {
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
float r;
float xi;
@@ -296,7 +296,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_hydro_loop2(
int k, l;
float Bi[3][3];
float Bj[3][3];
- GFLOAT Wi[5], Wj[5];
+ float Wi[5], Wj[5];
/* Initialize local variables */
for (k = 0; k < 3; k++) {
@@ -497,7 +497,7 @@ runner_iact_nonsym_hydro_loop2(float r2, float *dx, float hi, float hj,
float wi, wi_dx;
int k, l;
float Bi[3][3];
- GFLOAT Wi[5], Wj[5];
+ float Wi[5], Wj[5];
/* Initialize local variables */
for (k = 0; k < 3; k++) {
@@ -619,13 +619,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float xij_i[3], xfac, xijdotdx;
float vmax, dvdotdx;
float vi[3], vj[3], vij[3];
- GFLOAT Wi[5], Wj[5]; //, Whalf[5];
+ float Wi[5], Wj[5]; //, Whalf[5];
#ifdef USE_GRADIENTS
- GFLOAT dWi[5], dWj[5];
+ float dWi[5], dWj[5];
float xij_j[3];
#endif
- // GFLOAT rhoe;
- // GFLOAT flux[5][3];
+ // float rhoe;
+ // float flux[5][3];
float dti, dtj, mindt;
float n_unit[3];
@@ -658,8 +658,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
// }
/* calculate the maximal signal velocity */
- vmax = sqrtf(const_hydro_gamma * Wi[4] / Wi[0]) +
- sqrtf(const_hydro_gamma * Wj[4] / Wj[0]);
+ vmax =
+ sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
(Wi[3] - Wj[3]) * dx[2];
if (dvdotdx > 0.) {
@@ -790,14 +790,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
#ifdef PER_FACE_LIMITER
float xij_i_norm;
- GFLOAT phi_i, phi_j;
- GFLOAT delta1, delta2;
- GFLOAT phiminus, phiplus;
- GFLOAT phimin, phimax;
- GFLOAT phibar;
+ float phi_i, phi_j;
+ float delta1, delta2;
+ float phiminus, phiplus;
+ float phimin, phimax;
+ float phibar;
/* free parameters, values from Hopkins */
- GFLOAT psi1 = 0.5, psi2 = 0.25;
- GFLOAT phi_mid0, phi_mid;
+ float psi1 = 0.5, psi2 = 0.25;
+ float phi_mid0, phi_mid;
for (k = 0; k < 10; k++) {
if (k < 5) {
@@ -877,13 +877,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Wi[2] * pi->primitives.gradients.v[2][1] +
Wi[3] * pi->primitives.gradients.v[2][2] +
pi->primitives.gradients.P[2] / Wi[0]);
- dWi[4] -= 0.5 * mindt *
- (Wi[1] * pi->primitives.gradients.P[0] +
- Wi[2] * pi->primitives.gradients.P[1] +
- Wi[3] * pi->primitives.gradients.P[2] +
- const_hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
- pi->primitives.gradients.v[1][1] +
- pi->primitives.gradients.v[2][2]));
+ dWi[4] -=
+ 0.5 * mindt * (Wi[1] * pi->primitives.gradients.P[0] +
+ Wi[2] * pi->primitives.gradients.P[1] +
+ Wi[3] * pi->primitives.gradients.P[2] +
+ hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
+ pi->primitives.gradients.v[1][1] +
+ pi->primitives.gradients.v[2][2]));
dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
Wj[2] * pj->primitives.gradients.rho[1] +
@@ -903,13 +903,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Wj[2] * pj->primitives.gradients.v[2][1] +
Wj[3] * pj->primitives.gradients.v[2][2] +
pj->primitives.gradients.P[2] / Wj[0]);
- dWj[4] -= 0.5 * mindt *
- (Wj[1] * pj->primitives.gradients.P[0] +
- Wj[2] * pj->primitives.gradients.P[1] +
- Wj[3] * pj->primitives.gradients.P[2] +
- const_hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
- pj->primitives.gradients.v[1][1] +
- pj->primitives.gradients.v[2][2]));
+ dWj[4] -=
+ 0.5 * mindt * (Wj[1] * pj->primitives.gradients.P[0] +
+ Wj[2] * pj->primitives.gradients.P[1] +
+ Wj[3] * pj->primitives.gradients.P[2] +
+ hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
+ pj->primitives.gradients.v[1][1] +
+ pj->primitives.gradients.v[2][2]));
// printf("WL: %g %g %g %g %g\n", Wi[0], Wi[1], Wi[2], Wi[3], Wi[4]);
// printf("WR: %g %g %g %g %g\n", Wj[0], Wj[1], Wj[2], Wj[3], Wj[4]);
@@ -969,7 +969,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
error("Negative density or pressure!\n");
}
- GFLOAT totflux[5];
+ float totflux[5];
riemann_solve_for_flux(Wi, Wj, n_unit, vij, totflux);
/* Update conserved variables */
@@ -1014,16 +1014,32 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
}
/* this corresponds to task_subtype_fluxes */
-__attribute__((always_inline)) INLINE static void runner_iact_hydro_loop3(
+__attribute__((always_inline)) INLINE static void runner_iact_force(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1);
}
/* this corresponds to task_subtype_fluxes */
-__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_hydro_loop3(float r2, float *dx, float hi, float hj,
- struct part *pi, struct part *pj) {
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0);
}
+
+/* PLACEHOLDERS */
+__attribute__((always_inline)) INLINE static void runner_iact_vec_density(
+ float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
+ struct part **pj) {}
+
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
+ struct part **pi, struct part **pj) {}
+
+__attribute__((always_inline)) INLINE static void runner_iact_vec_force(
+ float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
+ struct part **pj) {}
+
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
+ float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
+ struct part **pj) {}
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index 3c51653d994bd9f01864bcc24c6886eba25d1d05..b5738df8adca15ebd60da394aa88f5e17ff25362 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -1,6 +1,6 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Coypright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
*
* 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
@@ -17,77 +17,75 @@
*
******************************************************************************/
+#include "adiabatic_index.h"
+#include "io_properties.h"
+
/**
- * @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
+ * @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.
*/
-__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) {
+void hydro_read_particles(struct part* parts, struct io_props* list,
+ int* num_fields) {
+
+ *num_fields = 8;
- /* 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);
+ /* 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, primitives.P);
+ 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, OPTIONAL,
+ UNIT_CONV_DENSITY, parts, primitives.rho);
+}
+
+float convert_u(struct engine* e, struct part* p) {
+ return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
}
/**
- * @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
+ * @brief Specifies which particle fields to write to a dataset
*
+ * @param parts The particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
*/
-__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) {
+void hydro_write_particles(struct part* parts, struct io_props* list,
+ int* num_fields) {
- /* 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);
+ *num_fields = 8;
+
+ /* List what we want to write */
+ list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
+ parts, x);
+ list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts,
+ primitives.v);
+ list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts,
+ conserved.mass);
+ list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+ parts, h);
+ list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS,
+ parts, primitives.P, convert_u);
+ list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
+ UNIT_CONV_NO_UNITS, parts, id);
+ list[6] = io_make_output_field("Acceleration", FLOAT, 3,
+ UNIT_CONV_ACCELERATION, parts, a_hydro);
+ list[7] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
+ primitives.rho);
}
/**
@@ -96,25 +94,27 @@ __attribute__((always_inline)) INLINE static void hydro_write_particles(
*/
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)");
+ "Price (2008) without switch");
+ writeAttribute_f(h_grpsph, "Thermal Conductivity alpha",
+ const_conductivity_alpha);
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);
+ "Morris & Monaghan (1997), Rosswog, Davies, Thielemann & "
+ "Piran (2000) with additional Balsara (1995) switch");
+ writeAttribute_f(h_grpsph, "Viscosity alpha_min", const_viscosity_alpha_min);
+ writeAttribute_f(h_grpsph, "Viscosity alpha_max", const_viscosity_alpha_max);
+ writeAttribute_f(h_grpsph, "Viscosity beta", 2.f);
+ writeAttribute_f(h_grpsph, "Viscosity decay length", const_viscosity_length);
/* 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));
+ writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",
+ const_max_u_change);
}
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+int writeEntropyFlag() { return 0; }
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index 9e5f32f758248d1d1616f4556c81fc8e0b52e83b..1dea99a549a78f7902137a6fa0f323ab27d28450 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/Gizmo/hydro_part.h
@@ -19,23 +19,15 @@
*
******************************************************************************/
-/* 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];
+ /* Offset between current position and position at last tree rebuild. */
+ float x_diff[3];
/* Velocity at the last full step. */
float v_full[3];
- /* Entropy at the half-step. */
- float u_hdt;
-
/* Old density. */
float omega;
@@ -47,17 +39,12 @@ struct part {
/* Particle position. */
double x[3];
- /* Particle velocity. */
+ /* Particle predicted 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;
@@ -67,76 +54,86 @@ struct part {
/* Particle time of end of time-step. */
int ti_end;
- /* The primitive hydrodynamical variables */
+ /* The primitive hydrodynamical variables. */
struct {
- /* fluid velocity */
- GFLOAT v[3];
+ /* Fluid velocity. */
+ float v[3];
- /* density */
- GFLOAT rho;
+ /* Density. */
+ float rho;
- /* pressure */
- GFLOAT P;
+ /* Pressure. */
+ float P;
+ /* Gradients of the primitive variables. */
struct {
- GFLOAT rho[3];
+ /* Density gradients. */
+ float rho[3];
- GFLOAT v[3][3];
+ /* Fluid velocity gradients. */
+ float v[3][3];
- GFLOAT P[3];
+ /* Pressure gradients. */
+ float P[3];
} gradients;
+ /* Quantities needed by the slope limiter. */
struct {
- /* extreme values among the neighbours */
- GFLOAT rho[2];
+ /* Extreme values of the density among the neighbours. */
+ float rho[2];
- GFLOAT v[3][2];
+ /* Extreme values of the fluid velocity among the neighbours. */
+ float v[3][2];
- GFLOAT P[2];
+ /* Extreme values of the pressure among the neighbours. */
+ float P[2];
- /* maximal distance to all neighbouring faces */
+ /* Maximal distance to all neighbouring faces. */
float maxr;
} limiter;
} primitives;
- /* The conserved hydrodynamical variables */
+ /* The conserved hydrodynamical variables. */
struct {
- /* fluid momentum */
- GFLOAT momentum[3];
+ /* Fluid momentum. */
+ float momentum[3];
- /* fluid mass */
- GFLOAT mass;
+ /* Fluid mass (this field already exists outside of this struct as well). */
+ float mass;
- /* fluid energy */
- GFLOAT energy;
+ /* Fluid total energy (= internal energy + fluid kinetic energy). */
+ float energy;
} conserved;
- /* Geometrical quantities used for hydro */
+ /* Geometrical quantities used for hydro. */
struct {
- /* volume of the particle */
+ /* Volume of the particle. */
float volume;
- /* gradient matrix */
+ /* Geometrical shear matrix used to calculate second order accurate
+ gradients */
float matrix_E[3][3];
} geometry;
+ /* Variables used for timestep calculation (currently not used). */
struct {
+ /* Maximum fluid velocity among all neighbours. */
float vmax;
} timestepvars;
- /* Quantities used during the density loop */
+ /* Quantities used during the volume (=density) loop. */
struct {
/* Particle velocity divergence. */
@@ -153,10 +150,21 @@ struct part {
} density;
+ /* Particle mass (this field is also part of the conserved quantities...). */
+ float mass;
+
/* Particle ID. */
unsigned long long id;
/* Associated gravitas. */
struct gpart *gpart;
+ /* Variables needed for the code to compile (should be removed/replaced). */
+ float rho;
+ struct {
+ float h_dt;
+ float v_sig;
+ } force;
+ float rho_dh;
+
} __attribute__((aligned(part_align)));
diff --git a/src/hydro_io.h b/src/hydro_io.h
index 30d663f647c9b763e9b19177e9ba8ef374855768..f0a619b90b774574c434007b1c01a0e55e75e464 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -28,6 +28,8 @@
#include "./hydro/Gadget2/hydro_io.h"
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_io.h"
+#elif defined(GIZMO_SPH)
+#include "./hydro/Gizmo/hydro_io.h"
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/part.h b/src/part.h
index efca7b6b5bef49f20df1e2c45b30f65ecbbf4960..e106be3962dbac9fe9c985e2e9696350b2f138c3 100644
--- a/src/part.h
+++ b/src/part.h
@@ -45,6 +45,8 @@
#include "./hydro/Gadget2/hydro_part.h"
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_part.h"
+#elif defined(GIZMO_SPH)
+#include "./hydro/Gizmo/hydro_part.h"
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/riemann/riemann_exact.h b/src/riemann/riemann_exact.h
index a2f3c30fb1daf5d53bf35abe4ca7e73eafba6018..ab5c95209dc0ed633c4ebeb80716802d1b3305df 100644
--- a/src/riemann/riemann_exact.h
+++ b/src/riemann/riemann_exact.h
@@ -50,12 +50,11 @@
* @param W The left or right state vector
* @param a The left or right sound speed
*/
-__attribute__((always_inline)) INLINE static GFLOAT riemann_fb(GFLOAT p,
- GFLOAT* W,
- GFLOAT a) {
+__attribute__((always_inline)) INLINE static float riemann_fb(float p, float* W,
+ float a) {
- GFLOAT fval = 0.;
- GFLOAT A, B;
+ float fval = 0.;
+ float A, B;
if (p > W[4]) {
A = const_riemann_tdgp1 / W[0];
B = const_riemann_gm1dgp1 * W[4];
@@ -78,9 +77,8 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_fb(GFLOAT p,
* @param aL The left sound speed
* @param aR The right sound speed
*/
-__attribute__((always_inline)) INLINE static GFLOAT riemann_f(
- GFLOAT p, GFLOAT* WL, GFLOAT* WR, GFLOAT vL, GFLOAT vR, GFLOAT aL,
- GFLOAT aR) {
+__attribute__((always_inline)) INLINE static float riemann_f(
+ float p, float* WL, float* WR, float vL, float vR, float aL, float aR) {
return riemann_fb(p, WL, aL) + riemann_fb(p, WR, aR) + (vR - vL);
}
@@ -92,12 +90,12 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_f(
* @param W The left or right state vector
* @param a The left or right sound speed
*/
-__attribute__((always_inline)) INLINE static GFLOAT riemann_fprimeb(GFLOAT p,
- GFLOAT* W,
- GFLOAT a) {
+__attribute__((always_inline)) INLINE static float riemann_fprimeb(float p,
+ float* W,
+ float a) {
- GFLOAT fval = 0.;
- GFLOAT A, B;
+ float fval = 0.;
+ float A, B;
if (p > W[4]) {
A = const_riemann_tdgp1 / W[0];
B = const_riemann_gm1dgp1 * W[4];
@@ -117,8 +115,8 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_fprimeb(GFLOAT p,
* @param aL The left sound speed
* @param aR The right sound speed
*/
-__attribute__((always_inline)) INLINE static GFLOAT riemann_fprime(
- GFLOAT p, GFLOAT* WL, GFLOAT* WR, GFLOAT aL, GFLOAT aR) {
+__attribute__((always_inline)) INLINE static float riemann_fprime(
+ float p, float* WL, float* WR, float aL, float aR) {
return riemann_fprimeb(p, WL, aL) + riemann_fprimeb(p, WR, aR);
}
@@ -129,10 +127,10 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_fprime(
* @param p The current guess for the pressure
* @param W The left or right state vector
*/
-__attribute__((always_inline)) INLINE static GFLOAT riemann_gb(GFLOAT p,
- GFLOAT* W) {
+__attribute__((always_inline)) INLINE static float riemann_gb(float p,
+ float* W) {
- GFLOAT A, B;
+ float A, B;
A = const_riemann_tdgp1 / W[0];
B = const_riemann_gm1dgp1 * W[4];
return sqrtf(A / (p + B));
@@ -151,11 +149,11 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_gb(GFLOAT p,
* @param aL The left sound speed
* @param aR The right sound speed
*/
-__attribute__((always_inline)) INLINE static GFLOAT riemann_guess_p(
- GFLOAT* WL, GFLOAT* WR, GFLOAT vL, GFLOAT vR, GFLOAT aL, GFLOAT aR) {
+__attribute__((always_inline)) INLINE static float riemann_guess_p(
+ float* WL, float* WR, float vL, float vR, float aL, float aR) {
- GFLOAT pguess, pmin, pmax, qmax;
- GFLOAT ppv;
+ float pguess, pmin, pmax, qmax;
+ float ppv;
pmin = fminf(WL[4], WR[4]);
pmax = fmaxf(WL[4], WR[4]);
@@ -202,14 +200,14 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_guess_p(
* @param aL The left sound speed
* @param aR The right sound speed
*/
-__attribute__((always_inline)) INLINE static GFLOAT riemann_solve_brent(
- GFLOAT lower_limit, GFLOAT upper_limit, GFLOAT lowf, GFLOAT upf,
- GFLOAT error_tol, GFLOAT* WL, GFLOAT* WR, GFLOAT vL, GFLOAT vR, GFLOAT aL,
- GFLOAT aR) {
-
- GFLOAT a, b, c, d, s;
- GFLOAT fa, fb, fc, fs;
- GFLOAT tmp, tmp2;
+__attribute__((always_inline)) INLINE static float riemann_solve_brent(
+ float lower_limit, float upper_limit, float lowf, float upf,
+ float error_tol, float* WL, float* WR, float vL, float vR, float aL,
+ float aR) {
+
+ float a, b, c, d, s;
+ float fa, fb, fc, fs;
+ float tmp, tmp2;
int mflag;
int i;
@@ -309,11 +307,11 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_solve_brent(
* @param n_unit Normal vector of the interface
*/
__attribute__((always_inline)) INLINE static void riemann_solve_vacuum(
- GFLOAT* WL, GFLOAT* WR, GFLOAT vL, GFLOAT vR, GFLOAT aL, GFLOAT aR,
- GFLOAT* Whalf, float* n_unit) {
+ float* WL, float* WR, float vL, float vR, float aL, float aR, float* Whalf,
+ float* n_unit) {
- GFLOAT SL, SR;
- GFLOAT vhalf;
+ float SL, SR;
+ float vhalf;
if (!WR[0] && !WL[0]) {
/* if both states are vacuum, the solution is also vacuum */
@@ -456,20 +454,20 @@ __attribute__((always_inline)) INLINE static void riemann_solve_vacuum(
* @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) {
+ float* WL, float* WR, float* Whalf, float* n_unit) {
/* velocity of the left and right state in a frame aligned with n_unit */
- GFLOAT vL, vR, vhalf;
+ float vL, vR, vhalf;
/* sound speeds */
- GFLOAT aL, aR;
+ float aL, aR;
/* variables used for finding pstar */
- GFLOAT p, pguess, fp, fpguess;
+ float p, pguess, fp, fpguess;
/* variables used for sampling the solution */
- GFLOAT u;
- GFLOAT pdpR, SR;
- GFLOAT SHR, STR;
- GFLOAT pdpL, SL;
- GFLOAT SHL, STL;
+ float u;
+ float pdpR, SR;
+ float SHR, STR;
+ float pdpL, SL;
+ float SHL, STL;
int errorFlag = 0;
/* sanity checks */
@@ -661,10 +659,10 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
}
__attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
- GFLOAT* Wi, GFLOAT* Wj, float* n_unit, float* vij, GFLOAT* totflux) {
+ float* Wi, float* Wj, float* n_unit, float* vij, float* totflux) {
- GFLOAT Whalf[5];
- GFLOAT flux[5][3];
+ float Whalf[5];
+ float flux[5][3];
float vtot[3];
float rhoe;
diff --git a/src/riemann/riemann_hllc.h b/src/riemann/riemann_hllc.h
index 6c583f6410f53ed64d630082926d816129768fab..af83bed6c3aabfb923ed465cebb0a54461702094 100644
--- a/src/riemann/riemann_hllc.h
+++ b/src/riemann/riemann_hllc.h
@@ -20,13 +20,15 @@
#ifndef SWIFT_RIEMANN_HLLC_H
#define SWIFT_RIEMANN_HLLC_H
+#include "adiabatic_index.h"
+
__attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
- GFLOAT *WL, GFLOAT *WR, float *n, float *vij, GFLOAT *totflux) {
+ float *WL, float *WR, float *n, float *vij, float *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];
+ float uL, uR, aL, aR;
+ float rhobar, abar, pPVRS, pstar, qL, qR, SL, SR, Sstar;
+ float v2, eL, eR;
+ float UstarL[5], UstarR[5];
/* Handle pure vacuum */
if (!WL[0] && !WR[0]) {
@@ -41,14 +43,14 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
/* 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]);
+ aL = sqrtf(hydro_gamma * WL[4] / WL[0]);
+ aR = sqrtf(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.) <
+ if (2. * aL / hydro_gamma_minus_one + 2. * aR / hydro_gamma_minus_one <
fabs(uL - uR)) {
error("Vacuum not yet supported");
}
@@ -64,14 +66,12 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
qL = 1.;
if (pstar > WL[4]) {
qL = sqrtf(1. +
- 0.5 * (const_hydro_gamma + 1.) / const_hydro_gamma *
- (pstar / WL[4] - 1.));
+ 0.5 * (hydro_gamma + 1.) / 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.));
+ 0.5 * (hydro_gamma + 1.) / hydro_gamma * (pstar / WR[4] - 1.));
}
SL = uL - aL * qL;
SR = uR + aR * qR;
@@ -88,7 +88,7 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
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;
+ eL = WL[4] / hydro_gamma_minus_one / WL[0] + 0.5 * v2;
totflux[4] = WL[0] * eL * uL + WL[4] * uL;
if (SL < 0.) {
/* add flux FstarL */
@@ -118,7 +118,7 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
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;
+ eR = WR[4] / hydro_gamma_minus_one / WR[0] + 0.5 * v2;
totflux[4] = WR[0] * eR * uR + WR[4] * uR;
if (SR > 0.) {
/* add flux FstarR */
diff --git a/src/riemann/riemann_trrs.h b/src/riemann/riemann_trrs.h
index efdbfb59877c09a59d535a4785ad74620c0f3651..4c031997d821790e0fd318004722b75d9e8ec8c2 100644
--- a/src/riemann/riemann_trrs.h
+++ b/src/riemann/riemann_trrs.h
@@ -50,14 +50,14 @@
* @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;
+ float* WL, float* WR, float* Whalf, float* n_unit) {
+ float aL, aR;
+ float PLR;
+ float vL, vR;
+ float ustar, pstar;
+ float vhalf;
+ float pdpR, SHR, STR;
+ float 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];
diff --git a/src/timestep.h b/src/timestep.h
index 99747484b2fad2d1dfadad749232f77848e56f35..569120cf9cf989b633da35529f7693c8f62a1910 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -108,7 +108,7 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
e->external_potential, e->physical_constants, p->gpart);
/* const float new_dt_self = */
/* gravity_compute_timestep_self(e->physical_constants, p->gpart); */
- const float new_dt_self = FLT_MAX; // MATTHIEU
+ const float new_dt_self = FLT_MAX; // MATTHIEU
new_dt_grav = fminf(new_dt_external, new_dt_self);
}
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 4a045307a7d731a71dabaa96ed39b66d013d4627..db228a3939621e383a2c63659113b104fc2c1bdb 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -99,6 +99,8 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(MINIMAL_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
+#elif defined(GIZMO_SPH)
+ part->primitives.P = pressure;
#else
error("Need to define pressure here !");
#endif