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

Leaner version of Gadget2 particles.

parent 1a1dd00b
No related branches found
No related tags found
2 merge requests!136Master,!90Improved multi-timestep SPH
Hide whitespace changes
Inline Side-by-side
src/hydro/Gadget2/hydro.h
+ 32
− 30
View file @ 7a607502
......@@ -36,7 +36,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const float dt_accel = sqrtf(2.f); // MATTHIEU
/* CFL condition */
const float dt_cfl = 2.f * const_cfl * kernel_gamma * p->h / p->v_sig;
const float dt_cfl = 2.f * const_cfl * kernel_gamma * p->h / p->force.v_sig;
return fminf(dt_cfl, dt_accel);
}
......@@ -56,10 +56,9 @@ __attribute__((always_inline))
p->rho = 0.f;
p->rho_dh = 0.f;
p->div_v = 0.f;
p->curl_v = 0.f;
p->rot_v[0] = 0.f;
p->rot_v[1] = 0.f;
p->rot_v[2] = 0.f;
p->density.rot_v[0] = 0.f;
p->density.rot_v[1] = 0.f;
p->density.rot_v[2] = 0.f;
}
/**
......@@ -94,23 +93,13 @@ __attribute__((always_inline))
/* Compute the derivative term */
p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh / p->rho);
/* Finish calculation of the velocity curl */
p->rot_v[0] *= ih4;
p->rot_v[1] *= ih4;
p->rot_v[2] *= ih4;
/* And compute the norm of the curl */
p->curl_v = sqrtf(p->rot_v[0] * p->rot_v[0] + p->rot_v[1] * p->rot_v[1] +
p->rot_v[2] * p->rot_v[2]) /
p->rho;
/* Finish calculation of the velocity curl components */
p->density.rot_v[0] *= ih4 / p->rho;
p->density.rot_v[1] *= ih4 / p->rho;
p->density.rot_v[2] *= ih4 / p->rho;
/* Finish calculation of the velocity divergence */
p->div_v *= ih4 / p->rho;
/* Compute the pressure */
const float dt = time - 0.5f * (p->t_begin + p->t_end);
p->pressure =
(p->entropy + p->entropy_dt * dt) * powf(p->rho, const_hydro_gamma);
}
/**
......@@ -120,10 +109,23 @@ __attribute__((always_inline))
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
* @param time The current time
*/
__attribute__((always_inline))
INLINE static void hydro_prepare_force(struct part* p, struct xpart* xp) {
INLINE static void hydro_prepare_force(struct part* p, struct xpart* xp, float time) {
/* Compute the norm of the curl */
p->force.curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
p->density.rot_v[1] * p->density.rot_v[1] +
p->density.rot_v[2] * p->density.rot_v[2]);
/* Compute the pressure */
const float dt = time - 0.5f * (p->t_begin + p->t_end);
p->force.pressure =
(p->entropy + p->force.entropy_dt * dt) * powf(p->rho, const_hydro_gamma);
/* Some smoothing length multiples. */
/* const float h = p->h; */
/* const float ih = 1.0f / h; */
......@@ -173,10 +175,10 @@ __attribute__((always_inline))
p->force.h_dt = 0.0f;
/* Reset the time derivatives. */
p->entropy_dt = 0.0f;
p->force.entropy_dt = 0.0f;
/* Reset maximal signal velocity */
p->v_sig = 0.0f;
p->force.v_sig = 0.0f;
}
/**
......@@ -194,8 +196,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
// p->h *= expf(0.33333333f * p->div_v * dt)
const float dt_entr = t1 - 0.5f * (p->t_begin + p->t_end);
p->pressure =
(p->entropy + p->entropy_dt * dt_entr) * powf(p->rho, const_hydro_gamma);
p->force.pressure =
(p->entropy + p->force.entropy_dt * dt_entr) * powf(p->rho, const_hydro_gamma);
}
/**
......@@ -208,7 +210,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
__attribute__((always_inline))
INLINE static void hydro_end_force(struct part* p) {
p->entropy_dt *=
p->force.entropy_dt *=
(const_hydro_gamma - 1.f) * powf(p->rho, -(const_hydro_gamma - 1.f));
}
......@@ -221,15 +223,15 @@ __attribute__((always_inline))
INLINE static void hydro_kick_extra(struct part* p, float dt) {
/* Do not decrease the entropy (temperature) by more than a factor of 2*/
const float entropy_change = p->entropy_dt * dt;
const float entropy_change = p->force.entropy_dt * dt;
if (entropy_change > -0.5f * p->entropy)
p->entropy += entropy_change;
else
p->entropy *= 0.5f;
p->entropy *= 0.5f;
/* Do not 'overcool' when timestep increases */
if (p->entropy + 0.5f * p->entropy_dt * dt < 0.5f * p->entropy)
p->entropy_dt = -0.5f * p->entropy / dt;
if (p->entropy + 0.5f * p->force.entropy_dt * dt < 0.5f * p->entropy)
p->force.entropy_dt = -0.5f * p->entropy / dt;
}
/**
......@@ -243,5 +245,5 @@ __attribute__((always_inline))
INLINE static void hydro_convert_quantities(struct part* p) {
p->entropy = (const_hydro_gamma - 1.f) * p->entropy *
powf(p->rho, -(const_hydro_gamma - 1.f));
powf(p->rho, -(const_hydro_gamma - 1.f));
}
src/hydro/Gadget2/hydro_debug.h
+ 3
− 3
View file @ 7a607502
......@@ -29,7 +29,7 @@ __attribute__((always_inline))
"v_sig=%e t_begin=%.3e, t_end=%.3e\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->pressure,
p->entropy, p->entropy_dt, p->div_v, p->curl_v, p->rot_v[0], p->rot_v[1],
p->rot_v[2], p->v_sig, p->t_begin, p->t_end);
p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->force.pressure,
p->entropy, p->force.entropy_dt, p->div_v, p->force.curl_v, p->density.rot_v[0],
p->density.rot_v[1], p->density.rot_v[2], p->force.v_sig, p->t_begin, p->t_end);
}
src/hydro/Gadget2/hydro_iact.h
+ 23
− 23
View file @ 7a607502
......@@ -103,13 +103,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
pi->rot_v[0] += faci * curlvr[0];
pi->rot_v[1] += faci * curlvr[1];
pi->rot_v[2] += faci * curlvr[2];
pi->density.rot_v[0] += faci * curlvr[0];
pi->density.rot_v[1] += faci * curlvr[1];
pi->density.rot_v[2] += faci * curlvr[2];
pj->rot_v[0] += facj * curlvr[0];
pj->rot_v[1] += facj * curlvr[1];
pj->rot_v[2] += facj * curlvr[2];
pj->density.rot_v[0] += facj * curlvr[0];
pj->density.rot_v[1] += facj * curlvr[1];
pj->density.rot_v[2] += facj * curlvr[2];
}
/**
......@@ -177,9 +177,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
pi->rot_v[0] += fac * curlvr[0];
pi->rot_v[1] += fac * curlvr[1];
pi->rot_v[2] += fac * curlvr[2];
pi->density.rot_v[0] += fac * curlvr[0];
pi->density.rot_v[1] += fac * curlvr[1];
pi->density.rot_v[2] += fac * curlvr[2];
}
/**
......@@ -201,8 +201,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float mj = pj->mass;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
const float pressurei = pi->pressure;
const float pressurej = pj->pressure;
const float pressurei = pi->force.pressure;
const float pressurej = pj->force.pressure;
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
......@@ -238,9 +238,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float mu_ij = fac_mu * dvdr * r_inv;
v_sig -= 3.f * mu_ij;
const float rho_ij = 0.5f * (rhoi + rhoj);
const float balsara_i = fabsf(pi->div_v) / (fabsf(pi->div_v) + pi->curl_v +
const float balsara_i = fabsf(pi->div_v) / (fabsf(pi->div_v) + pi->force.curl_v +
0.0001 * ci / fac_mu / hi);
const float balsara_j = fabsf(pj->div_v) / (fabsf(pj->div_v) + pj->curl_v +
const float balsara_j = fabsf(pj->div_v) / (fabsf(pj->div_v) + pj->force.curl_v +
0.0001 * cj / fac_mu / hj);
visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij / rho_ij *
(balsara_i + balsara_j);
......@@ -285,12 +285,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->force.h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
/* Update the signal velocity. */
pi->v_sig = fmaxf(pi->v_sig, v_sig);
pj->v_sig = fmaxf(pj->v_sig, v_sig);
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
pj->force.v_sig = fmaxf(pj->force.v_sig, v_sig);
/* Change in entropy */
pi->entropy_dt += 0.5f * visc_term * dvdr;
pj->entropy_dt -= 0.5f * visc_term * dvdr;
pi->force.entropy_dt += 0.5f * visc_term * dvdr;
pj->force.entropy_dt -= 0.5f * visc_term * dvdr;
}
/**
......@@ -312,8 +312,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float mj = pj->mass;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
const float pressurei = pi->pressure;
const float pressurej = pj->pressure;
const float pressurei = pi->force.pressure;
const float pressurej = pj->force.pressure;
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
......@@ -349,9 +349,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float mu_ij = fac_mu * dvdr * r_inv;
v_sig -= 3.f * mu_ij;
const float rho_ij = 0.5f * (rhoi + rhoj);
const float balsara_i = fabsf(pi->div_v) / (fabsf(pi->div_v) + pi->curl_v +
const float balsara_i = fabsf(pi->div_v) / (fabsf(pi->div_v) + pi->force.curl_v +
0.0001 * ci / fac_mu / hi);
const float balsara_j = fabsf(pj->div_v) / (fabsf(pj->div_v) + pj->curl_v +
const float balsara_j = fabsf(pj->div_v) / (fabsf(pj->div_v) + pj->force.curl_v +
0.0001 * cj / fac_mu / hj);
visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij / rho_ij *
(balsara_i + balsara_j);
......@@ -374,10 +374,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
/* Update the signal velocity. */
pi->v_sig = fmaxf(pi->v_sig, v_sig);
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
/* Change in entropy */
pi->entropy_dt += 0.5f * visc_term * dvdr;
pi->force.entropy_dt += 0.5f * visc_term * dvdr;
}
#endif /* SWIFT_RUNNER_IACT_LEGACY_H */
src/hydro/Gadget2/hydro_part.h
+ 42
− 38
View file @ 7a607502
......@@ -17,10 +17,8 @@
*
******************************************************************************/
/* Some standard headers. */
#include <stdlib.h>
/* Extra particle data not needed during the computation. */
/* Extra particle data not needed during the SPH loops over neighbours. */
struct xpart {
/* Old position, at last tree rebuild. */
......@@ -28,7 +26,7 @@ struct xpart {
/* Velocity at the last full step. */
float v_full[3];
} __attribute__((aligned(xpart_align)));
/* Data of a single particle. */
......@@ -52,26 +50,6 @@ struct part {
/* Particle time of end of time-step. */
float t_end;
struct {
/* Number of neighbours */
float wcount;
/* Number of neighbours spatial derivative */
float wcount_dh;
} density;
struct {
/* Time derivative of the smoothing length */
float h_dt;
} force;
/* Particle entropy. */
float entropy;
/* Particle density. */
float rho;
......@@ -79,25 +57,51 @@ struct part {
*/
float rho_dh;
/* Particle entropy. */
float entropy;
/* Particle mass. */
float mass;
/* Particle pressure */
float pressure;
/* Entropy time derivative */
float entropy_dt;
union{
struct {
/* Number of neighbours */
float wcount;
/* Number of neighbours spatial derivative */
float wcount_dh;
/* Velocity curl components */
float rot_v[3];
} density;
struct {
/* Time derivative of the smoothing length */
float h_dt;
/* Velocity curl norm*/
float curl_v;
/* Signal velocity */
float v_sig;
/* Entropy time derivative */
float entropy_dt;
/* Particle pressure */
float pressure;
} force;
};
/* Velocity divergence */
float div_v;
/* Velocity curl */
float curl_v;
float rot_v[3];
/* Signal velocity */
float v_sig;
/* Particle ID. */
unsigned long long id;
......
src/runner.c
+ 6
− 6
View file @ 7a607502
......@@ -620,15 +620,15 @@ void runner_doghost(struct runner *r, struct cell *c) {
}
/* We now have a particle whose smoothing length has converged */
// if(p->id == 1000)
// printParticle(parts, 1000, count);
/* As of here, particle force variables will be set. Do _NOT_
try to read any particle density variables! */
/* As of here, particle force variables will be set. */
/* Compute variables required for the force loop */
hydro_prepare_force(p, xp);
hydro_prepare_force(p, xp, t_end);
/* The particle force values are now set. Do _NOT_
try to read any particle density variables! */
/* Prepare the particle for the force loop over neighbours */
hydro_reset_acceleration(p);
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment