Merge branch 'comoving_energy_sph' into 'master'

Comoving energy sph

See merge request !597

examples/ZeldovichPancake_3D/plotSolution.py

@@ -69,6 +69,7 @@ scheme = sim["/HydroScheme"].attrs["Scheme"]
 kernel = sim["/HydroScheme"].attrs["Kernel function"]
 neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
 eta = sim["/HydroScheme"].attrs["Kernel eta"]
+alpha = sim["/HydroScheme"].attrs["Alpha viscosity"]
 git = sim["Code"].attrs["Git Revision"]
 
 # Cosmological parameters
@@ -178,8 +179,8 @@ ylabel("${\\rm{Temperature}}~T$", labelpad=0)
 # Information -------------------------------------
 subplot(236, frameon=False)
 
-text(-0.49, 0.9, "Zeldovich pancake with  $\\gamma=%.3f$ in 1D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
-text(-0.49, 0.8, "$z={0:.2f}$".format(redshift))
+text(-0.49, 0.9, "Zeldovich pancake at z=%.2f "%(redshift), fontsize=10)
+text(-0.49, 0.8, "adiabatic index $\\gamma=%.2f$, viscosity $\\alpha=%.2f$"%(gas_gamma, alpha), fontsize=10)
 plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
 text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
 text(-0.49, 0.4, scheme, fontsize=10)

examples/ZeldovichPancake_3D/zeldovichPancake.yml

@@ -24,7 +24,8 @@ Snapshots:
   basename:            zeldovichPancake # Common part of the name of output files
   time_first:          0.       # Time of the first output (in internal units)
   delta_time:          1.04     # Time difference between consecutive outputs (in internal units)
-  scale_factor_first: 0.00991
+  scale_factor_first:  0.00991
+  compression:         4
 
 # Parameters governing the conserved quantities statistics
 Statistics:

src/hydro/Minimal/hydro.h

@@ -570,6 +570,12 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
     struct part *restrict p, struct xpart *restrict xp,
     const struct cosmology *cosmo) {
 
+  /* 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 *= factor;
+
   /* Compute the pressure */
   const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
 

src/hydro/Minimal/hydro_iact.h

@@ -172,10 +172,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];
+
+  /* Add Hubble flow */
+  const float dvdr_Hubble = dvdr + a2_Hubble * r2;
 
   /* Are the particles 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 */
@@ -285,10 +288,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];
+
+  /* Add Hubble flow */
+  const float dvdr_Hubble = dvdr + a2_Hubble * r2;
 
   /* Are the particles 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_properties.c

@@ -215,6 +215,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
                        p->initial_internal_energy);
   io_write_attribute_f(h_grpsph, "Hydrogen mass fraction",
                        p->hydrogen_mass_fraction);
+  io_write_attribute_f(h_grpsph, "Alpha viscosity", const_viscosity_alpha);
 }
 #endif
 

theory/Cosmology/coordinates.tex

@@ -88,13 +88,13 @@ gravitational terms. SPH flavours that evolve the internal energy $u$ instead of
 entropy require the additional equation of motion describing the evolution of
 $u'$:
 \begin{equation}
-  \dot{u}_i' = \frac{P_i'}{\rho_i'^2}\left[3H\rho_i' + \frac{1}{a^2}f_i'\sum_jm_j\left(\mathbf{v}_i' -
-    \mathbf{v}_j'\right)\cdot\mathbf{\nabla}_i'W_{ij}'(h_i)\right],
+  \dot{u}_i' = \frac{1}{a^2}\frac{P_i'}{\rho_i'^2} f_i'\sum_jm_j\left(\mathbf{v}_i' -
+    \mathbf{v}_j'\right)\cdot\mathbf{\nabla}_i'W_{ij}'(h_i).
   \label{eq:cosmo_eom_u}
 \end{equation}
-where the first term in the brackets accounts for the change in energy
-due to the expansion of the Universe. The scale-factors appearing in
-the equations are later absorbed in the time-integration operators
+
+In all these cases, the scale-factors appearing in the equations are
+later absorbed in the time-integration operators
 (Sec.~\ref{ssec:operators}) such that the RHS of the equations of
 motions is identical for the primed quantities to the ones obtained in
 the non-cosmological case for the physical quantities.