Merge branch 'cosmology-pressure-energy' into 'master'

Cosmology pressure energy

See merge request !611

examples/ConstantCosmoVolume/makeIC.py

@@ -65,7 +65,7 @@ glass.close()
 # Total number of particles
 numPart = size(pos)/3
 if numPart != numPart_1D**3:
-  print "Non-matching glass file"
+  print("Non-matching glass file")
 numPart = numPart_1D**3
 
 # Set box size and interparticle distance

examples/ZeldovichPancake_3D/plotSolution.py

@@ -83,7 +83,12 @@ S = sim["/PartType0/Entropy"][:]
 P = sim["/PartType0/Pressure"][:]
 rho = sim["/PartType0/Density"][:]
 m = sim["/PartType0/Masses"][:]
-phi = sim["/PartType0/Potential"][:]
+try:
+    phi = sim["/PartType0/Potential"][:]
+except KeyError:
+    # We didn't write the potential, try to go on without
+    print("Couldn't find potential in your output file")
+    phi = np.zeros_like(m)
 
 x -= 0.5 * boxSize
 
@@ -97,7 +102,7 @@ if os.path.exists(filename_g):
     rho_g = sim_g["/PartType0/Density"][:]
     phi_g = sim_g["/PartType0/Potential"][:]
     a_g = sim_g["/Header"].attrs["Time"]
-    print "Gadget Scale-factor:", a_g, "redshift:", 1/a_g - 1.
+    print("Gadget Scale-factor:", a_g, "redshift:", 1/a_g - 1.)
     
     x_g -= 0.5 * boxSize
 else:
@@ -168,7 +173,7 @@ u /= a**(3 * (gas_gamma - 1.))
 u_g /= a**(3 * (gas_gamma - 1.))
 T = (gas_gamma - 1.) * u * mH_in_kg / k_in_J_K
 T_g = (gas_gamma - 1.) * u_g * mH_in_kg / k_in_J_K
-print "z = {0:.2f}, T_avg = {1:.2f}".format(redshift, T.mean())
+print("z = {0:.2f}, T_avg = {1:.2f}".format(redshift, T.mean()))
 if np.size(x_g) > 1:
     plot(x_g, T_g, 's', color='g', alpha=0.8, lw=1.2, ms=4)
 plot(x, T, '.', color='r', ms=4.0)

src/hydro/PressureEnergy/hydro.h

@@ -580,12 +580,16 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
  * @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_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
  */
 __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_hydro, float dt_kick_corr,
-    const struct cosmology *cosmo,
-    const struct hydro_props *restrict hydro_properties) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
 
   /* Do not decrease the energy by more than a factor of 2*/
   if (dt_therm > 0. && p->u_dt * dt_therm < -0.5f * xp->u_full) {
@@ -593,6 +597,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
   }
   xp->u_full += p->u_dt * dt_therm;
 
+  /* Apply the minimal energy limit */
+  const float min_energy =
+      hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+  if (xp->u_full < min_energy) {
+    xp->u_full = min_energy;
+    p->u_dt = 0.f;
+  }
+
   /* Compute the sound speed */
   const float soundspeed = hydro_get_comoving_soundspeed(p);
 
@@ -609,10 +621,31 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
  *
  * @param p The particle to act upon
  * @param xp The extended particle to act upon
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme.
  */
 __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 cosmology *cosmo, const struct hydro_props *hydro_props) {
+
+  /* Convert the physcial internal energy to the comoving one. */
+  /* u' = a^(3(g-1)) u */
+  const float factor = 1.f / cosmo->a_factor_internal_energy;
+  p->u *= factor;
+  xp->u_full = p->u;
+
+  /* Apply the minimal energy limit */
+  const float min_energy =
+      hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
+  if (xp->u_full < min_energy) {
+    xp->u_full = min_energy;
+    p->u = min_energy;
+    p->u_dt = 0.f;
+  }
+
+  /* Note that unlike Minimal the pressure and sound speed cannot be calculated
+   * here because they are smoothed properties in this scheme. */
+}
 
 /**
  * @brief Initialises the particles for the first time

src/hydro/PressureEnergy/hydro_iact.h

@@ -232,10 +232,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   /* Compute dv dot r. */
   const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
                      (pi->v[1] - pj->v[1]) * dx[1] +
-                     (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
+                     (pi->v[2] - pj->v[2]) * dx[2];
+
+  /* Includes the hubble flow term; not used for du/dt */
+  const float dvdr_Hubble = dvdr + a2_Hubble * r2;
 
   /* Are the part*icles moving towards each others ? */
-  const float omega_ij = min(dvdr, 0.f);
+  const float omega_ij = min(dvdr_Hubble, 0.f);
   const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
 
   /* Compute sound speeds and signal velocity */
@@ -357,10 +360,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   /* Compute dv dot r. */
   const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
                      (pi->v[1] - pj->v[1]) * dx[1] +
-                     (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
+                     (pi->v[2] - pj->v[2]) * dx[2];
+
+  /* Includes the hubble flow term; not used for du/dt */
+  const float dvdr_Hubble = dvdr + a2_Hubble * r2;
 
   /* Are the part*icles moving towards each others ? */
-  const float omega_ij = min(dvdr, 0.f);
+  const float omega_ij = min(dvdr_Hubble, 0.f);
   const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
 
   /* Compute sound speeds and signal velocity */

src/hydro/PressureEnergy/hydro_io.h

@@ -136,6 +136,15 @@ INLINE static void convert_part_vel(const struct engine* e,
   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
  *
@@ -148,7 +157,7 @@ INLINE static void hydro_write_particles(const struct part* parts,
                                          struct io_props* list,
                                          int* num_fields) {
 
-  *num_fields = 9;
+  *num_fields = 10;
 
   /* List what we want to write */
   list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
@@ -172,6 +181,9 @@ INLINE static void hydro_write_particles(const struct part* parts,
   list[8] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
                                               UNIT_CONV_ENTROPY_PER_UNIT_MASS,
                                               parts, xparts, convert_S);
+  list[9] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+		                              UNIT_CONV_POTENTIAL, parts,
+					      xparts, convert_part_potential);
 }
 
 /**