Make sure the fluid quantities are converted correctly from ICs to internal stuff in GIZMO.

src/hydro/Gizmo/hydro.h

@@ -124,14 +124,12 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 
   const float ihdim = pow_dimension(ih);
 
-  float volume;
-  float m, momentum[3], energy;
-
   /* Final operation on the geometry. */
   /* we multiply with the smoothing kernel normalization ih3 and calculate the
    * volume */
-  volume = ihdim * (p->geometry.volume + kernel_root);
-  p->geometry.volume = volume = 1. / volume;
+  const float volume = 1.f / (ihdim * (p->geometry.volume + kernel_root));
+  p->geometry.volume = volume;
+
   /* we multiply with the smoothing kernel normalization */
   p->geometry.matrix_E[0][0] = ihdim * p->geometry.matrix_E[0][0];
   p->geometry.matrix_E[0][1] = ihdim * p->geometry.matrix_E[0][1];
@@ -149,8 +147,9 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 
   /* compute primitive variables */
   /* eqns (3)-(5) */
-  m = p->conserved.mass;
-  if (m) {
+  const float m = p->conserved.mass;
+  if (m > 0.f) {
+    float momentum[3];
     momentum[0] = p->conserved.momentum[0];
     momentum[1] = p->conserved.momentum[1];
     momentum[2] = p->conserved.momentum[2];
@@ -158,7 +157,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
     p->primitives.v[0] = momentum[0] / m;
     p->primitives.v[1] = momentum[1] / m;
     p->primitives.v[2] = momentum[2] / m;
-    energy = p->conserved.energy;
+    const float energy = p->conserved.energy;
     p->primitives.P = hydro_gamma_minus_one * energy / volume;
   }
 }
@@ -254,18 +253,17 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
 __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
     struct part* p) {
 
-  float momentum[3];
   const float volume = p->geometry.volume;
   const float m = p->conserved.mass;
   p->primitives.rho = m / volume;
 
-  /* P actually contains internal energy at this point */
-  p->primitives.P *= hydro_gamma_minus_one * p->primitives.rho;
+  p->conserved.momentum[0] = m * p->primitives.v[0];
+  p->conserved.momentum[1] = m * p->primitives.v[1];
+  p->conserved.momentum[2] = m * p->primitives.v[2];
+  p->conserved.energy *= m;
 
-  p->conserved.momentum[0] = momentum[0] = m * p->primitives.v[0];
-  p->conserved.momentum[1] = momentum[1] = m * p->primitives.v[1];
-  p->conserved.momentum[2] = momentum[2] = m * p->primitives.v[2];
-  p->conserved.energy = p->primitives.P / hydro_gamma_minus_one * volume;
+  p->primitives.P =
+      hydro_gamma_minus_one * p->conserved.energy * p->primitives.rho;
 }
 
 /**

src/hydro/Gizmo/hydro_io.h

@@ -44,9 +44,9 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
                                 parts, conserved.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[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
+                                UNIT_CONV_ENERGY_PER_UNIT_MASS, parts,
+                                conserved.energy);
   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,

src/hydro/Gizmo/hydro_part.h

@@ -103,7 +103,7 @@ struct part {
     /* Fluid momentum. */
     float momentum[3];
 
-    /* Fluid mass (this field already exists outside of this struct as well). */
+    /* Fluid mass */
     float mass;
 
     /* Fluid thermal energy (not per unit mass!). */

[Bert Vandenbroucke]

This should not be done until after the initial pressure calculation.

[Matthieu Schaller]

Well, you were doing it. I have only removed the temporary variable momentum[3] here that was not used elsewhere.

[Bert Vandenbroucke]

No, in the last version that I have the internal energy (not per unit mass) is calculated from the pressure. The pressure initially contains the internal energy per unit mass and is then converted into the pressure.

If u is the internal energy per unit mass, then P = (gamma-1) * rhou; conserved.energy = P / (gamma-1) * volume = um.

What you have now is conserved.energy = um (which is correct), and P = (gamma-1)rhomu (which is not correct). If you put the p->conserved.energy *= m below the computation of P, it will be correct.

[Bert Vandenbroucke]

Not that it really matters though. The pressure you compute here is never actually used and is recalculated after the second density loop during the first real time step.

[Matthieu Schaller]

Ok, I see what you mean now. Thought you were referring to the momentum lines. Indeed, these values are unused but I'll make the change nevertheless in order to be clean. Thanks !