Added option to use densities rather than masses from the IC for GizmoMFV.

src/const.h

@@ -96,6 +96,9 @@
    information. */
 #define GIZMO_FLAG_VARIABLE
 
+/* Initialize densities instead of masses to reduce particle noise. */
+//#define GIZMO_USE_IC_DENSITY
+
 /* Types of gradients to use for SHADOWFAX_SPH */
 /* If no option is chosen, no gradients are used (first order scheme) */
 #define SHADOWFAX_GRADIENTS

src/hydro/GizmoMFV/hydro.h

@@ -123,6 +123,17 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
 __attribute__((always_inline)) INLINE static void hydro_first_init_part(
     struct part* p, struct xpart* xp) {
 
+#ifdef GIZMO_USE_IC_DENSITY
+  /* we temporarily store the primitive variables in the flux variables.
+     This is a hack: we will restore them in hydro_convert_quantities after
+     the volumes have been computed. */
+  p->conserved.flux.mass = p->primitives.rho;
+  p->conserved.flux.momentum[0] = p->v[0];
+  p->conserved.flux.momentum[1] = p->v[1];
+  p->conserved.flux.momentum[2] = p->v[2];
+  p->conserved.flux.energy = p->conserved.energy;
+#endif
+
   const float mass = p->conserved.mass;
 
   p->time_bin = 0;
@@ -372,7 +383,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
      this is why we divide by the volume, and not by the density */
   p->primitives.P = hydro_gamma_minus_one * energy / volume;
   p->primitives.A = p->conserved.entropy / m;
-  p->primitives.P = p->primitives.A * pow_gamma(p->primitives.rho);
 #endif
 
   /* sanity checks */
@@ -580,6 +590,39 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
     struct part* p, struct xpart* xp, const struct cosmology* cosmo,
     const struct hydro_props* hydro_props) {
 
+#ifdef GIZMO_USE_IC_DENSITY
+  /* read the temporarily stored initial density and internal energy from the
+     snapshot */
+  const float rho = p->conserved.flux.mass;
+  const float v[3] = {p->conserved.flux.momentum[0],
+                      p->conserved.flux.momentum[1],
+                      p->conserved.flux.momentum[2]};
+  const float u = p->conserved.flux.energy;
+
+  /* get the particle volume */
+  const float V = p->geometry.volume;
+
+  /* now compute all derived quantities */
+  p->conserved.mass = rho * V;
+  p->conserved.momentum[0] = p->conserved.mass * v[0];
+  p->conserved.momentum[1] = p->conserved.mass * v[1];
+  p->conserved.momentum[2] = p->conserved.mass * v[2];
+  p->conserved.energy = p->conserved.mass * u;
+
+  p->primitives.rho = rho;
+  p->primitives.v[0] = v[0];
+  p->primitives.v[1] = v[1];
+  p->primitives.v[2] = v[2];
+  p->primitives.P = gas_pressure_from_internal_energy(rho, u);
+  p->primitives.A = p->primitives.P * pow_minus_gamma(rho);
+
+  p->v[0] = v[0];
+  p->v[1] = v[1];
+  p->v[2] = v[2];
+
+  hydro_velocities_init(p, xp);
+#endif
+
   p->conserved.energy /= cosmo->a_factor_internal_energy;
 
   /* initialize the entropy */

src/hydro/GizmoMFV/hydro_io.h

@@ -62,8 +62,13 @@ INLINE static void hydro_read_particles(struct part* parts,
                                 UNIT_CONV_NO_UNITS, parts, id);
   list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
                                 UNIT_CONV_ACCELERATION, parts, a_hydro);
+#ifdef GIZMO_USE_IC_DENSITY
+  list[7] = io_make_input_field("Density", FLOAT, 1, COMPULSORY,
+                                UNIT_CONV_DENSITY, parts, primitives.rho);
+#else
   list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
                                 UNIT_CONV_DENSITY, parts, primitives.rho);
+#endif
 }
 
 /**