Skip to content
Snippets Groups Projects
Commit 63a3cf62 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Comoving energy sph

parent cfd6f2c0
No related branches found
No related tags found
No related merge requests found
...@@ -69,6 +69,7 @@ scheme = sim["/HydroScheme"].attrs["Scheme"] ...@@ -69,6 +69,7 @@ scheme = sim["/HydroScheme"].attrs["Scheme"]
kernel = sim["/HydroScheme"].attrs["Kernel function"] kernel = sim["/HydroScheme"].attrs["Kernel function"]
neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"] neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
eta = sim["/HydroScheme"].attrs["Kernel eta"] eta = sim["/HydroScheme"].attrs["Kernel eta"]
alpha = sim["/HydroScheme"].attrs["Alpha viscosity"]
git = sim["Code"].attrs["Git Revision"] git = sim["Code"].attrs["Git Revision"]
# Cosmological parameters # Cosmological parameters
...@@ -178,8 +179,8 @@ ylabel("${\\rm{Temperature}}~T$", labelpad=0) ...@@ -178,8 +179,8 @@ ylabel("${\\rm{Temperature}}~T$", labelpad=0)
# Information ------------------------------------- # Information -------------------------------------
subplot(236, frameon=False) 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.9, "Zeldovich pancake at z=%.2f "%(redshift), fontsize=10)
text(-0.49, 0.8, "$z={0:.2f}$".format(redshift)) 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) 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.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
text(-0.49, 0.4, scheme, fontsize=10) text(-0.49, 0.4, scheme, fontsize=10)
......
...@@ -24,7 +24,8 @@ Snapshots: ...@@ -24,7 +24,8 @@ Snapshots:
basename: zeldovichPancake # Common part of the name of output files basename: zeldovichPancake # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units) time_first: 0. # Time of the first output (in internal units)
delta_time: 1.04 # Time difference between consecutive outputs (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 # Parameters governing the conserved quantities statistics
Statistics: Statistics:
......
...@@ -570,6 +570,12 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( ...@@ -570,6 +570,12 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p, struct xpart *restrict xp, struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) { 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 */ /* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
......
...@@ -172,10 +172,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -172,10 +172,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Compute dv dot r. */ /* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] + (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 ? */ /* 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 */ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */ /* Compute sound speeds and signal velocity */
...@@ -285,10 +288,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -285,10 +288,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Compute dv dot r. */ /* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] + (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 ? */ /* 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 */ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */ /* Compute sound speeds and signal velocity */
......
...@@ -215,6 +215,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) { ...@@ -215,6 +215,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
p->initial_internal_energy); p->initial_internal_energy);
io_write_attribute_f(h_grpsph, "Hydrogen mass fraction", io_write_attribute_f(h_grpsph, "Hydrogen mass fraction",
p->hydrogen_mass_fraction); p->hydrogen_mass_fraction);
io_write_attribute_f(h_grpsph, "Alpha viscosity", const_viscosity_alpha);
} }
#endif #endif
......
...@@ -88,13 +88,13 @@ gravitational terms. SPH flavours that evolve the internal energy $u$ instead of ...@@ -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 entropy require the additional equation of motion describing the evolution of
$u'$: $u'$:
\begin{equation} \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' - \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)\right], \mathbf{v}_j'\right)\cdot\mathbf{\nabla}_i'W_{ij}'(h_i).
\label{eq:cosmo_eom_u} \label{eq:cosmo_eom_u}
\end{equation} \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 In all these cases, the scale-factors appearing in the equations are
the equations are later absorbed in the time-integration operators later absorbed in the time-integration operators
(Sec.~\ref{ssec:operators}) such that the RHS of the equations of (Sec.~\ref{ssec:operators}) such that the RHS of the equations of
motions is identical for the primed quantities to the ones obtained in motions is identical for the primed quantities to the ones obtained in
the non-cosmological case for the physical quantities. the non-cosmological case for the physical quantities.
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment