Commit 095b5eda authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Default SPH also moved to integer timeline

parent 9ec53b11
......@@ -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
......
......@@ -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);
}
......@@ -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;
......
......@@ -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;
......
......@@ -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;
......
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