diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h index 466bc72757be625ab71c90003269823845733291..fca4a346047d7dce0741924a69e95fdad5a5ce45 100644 --- a/src/hydro/Default/hydro.h +++ b/src/hydro/Default/hydro.h @@ -17,6 +17,8 @@ * ******************************************************************************/ +#include "approx_math.h" + /** * @brief Computes the hydro time-step of a given particle * @@ -108,7 +110,7 @@ __attribute__((always_inline)) * @param time The current time */ __attribute__((always_inline)) INLINE static void hydro_prepare_force( - struct part* p, struct xpart* xp, float time) { + struct part* p, struct xpart* xp, int ti_current, double timeBase) { /* Some smoothing length multiples. */ const float h = p->h; @@ -146,7 +148,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( (const_viscosity_alpha_max - p->alpha) * S; /* Update particle's viscosity paramter */ - p->alpha += alpha_dot * (p->t_end - p->t_begin); + p->alpha += alpha_dot * (p->ti_end - p->ti_begin) * timeBase; } /** @@ -178,18 +180,18 @@ __attribute__((always_inline)) * @param xp The extended data of the particle * @param t0 The time at the start of the drift * @param t1 The time at the end of the drift + * @param timeBase The minimal time-step size */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( - struct part* p, struct xpart* xp, float t0, float t1) { + struct part* p, struct xpart* xp, int t0, int t1, double timeBase) { float u, w; - const float dt = t1 - t0; + const float dt = (t1 - t0) * timeBase; /* Predict internal energy */ w = p->force.u_dt / p->u * dt; - if (fabsf(w) < 0.01f) /* 1st order expansion of exp(w) */ - u = p->u *= - 1.0f + w * (1.0f + w * (0.5f + w * (1.0f / 6.0f + 1.0f / 24.0f * w))); + if (fabsf(w) < 0.2f) + u = p->u *= approx_expf(w); else u = p->u *= expf(w); @@ -216,7 +218,7 @@ __attribute__((always_inline)) * @param half_dt The half time-step for this kick */ __attribute__((always_inline)) INLINE static void hydro_kick_extra( - struct part* p, struct xpart* xp, float dt, float half_dt) { } + struct part* p, struct xpart* xp, float dt, float half_dt) {} /** * @brief Converts hydro quantity of a particle diff --git a/src/hydro/Default/hydro_debug.h b/src/hydro/Default/hydro_debug.h index 5b1648e222c36dc142362f26ad5188a2af397192..2e7c3d683aa18eee7a550ac89517e3bd01e42107 100644 --- a/src/hydro/Default/hydro_debug.h +++ b/src/hydro/Default/hydro_debug.h @@ -23,9 +23,9 @@ __attribute__((always_inline)) "x=[%.3e,%.3e,%.3e], " "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n " "h=%.3e, " - "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, t_begin=%.3e, t_end=%.3e\n", + "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, t_begin=%d, t_end=%d\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->t_begin, - p->t_end); + p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin, + p->ti_end); } diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h index b84fc2700b353e1aadc597b6fa7fe1e30676a26d..4dbc0388ac20cd7f354007e2e80aa6fccccdc5ad 100644 --- a/src/hydro/Default/hydro_part.h +++ b/src/hydro/Default/hydro_part.h @@ -50,10 +50,10 @@ struct part { float h_dt; /* Particle time of beginning of time-step. */ - float t_begin; + int ti_begin; /* Particle time of end of time-step. */ - float t_end; + int ti_end; /* Particle internal energy. */ float u; diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index b135bed20c2642d18e4d1ec121869865c84e87ad..8cc553363122099c748e3e3e1941611e986c8581 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -81,10 +81,10 @@ __attribute__((always_inline)) * and add the self-contribution term. * * @param p The particle to act upon - * @param time The current time + * @param ti_current The current time (on the integer timeline) */ __attribute__((always_inline)) - INLINE static void hydro_end_density(struct part* p, float time) { + INLINE static void hydro_end_density(struct part* p, int ti_current) { /* Some smoothing length multiples. */ const float h = p->h; diff --git a/src/runner.c b/src/runner.c index 4311cc2b0b81e47aee8cf9517aa7b1dd07225b4e..202a62d91672109030c65c4dda4e2554e3ef2da1 100644 --- a/src/runner.c +++ b/src/runner.c @@ -583,7 +583,7 @@ void runner_doghost(struct runner *r, struct cell *c) { if (p->ti_end <= ti_current) { /* Finish the density calculation */ - hydro_end_density(p, ti_current); // MATTHIEU + hydro_end_density(p, ti_current); /* If no derivative, double the smoothing length. */ if (p->density.wcount_dh == 0.0f) h_corr = p->h; @@ -800,7 +800,6 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) { void runner_dokick(struct runner *r, struct cell *c, int timer) { - // const float dt_max_timeline = r->e->timeEnd - r->e->timeBegin; const float global_dt_min = r->e->dt_min; const float global_dt_max = r->e->dt_max; const float ti_current = r->e->ti_current;