Commit d5e142f0 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

h_dt not part of the force sub-structure anymore

parent aaf537a9
...@@ -1731,13 +1731,15 @@ void engine_init_particles(struct engine *e) { ...@@ -1731,13 +1731,15 @@ void engine_init_particles(struct engine *e) {
TIMER_TOC(timer_runners); TIMER_TOC(timer_runners);
// message("\n0th ENTROPY CONVERSION\n"); /* message("\n0th ENTROPY CONVERSION\n"); */
space_map_cells_pre(s, 1, cell_convert_hydro, NULL); space_map_cells_pre(s, 1, cell_convert_hydro, NULL);
// printParticle(e->s->parts, e->s->xparts,1000, e->s->nr_parts); /* printParticle(e->s->parts, e->s->xparts,1000, e->s->nr_parts); */
// printParticle(e->s->parts, e->s->xparts,515050, e->s->nr_parts); /* printParticle(e->s->parts, e->s->xparts,515050, e->s->nr_parts); */
// exit(0);
/* Ready to go */ /* Ready to go */
e->step = -1; e->step = -1;
} }
...@@ -1808,7 +1810,7 @@ if ( e->nodeID == 0 ) ...@@ -1808,7 +1810,7 @@ if ( e->nodeID == 0 )
message( "nr_parts=%i." , nr_parts ); */ message( "nr_parts=%i." , nr_parts ); */
#endif #endif
// message("\nDRIFT\n"); /* message("\nDRIFT\n"); */
/* Move forward in time */ /* Move forward in time */
e->timeOld = e->time; e->timeOld = e->time;
...@@ -1819,12 +1821,12 @@ if ( e->nodeID == 0 ) ...@@ -1819,12 +1821,12 @@ if ( e->nodeID == 0 )
/* Drift everybody */ /* Drift everybody */
engine_launch(e, e->nr_threads, 1 << task_type_drift, 0); engine_launch(e, e->nr_threads, 1 << task_type_drift, 0);
// printParticle(e->s->parts, e->s->xparts, 1000, e->s->nr_parts); /* printParticle(e->s->parts, e->s->xparts, 1000, e->s->nr_parts); */
// printParticle(e->s->parts, e->s->xparts, 515050, e->s->nr_parts); /* printParticle(e->s->parts, e->s->xparts, 515050, e->s->nr_parts); */
// if(e->step == 2) exit(0); // if(e->step == 2) exit(0);
// message("\nACCELERATION AND KICK\n"); /* message("\nACCELERATION AND KICK\n"); */
/* Re-distribute the particles amongst the nodes? */ /* Re-distribute the particles amongst the nodes? */
if (e->forcerepart) engine_repartition(e); if (e->forcerepart) engine_repartition(e);
...@@ -1850,8 +1852,11 @@ if ( e->nodeID == 0 ) ...@@ -1850,8 +1852,11 @@ if ( e->nodeID == 0 )
((double)timers[timer_count - 1]) / CPU_TPS * 1000); ((double)timers[timer_count - 1]) / CPU_TPS * 1000);
fflush(stdout); fflush(stdout);
// printParticle(e->s->parts, e->s->xparts,1000, e->s->nr_parts); /* printParticle(e->s->parts, e->s->xparts,1000, e->s->nr_parts); */
// printParticle(e->s->parts, e->s->xparts,515050, e->s->nr_parts); /* printParticle(e->s->parts, e->s->xparts,515050, e->s->nr_parts); */
/* if(e->step == 2) */
/* exit(0); */
} }
/** /**
......
...@@ -31,8 +31,8 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep( ...@@ -31,8 +31,8 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
float dt_cfl = const_cfl * p->h / p->force.v_sig; float dt_cfl = const_cfl * p->h / p->force.v_sig;
/* Limit change in h */ /* Limit change in h */
float dt_h_change = (p->force.h_dt != 0.0f) float dt_h_change = (p->h_dt != 0.0f)
? fabsf(const_ln_max_h_change * p->h / p->force.h_dt) ? fabsf(const_ln_max_h_change * p->h / p->h_dt)
: FLT_MAX; : FLT_MAX;
/* Limit change in u */ /* Limit change in u */
...@@ -159,7 +159,7 @@ __attribute__((always_inline)) ...@@ -159,7 +159,7 @@ __attribute__((always_inline))
/* Reset the time derivatives. */ /* Reset the time derivatives. */
p->force.u_dt = 0.0f; p->force.u_dt = 0.0f;
p->force.h_dt = 0.0f; p->h_dt = 0.0f;
p->force.v_sig = 0.0f; p->force.v_sig = 0.0f;
} }
...@@ -197,10 +197,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( ...@@ -197,10 +197,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
* @param p The particle to act upon * @param p The particle to act upon
*/ */
__attribute__((always_inline)) __attribute__((always_inline))
INLINE static void hydro_end_force(struct part* p) { INLINE static void hydro_end_force(struct part* p) {}
p->force.h_dt *= p->h * 0.333333333f;
}
/** /**
* @brief Kick the additional variables * @brief Kick the additional variables
......
...@@ -447,8 +447,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -447,8 +447,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->force.u_dt += mi * tc * (pj->u - pi->u); pj->force.u_dt += mi * tc * (pj->u - pi->u);
/* Get the time derivative for h. */ /* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr / rhoj * wi_dr; pi->h_dt -= mj * dvdr / rhoj * wi_dr;
pj->force.h_dt -= mi * dvdr / rhoi * wj_dr; pj->h_dt -= mi * dvdr / rhoi * wj_dr;
/* Update the signal velocity. */ /* Update the signal velocity. */
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig); pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
...@@ -656,8 +656,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force( ...@@ -656,8 +656,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
for (k = 0; k < VEC_SIZE; k++) { for (k = 0; k < VEC_SIZE; k++) {
pi[k]->force.u_dt += piu_dt.f[k]; pi[k]->force.u_dt += piu_dt.f[k];
pj[k]->force.u_dt += pju_dt.f[k]; pj[k]->force.u_dt += pju_dt.f[k];
pi[k]->force.h_dt -= pih_dt.f[k]; pi[k]->h_dt -= pih_dt.f[k];
pj[k]->force.h_dt -= pjh_dt.f[k]; pj[k]->h_dt -= pjh_dt.f[k];
pi[k]->force.v_sig = vi_sig.f[k]; pi[k]->force.v_sig = vi_sig.f[k];
pj[k]->force.v_sig = vj_sig.f[k]; pj[k]->force.v_sig = vj_sig.f[k];
for (j = 0; j < 3; j++) { for (j = 0; j < 3; j++) {
...@@ -758,7 +758,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -758,7 +758,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.u_dt += mj * tc * (pi->u - pj->u); pi->force.u_dt += mj * tc * (pi->u - pj->u);
/* Get the time derivative for h. */ /* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr / rhoj * wi_dr; pi->h_dt -= mj * dvdr / rhoj * wi_dr;
/* Update the signal velocity. */ /* Update the signal velocity. */
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig); pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
...@@ -957,7 +957,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force( ...@@ -957,7 +957,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
/* Store the forces back on the particles. */ /* Store the forces back on the particles. */
for (k = 0; k < VEC_SIZE; k++) { for (k = 0; k < VEC_SIZE; k++) {
pi[k]->force.u_dt += piu_dt.f[k]; pi[k]->force.u_dt += piu_dt.f[k];
pi[k]->force.h_dt -= pih_dt.f[k]; pi[k]->h_dt -= pih_dt.f[k];
pi[k]->force.v_sig = vi_sig.f[k]; pi[k]->force.v_sig = vi_sig.f[k];
pj[k]->force.v_sig = vj_sig.f[k]; pj[k]->force.v_sig = vj_sig.f[k];
for (j = 0; j < 3; j++) pi[k]->a[j] -= pia[j].f[k]; for (j = 0; j < 3; j++) pi[k]->a[j] -= pia[j].f[k];
......
...@@ -46,6 +46,9 @@ struct part { ...@@ -46,6 +46,9 @@ struct part {
/* Particle cutoff radius. */ /* Particle cutoff radius. */
float h; float h;
/* Change in smoothing length over time. */
float h_dt;
/* Particle time of beginning of time-step. */ /* Particle time of beginning of time-step. */
float t_begin; float t_begin;
...@@ -95,9 +98,6 @@ struct part { ...@@ -95,9 +98,6 @@ struct part {
/* Change in particle energy over time. */ /* Change in particle energy over time. */
float u_dt; float u_dt;
/* Change in smoothing length over time. */
float h_dt;
/* Signal velocity */ /* Signal velocity */
float v_sig; float v_sig;
......
...@@ -147,7 +147,7 @@ __attribute__((always_inline)) ...@@ -147,7 +147,7 @@ __attribute__((always_inline))
p->a_hydro[1] = 0.0f; p->a_hydro[1] = 0.0f;
p->a_hydro[2] = 0.0f; p->a_hydro[2] = 0.0f;
p->force.h_dt = 0.0f; p->h_dt = 0.0f;
/* Reset the time derivatives. */ /* Reset the time derivatives. */
p->entropy_dt = 0.0f; p->entropy_dt = 0.0f;
......
...@@ -32,6 +32,6 @@ __attribute__((always_inline)) ...@@ -32,6 +32,6 @@ __attribute__((always_inline))
p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->force.pressure, p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->force.pressure,
p->entropy, p->entropy_dt, p->force.soundspeed, p->entropy, p->entropy_dt, p->force.soundspeed,
p->div_v, p->force.curl_v, p->density.rot_v[0], 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->force.h_dt, p->density.rot_v[1], p->density.rot_v[2], p->force.v_sig, p->h_dt,
p->t_begin, p->t_end); p->t_begin, p->t_end);
} }
...@@ -281,8 +281,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -281,8 +281,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->a_hydro[2] += acc * dx[2]; pj->a_hydro[2] += acc * dx[2];
/* Get the time derivative for h. */ /* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr; pi->h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
pj->force.h_dt -= mi * dvdr * r_inv / rhoi * wj_dr; pj->h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
/* Update the signal velocity. */ /* Update the signal velocity. */
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig); pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
...@@ -371,7 +371,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -371,7 +371,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->a_hydro[2] -= acc * dx[2]; pi->a_hydro[2] -= acc * dx[2];
/* Get the time derivative for h. */ /* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr; pi->h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
/* Update the signal velocity. */ /* Update the signal velocity. */
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig); pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
......
...@@ -44,6 +44,9 @@ struct part { ...@@ -44,6 +44,9 @@ struct part {
/* Particle cutoff radius. */ /* Particle cutoff radius. */
float h; float h;
/* Time derivative of the smoothing length */
float h_dt;
/* Particle time of beginning of time-step. */ /* Particle time of beginning of time-step. */
float t_begin; float t_begin;
...@@ -82,10 +85,7 @@ struct part { ...@@ -82,10 +85,7 @@ struct part {
} density; } density;
struct { struct {
/* Time derivative of the smoothing length */
float h_dt;
/* Velocity curl norm*/ /* Velocity curl norm*/
float curl_v; float curl_v;
......
...@@ -734,14 +734,14 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) { ...@@ -734,14 +734,14 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
p->v[2] += p->a_hydro[2] * dt; p->v[2] += p->a_hydro[2] * dt;
/* Predict smoothing length */ /* Predict smoothing length */
w = p->force.h_dt * h_inv * dt; w = p->h_dt * h_inv * dt;
if (fabsf(w) < 0.2f) if (fabsf(w) < 0.2f)
p->h *= approx_expf(w); /* 4th order expansion of exp(w) */ p->h *= approx_expf(w); /* 4th order expansion of exp(w) */
else else
p->h *= expf(w); p->h *= expf(w);
/* Predict density */ /* Predict density */
w = -3.0f * p->force.h_dt * h_inv * dt; w = -3.0f * p->h_dt * h_inv * dt;
if (fabsf(w) < 0.2f) if (fabsf(w) < 0.2f)
p->rho *= approx_expf(w); /* 4th order expansion of exp(w) */ p->rho *= approx_expf(w); /* 4th order expansion of exp(w) */
else else
...@@ -759,7 +759,7 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) { ...@@ -759,7 +759,7 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
/* p->v[0], */ /* p->v[0], */
/* p->v[1], */ /* p->v[1], */
/* p->v[2], */ /* p->v[2], */
/* p->force.h_dt * h_inv * dt, */ /* p->h_dt * h_inv * dt, */
/* 0.333333f * p->div_v * dt); */ /* 0.333333f * p->div_v * dt); */
/* Compute motion since last cell construction */ /* Compute motion since last cell construction */
...@@ -852,7 +852,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) { ...@@ -852,7 +852,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
if (is_fixdt || p->t_end <= t_current) { if (is_fixdt || p->t_end <= t_current) {
/* First, finish the force loop */ /* First, finish the force loop */
p->force.h_dt *= p->h * 0.333333333f; p->h_dt *= p->h * 0.333333333f;
/* And do the same of the extra variable */ /* And do the same of the extra variable */
hydro_end_force(p); hydro_end_force(p);
......
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