### Merge branch 'cleanup_rho_dh' into 'master'

Cleanup rho dh and use of the EoS functions

A few more rather unimportant cosmetic changes:

- moved rho_dh to density.rho_dh everywhere.
- added a force.f term to contain the transformed value of rho_dh.
- When computing the sound speed, I am now using the functions defined in equation_of_state.h everywhere.

The particles are the same size in memory as before despite the change. I had to update the tolerances for the tests as rho_dh changed meaning.

See merge request !255
parents 61111a27 54649eb4
 ... ... @@ -126,6 +126,21 @@ gas_soundspeed_from_internal_energy(float density, float u) { return sqrtf(u * hydro_gamma * hydro_gamma_minus_one); } /** * @brief Returns the sound speed given density and pressure * * Computes \f$c = \sqrt{\frac{\gamma P}{\rho} }\f$. * * @param density The density \f$\rho\f$ * @param P The pressure \f$P\f$ */ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure( float density, float P) { const float density_inv = 1.f / density; return sqrtf(hydro_gamma * P * density_inv); } /* ------------------------------------------------------------------------- */ #elif defined(EOS_ISOTHERMAL_GAS) ... ... @@ -221,6 +236,22 @@ gas_soundspeed_from_internal_energy(float density, float u) { hydro_gamma_minus_one); } /** * @brief Returns the sound speed given density and pressure * * Since we are using an isothermal EoS, the pressure value is ignored * Computes \f$c = \sqrt{u_{cst} \gamma \rho^{\gamma-1}}\f$. * * @param density The density \f$\rho\f$ * @param P The pressure \f$P\f$ */ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure( float density, float P) { return sqrtf(const_isothermal_internal_energy * hydro_gamma * hydro_gamma_minus_one); } /* ------------------------------------------------------------------------- */ #else ... ...
 ... ... @@ -126,6 +126,17 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy( struct part *restrict p, float u) { p->entropy = gas_entropy_from_internal_energy(p->rho, u); /* Compute the new pressure */ const float pressure = gas_pressure_from_internal_energy(p->rho, u); /* Compute the new sound speed */ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); const float rho_inv = 1.f / p->rho; p->force.soundspeed = soundspeed; p->force.P_over_rho2 = pressure * rho_inv * rho_inv; } /** ... ... @@ -141,6 +152,17 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy( struct part *restrict p, float S) { p->entropy = S; /* Compute the pressure */ const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); /* Compute the new sound speed */ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); const float rho_inv = 1.f / p->rho; p->force.soundspeed = soundspeed; p->force.P_over_rho2 = pressure * rho_inv * rho_inv; } /** ... ... @@ -195,7 +217,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part( p->density.wcount = 0.f; p->density.wcount_dh = 0.f; p->rho = 0.f; p->rho_dh = 0.f; p->density.rho_dh = 0.f; p->density.div_v = 0.f; p->density.rot_v = 0.f; p->density.rot_v = 0.f; ... ... @@ -222,20 +244,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( /* Final operation on the density (add self-contribution). */ p->rho += p->mass * kernel_root; p->rho_dh -= hydro_dimension * p->mass * kernel_root; p->density.rho_dh -= hydro_dimension * p->mass * kernel_root; p->density.wcount += kernel_root; /* Finish the calculation by inserting the missing h-factors */ p->rho *= h_inv_dim; p->rho_dh *= h_inv_dim_plus_one; p->density.rho_dh *= h_inv_dim_plus_one; p->density.wcount *= kernel_norm; p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm; const float rho_inv = 1.f / p->rho; /* Compute the derivative term */ p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * rho_inv); /* Finish calculation of the velocity curl components */ p->density.rot_v *= h_inv_dim_plus_one * rho_inv; p->density.rot_v *= h_inv_dim_plus_one * rho_inv; ... ... @@ -273,19 +292,23 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase; const float pressure = hydro_get_pressure(p, half_dt); const float rho_inv = 1.f / p->rho; /* Divide the pressure by the density and density gradient */ const float P_over_rho2 = pressure * rho_inv * rho_inv * p->rho_dh; /* Compute the sound speed */ const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv); const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); /* Divide the pressure by the density squared to get the SPH term */ const float rho_inv = 1.f / p->rho; const float P_over_rho2 = pressure * rho_inv * rho_inv; /* Compute the Balsara switch */ const float balsara = abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h); /* Compute the "grad h" term */ const float grad_h_term = 1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv); /* Update variables. */ p->force.f = grad_h_term; p->force.P_over_rho2 = P_over_rho2; p->force.soundspeed = soundspeed; p->force.balsara = balsara; ... ... @@ -349,17 +372,16 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase; const float pressure = hydro_get_pressure(p, dt_entr); const float rho_inv = 1.f / p->rho; /* Divide the pressure by the density and density gradient */ const float P_over_rho2 = pressure * rho_inv * rho_inv * p->rho_dh; /* Compute the new sound speed */ const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv); const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); /* Divide the pressure by the density squared to get the SPH term */ const float rho_inv = 1.f / p->rho; const float P_over_rho2 = pressure * rho_inv * rho_inv; /* Update variables */ p->force.P_over_rho2 = P_over_rho2; p->force.soundspeed = soundspeed; p->force.P_over_rho2 = P_over_rho2; } /** ... ... @@ -400,6 +422,19 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( /* Do not 'overcool' when timestep increases */ if (p->entropy + p->entropy_dt * half_dt < 0.5f * p->entropy) p->entropy_dt = -0.5f * p->entropy / half_dt; /* Compute the pressure */ const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); /* Compute the new sound speed */ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); /* Divide the pressure by the density squared to get the SPH term */ const float rho_inv = 1.f / p->rho; const float P_over_rho2 = pressure * rho_inv * rho_inv; p->force.soundspeed = soundspeed; p->force.P_over_rho2 = P_over_rho2; } /** ... ... @@ -414,6 +449,19 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( /* We read u in the entropy field. We now get S from u */ p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy); /* Compute the pressure */ const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); /* Compute the sound speed */ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); /* Divide the pressure by the density squared to get the SPH term */ const float rho_inv = 1.f / p->rho; const float P_over_rho2 = pressure * rho_inv * rho_inv; p->force.soundspeed = soundspeed; p->force.P_over_rho2 = P_over_rho2; } #endif /* SWIFT_GADGET2_HYDRO_H */
 ... ... @@ -30,8 +30,8 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle( "v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n", p->x, p->x, p->x, p->v, p->v, p->v, xp->v_full, xp->v_full, xp->v_full, p->a_hydro, p->a_hydro, p->a_hydro, p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh, p->rho, hydro_get_pressure(p, 0.), p->force.P_over_rho2, p->entropy, p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->density.rho_dh, p->rho, hydro_get_pressure(p, 0.), p->force.P_over_rho2, p->entropy, p->entropy_dt, p->force.soundspeed, p->density.div_v, p->density.rot_v, p->density.rot_v, p->density.rot_v, p->force.balsara, p->force.v_sig, p->force.h_dt, p->ti_begin, p->ti_end); ... ...
 ... ... @@ -59,7 +59,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( /* Compute contribution to the density */ pi->rho += mj * wi; pi->rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx); pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx); /* Compute contribution to the number of neighbours */ pi->density.wcount += wi; ... ... @@ -72,7 +72,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( /* Compute contribution to the density */ pj->rho += mi * wj; pj->rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx); pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx); /* Compute contribution to the number of neighbours */ pj->density.wcount += wj; ... ... @@ -209,13 +209,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density( /* Update particles. */ for (k = 0; k < VEC_SIZE; k++) { pi[k]->rho += rhoi.f[k]; pi[k]->rho_dh -= rhoi_dh.f[k]; pi[k]->density.rho_dh -= rhoi_dh.f[k]; pi[k]->density.wcount += wcounti.f[k]; pi[k]->density.wcount_dh -= wcounti_dh.f[k]; pi[k]->density.div_v -= div_vi.f[k]; for (j = 0; j < 3; j++) pi[k]->density.rot_v[j] += curl_vi[j].f[k]; pj[k]->rho += rhoj.f[k]; pj[k]->rho_dh -= rhoj_dh.f[k]; pj[k]->density.rho_dh -= rhoj_dh.f[k]; pj[k]->density.wcount += wcountj.f[k]; pj[k]->density.wcount_dh -= wcountj_dh.f[k]; pj[k]->density.div_v -= div_vj.f[k]; ... ... @@ -254,7 +254,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( /* Compute contribution to the density */ pi->rho += mj * wi; pi->rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx); pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx); /* Compute contribution to the number of neighbours */ pi->density.wcount += wi; ... ... @@ -366,7 +366,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj, /* Update particles. */ for (k = 0; k < VEC_SIZE; k++) { pi[k]->rho += rhoi.f[k]; pi[k]->rho_dh -= rhoi_dh.f[k]; pi[k]->density.rho_dh -= rhoi_dh.f[k]; pi[k]->density.wcount += wcounti.f[k]; pi[k]->density.wcount_dh -= wcounti_dh.f[k]; pi[k]->density.div_v -= div_vi.f[k]; ... ... @@ -415,7 +415,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( kernel_deval(xj, &wj, &wj_dx); const float wj_dr = hjd_inv * wj_dx; /* Compute gradient terms */ /* Compute h-gradient terms */ const float f_i = pi->force.f; const float f_j = pj->force.f; /* Compute pressure terms */ const float P_over_rho2_i = pi->force.P_over_rho2; const float P_over_rho2_j = pj->force.P_over_rho2; ... ... @@ -447,7 +451,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( /* Now, convolve with the kernel */ const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv; const float sph_term = (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv; (f_i * P_over_rho2_i * wi_dr + f_j * P_over_rho2_j * wj_dr) * r_inv; /* Eventually got the acceleration */ const float acc = visc_term + sph_term; ... ... @@ -690,7 +694,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( kernel_deval(xj, &wj, &wj_dx); const float wj_dr = hjd_inv * wj_dx; /* Compute gradient terms */ /* Compute h-gradient terms */ const float f_i = pi->force.f; const float f_j = pj->force.f; /* Compute pressure terms */ const float P_over_rho2_i = pi->force.P_over_rho2; const float P_over_rho2_j = pj->force.P_over_rho2; ... ... @@ -722,7 +730,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( /* Now, convolve with the kernel */ const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv; const float sph_term = (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv; (f_i * P_over_rho2_i * wi_dr + f_j * P_over_rho2_j * wj_dr) * r_inv; /* Eventually got the acceleration */ const float acc = visc_term + sph_term; ... ...
 ... ... @@ -80,9 +80,6 @@ struct part { /* Entropy time derivative */ float entropy_dt; /* Derivative of the density with respect to smoothing length. */ float rho_dh; union { struct { ... ... @@ -93,6 +90,9 @@ struct part { /* Number of neighbours spatial derivative. */ float wcount_dh; /* Derivative of the density with respect to h. */ float rho_dh; /* Particle velocity curl. */ float rot_v; ... ... @@ -106,6 +106,9 @@ struct part { /* Balsara switch */ float balsara; /*! "Grad h" term */ float f; /* Pressure over density squared (including drho/dh term) */ float P_over_rho2; ... ...
 ... ... @@ -65,9 +65,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy( __attribute__((always_inline)) INLINE static float hydro_get_pressure( const struct part *restrict p, float dt) { const float u = p->u + p->u_dt * dt; return gas_pressure_from_internal_energy(p->rho, u); return p->force.pressure; } /** ... ... @@ -97,9 +95,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy( __attribute__((always_inline)) INLINE static float hydro_get_soundspeed( const struct part *restrict p, float dt) { const float u = p->u + p->u_dt * dt; return gas_soundspeed_from_internal_energy(p->rho, u); return p->force.soundspeed; } /** ... ... @@ -138,6 +134,15 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy( struct part *restrict p, float u) { p->u = u; /* Compute the new pressure */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); /* Compute the new sound speed */ const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u); p->force.soundspeed = soundspeed; p->force.pressure = pressure; } /** ... ... @@ -153,6 +158,15 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy( struct part *restrict p, float S) { p->u = gas_internal_energy_from_entropy(p->rho, S); /* Compute the pressure */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); /* Compute the new sound speed */ const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u); p->force.soundspeed = soundspeed; p->force.pressure = pressure; } /** ... ... @@ -214,7 +228,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part( p->density.wcount = 0.f; p->density.wcount_dh = 0.f; p->rho = 0.f; p->rho_dh = 0.f; p->density.rho_dh = 0.f; } /** ... ... @@ -240,19 +254,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( /* Final operation on the density (add self-contribution). */ p->rho += p->mass * kernel_root; p->rho_dh -= hydro_dimension * p->mass * kernel_root; p->density.rho_dh -= hydro_dimension * p->mass * kernel_root; p->density.wcount += kernel_root; /* Finish the calculation by inserting the missing h-factors */ p->rho *= h_inv_dim; p->rho_dh *= h_inv_dim_plus_one; p->density.rho_dh *= h_inv_dim_plus_one; p->density.wcount *= kernel_norm; p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm; const float irho = 1.f / p->rho; /* Compute the derivative term */ p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho); } /** ... ... @@ -274,10 +283,22 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( struct part *restrict p, struct xpart *restrict xp, int ti_current, double timeBase) { /* Compute the pressure */ const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase; const float pressure = hydro_get_pressure(p, half_dt); /* Compute the sound speed */ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); /* Compute the "grad h" term */ const float rho_inv = 1.f / p->rho; const float grad_h_term = 1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv); /* Update variables. */ p->force.f = grad_h_term; p->force.pressure = pressure; p->force.soundspeed = soundspeed; } /** ... ... @@ -337,7 +358,13 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( /* Drift the pressure */ const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase; p->force.pressure = hydro_get_pressure(p, dt_entr); const float pressure = hydro_get_pressure(p, dt_entr); /* Compute the new sound speed */ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); p->force.pressure = pressure; p->force.soundspeed = soundspeed; } /** ... ... @@ -379,6 +406,15 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( /* Do not 'overcool' when timestep increases */ if (p->u + p->u_dt * half_dt < 0.5f * p->u) p->u_dt = -0.5f * p->u / half_dt; /* Compute the pressure */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); /* Compute the sound speed */ const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u); p->force.pressure = pressure; p->force.soundspeed = soundspeed; } /** ... ... @@ -392,6 +428,16 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( * @param p The particle to act upon */ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( struct part *restrict p) {} struct part *restrict p) { /* Compute the pressure */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); /* Compute the sound speed */ const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u); p->force.pressure = pressure; p->force.soundspeed = soundspeed; } #endif /* SWIFT_MINIMAL_HYDRO_H */
 ... ... @@ -44,7 +44,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle( p->x, p->x, p->x, p->v, p->v, p->v, xp->v_full, xp->v_full, xp->v_full, p->a_hydro, p->a_hydro, p->a_hydro, p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h, p->force.h_dt, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho, p->ti_begin, p->ti_end); } ... ...
 ... ... @@ -55,7 +55,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( kernel_deval(xi, &wi, &wi_dx); pi->rho += mj * wi; pi->rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx); pi->density.rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx); pi->density.wcount += wi; pi->density.wcount_dh -= xi * wi_dx; ... ... @@ -65,7 +65,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( kernel_deval(xj, &wj, &wj_dx); pj->rho += mi * wj; pj->rho_dh -= mi * (hydro_dimension * wj + xj * wj_dx); pj->density.rho_dh -= mi * (hydro_dimension * wj + xj * wj_dx); pj->density.wcount += wj; pj->density.wcount_dh -= xj * wj_dx; } ... ... @@ -100,7 +100,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( kernel_deval(xi, &wi, &wi_dx); pi->rho += mj * wi; pi->rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx); pi->density.rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx); pi->density.wcount += wi; pi->density.wcount_dh -= xi * wi_dx; } ... ... @@ -152,8 +152,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( const float wj_dr = hjd_inv * wj_dx; /* Compute gradient terms */ const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh; const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh; const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f; const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f; /* Compute dv dot r. */ const float dvdr = (pi->v - pj->v) * dx + ... ... @@ -165,8 +165,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ /* Compute sound speeds and signal velocity */ const float ci = sqrtf(hydro_gamma * pressurei / rhoi); const float cj = sqrtf(hydro_gamma * pressurej / rhoj); const float ci = pi->force.soundspeed; const float cj = pj->force.soundspeed; const float v_sig = ci + cj - 3.f * mu_ij; /* Construct the full viscosity term */ ... ... @@ -263,8 +263,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( const float wj_dr = hjd_inv * wj_dx; /* Compute gradient terms */ const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh; const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh; const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f; const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f; /* Compute dv dot r. */ const float dvdr = (pi->v - pj->v) * dx + ... ... @@ -276,8 +276,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ /* Compute sound speeds and signal velocity */ const float ci = sqrtf(hydro_gamma * pressurei / rhoi); const float cj = sqrtf(hydro_gamma * pressurej / rhoj); const float ci = pi->force.soundspeed; const float cj = pj->force.soundspeed; const float v_sig = ci + cj - 3.f * mu_ij; /* Construct the full viscosity term */ ... ...
 ... ... @@ -93,9 +93,6 @@ struct part { /*! Particle density. */ float rho; /*! Derivative of density with respect to h */ float rho_dh; /* Store density/force specific stuff. */ union { ... ... @@ -114,6 +111,9 @@ struct part { /*! Derivative of the neighbour number with respect to h. */ float wcount_dh; /*! Derivative of density with respect to h */ float rho_dh; } density; /** ... ... @@ -125,9 +125,15 @@ struct part { */ struct { /*! "Grad h" term */ float f; /*! Particle pressure. */ float pressure; /*! Particle soundspeed. */ float soundspeed; /*! Particle signal velocity */ float v_sig;