From 2b522851a7bfa9f6ab858767f3ac8265b10ea2ee Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Tue, 9 Feb 2016 20:46:00 +0000
Subject: [PATCH] Ready to work on dh/dt
---
src/hydro/Gadget2/hydro.h | 6 ++---
src/hydro/Gadget2/hydro_part.h | 7 ++++++
src/runner.c | 41 ++++++++++++++++------------------
3 files changed, 28 insertions(+), 26 deletions(-)
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 7ac07f499e..53366bd84d 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -186,10 +186,8 @@ __attribute__((always_inline))
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part* p, struct xpart* xp, float t0, float t1) {
- const float dt = t1 - t0;
-
- p->rho *= expf(-p->div_v * dt);
- p->h *= expf(0.33333333f * p->div_v * dt);
+ // p->rho *= expf(-p->div_v * dt);
+ // p->h *= expf(0.33333333f * p->div_v * dt)
const float dt_entr = t1 - 0.5f * (p->t_begin + p->t_end);
p->pressure =
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 8ad6a6ef24..45c04ba75c 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -68,6 +68,13 @@ struct part {
} density;
+ struct {
+
+ /* Time derivative of the smoothing length */
+ float h_dt;
+
+ } force;
+
/* Particle entropy. */
float entropy;
diff --git a/src/runner.c b/src/runner.c
index bfdaba3c64..66c1e02deb 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -34,6 +34,7 @@
#include "runner.h"
/* Local headers. */
+#include "approx_math.h"
#include "atomic.h"
#include "const.h"
#include "debug.h"
@@ -705,6 +706,7 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
struct part *restrict p, *restrict parts = c->parts;
struct xpart *restrict xp, *restrict xparts = c->xparts;
float dx_max = 0.f, h_max = 0.f;
+ float w;
TIMER_TIC
@@ -718,37 +720,32 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
p = &parts[k];
xp = &xparts[k];
+ /* Useful quantity */
+ const float h_inv = 1.0f / p->h;
+
/* Drift... */
p->x[0] += xp->v_full[0] * dt;
p->x[1] += xp->v_full[1] * dt;
p->x[2] += xp->v_full[2] * dt;
- /* Predict velocities */
+ /* Predict velocities (for hydro terms) */
p->v[0] += p->a[0] * dt;
p->v[1] += p->a[1] * dt;
p->v[2] += p->a[2] * dt;
- /* /\* Predict smoothing length *\/ */
- /* w = p->force.h_dt * ih * dt; */
- /* if (fabsf(w) < 0.01f) /\* 1st order expansion of exp(w) *\/ */
- /* p->h *= */
- /* 1.0f + */
- /* w * (1.0f + w * (0.5f + w * (1.0f / 6.0f + 1.0f / 24.0f *
- * w))); */
- /* else */
- /* p->h *= expf(w); */
-
- // MATTHIEU
-
- /* /\* Predict density *\/ */
- /* w = -3.0f * p->force.h_dt * ih * dt; */
- /* if (fabsf(w) < 0.1f) */
- /* p->rho *= */
- /* 1.0f + */
- /* w * (1.0f + w * (0.5f + w * (1.0f / 6.0f + 1.0f / 24.0f *
- * w))); */
- /* else */
- /* p->rho *= expf(w); */
+ /* Predict smoothing length */
+ w = p->force.h_dt * h_inv * dt;
+ if (fabsf(w) < 0.2f)
+ p->h *= approx_expf(w); /* 4th order expansion of exp(w) */
+ else
+ p->h *= expf(w);
+
+ /* Predict density */
+ w = -3.0f * p->force.h_dt * h_inv * dt;
+ if (fabsf(w) < 0.2f)
+ p->rho *= approx_expf(w); /* 4th order expansion of exp(w) */
+ else
+ p->rho *= expf(w);
/* Predict the values of the extra fields */
hydro_predict_extra(p, xp, r->e->timeOld, r->e->time);
--
GitLab