Commit 8633109d authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Fixed the minimal SPH hydro version

parent 22f10e35
......@@ -100,7 +100,7 @@ void printgParticle(struct gpart *gparts, long long int id, size_t N) {
void printParticle_single(struct part *p, struct xpart *xp) {
printf("## Particle: id=%lld", p->id);
printf("## Particle: id=%lld ", p->id);
hydro_debug_particle(p, xp);
}
......
......@@ -1904,18 +1904,18 @@ void engine_init_particles(struct engine *e) {
if (e->nodeID == 0) message("Initialising particles");
engine_prepare(e);
/* Make sure all particles are ready to go */
/* i.e. clean-up any stupid state in the ICs */
if (e->policy & engine_policy_hydro) {
space_map_cells_pre(s, 1, cell_init_parts, NULL);
space_map_cells_pre(s, 0, cell_init_parts, NULL);
}
if ((e->policy & engine_policy_self_gravity) ||
(e->policy & engine_policy_external_gravity)) {
space_map_cells_pre(s, 1, cell_init_gparts, NULL);
space_map_cells_pre(s, 0, cell_init_gparts, NULL);
}
engine_prepare(e);
engine_marktasks(e);
/* Build the masks corresponding to the policy */
......
......@@ -172,14 +172,7 @@ __attribute__((always_inline))
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {
const float dt = t1 - t0;
/* Predict internal energy */
const float w = p->u_dt / p->u * dt;
if (fabsf(w) < 0.2f)
p->u *= approx_expf(w); /* 4th order expansion of exp(w) */
else
p->u *= expf(w);
p->u = xp->u_full;
/* Need to recompute the pressure as well */
p->force.pressure = p->rho * p->u * (const_hydro_gamma - 1.f);
......
......@@ -131,18 +131,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute dv dot r. */
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];
dvdr *= r_inv;
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];
/* Compute the relative velocity. (This is 0 if the particles move away from
* each other and negative otherwise) */
const float omega_ij = fminf(dvdr, 0.f);
/* Compute sound speeds */
const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
const float v_sig = ci + cj + 3.f * omega_ij;
const float v_sig = ci + cj;
/* SPH acceleration term */
const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
......@@ -157,8 +154,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->a_hydro[2] += mi * sph_term * dx[2];
/* Get the time derivative for u. */
pi->u_dt += P_over_rho_i * mj * dvdr * wi_dr;
pj->u_dt += P_over_rho_j * mi * dvdr * wj_dr;
pi->u_dt += P_over_rho_i * mj * dvdr * r_inv * wi_dr;
pj->u_dt += P_over_rho_j * mi * dvdr * r_inv * wj_dr;
/* Get the time derivative for h. */
pi->h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
......@@ -208,18 +205,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute dv dot r. */
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];
dvdr *= r_inv;
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];
/* Compute the relative velocity. (This is 0 if the particles move away from
* each other and negative otherwise) */
const float omega_ij = fminf(dvdr, 0.f);
/* Compute sound speeds */
const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
const float v_sig = ci + cj + 3.f * omega_ij;
const float v_sig = ci + cj;
/* SPH acceleration term */
const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
......@@ -230,7 +224,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->a_hydro[2] -= mj * sph_term * dx[2];
/* Get the time derivative for u. */
pi->u_dt += P_over_rho_i * mj * dvdr * wi_dr;
pi->u_dt += P_over_rho_i * mj * dvdr * r_inv * wi_dr;
/* Get the time derivative for h. */
pi->h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment