From abdb51c96f3092fcf35ac69a01cce308053df04a Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <bert.vandenbroucke@ugent.be>
Date: Sun, 14 Aug 2016 12:01:26 +0100
Subject: [PATCH] Tried implementing primitive variable drift. Does not work
(yet).
---
src/engine.c | 2 +-
src/hydro/Gizmo/hydro.h | 21 +++++++++++++++++----
src/hydro/Gizmo/hydro_iact.h | 14 ++++++++++++++
src/hydro/Gizmo/hydro_part.h | 2 +-
4 files changed, 33 insertions(+), 6 deletions(-)
diff --git a/src/engine.c b/src/engine.c
index 7d50272372..4a91c6be94 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2720,7 +2720,7 @@ void engine_step(struct engine *e) {
mask |= 1 << task_type_sub_pair;
mask |= 1 << task_type_ghost;
mask |= 1 << task_type_extra_ghost; /* Adding unnecessary things to the mask
- does not harm */
+ does not harm */
submask |= 1 << task_subtype_density;
submask |= 1 << task_subtype_gradient;
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index b90632b986..60f7cf55f3 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -19,6 +19,7 @@
#include <float.h>
#include "adiabatic_index.h"
+#include "approx_math.h"
#include "hydro_gradients.h"
/**
@@ -219,6 +220,9 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
p->a_hydro[0] = 0.0f;
p->a_hydro[1] = 0.0f;
p->a_hydro[2] = 0.0f;
+
+ /* Reset the time derivatives. */
+ p->force.h_dt = 0.0f;
}
/**
@@ -273,10 +277,17 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
return;
float dt = (t1 - t0) * timeBase;
-
- p->primitives.rho =
- (p->conserved.mass + p->conserved.flux.mass * dt / p->force.dt) /
- p->geometry.volume;
+ float h_inv = 1.0f / p->h;
+ float w = -hydro_dimension * p->force.h_dt * h_inv * dt;
+
+ // p->primitives.rho =
+ // (p->conserved.mass + p->conserved.flux.mass * dt / p->force.dt) /
+ // p->geometry.volume;
+ if (fabsf(w) < 0.2f) {
+ p->primitives.rho *= approx_expf(w);
+ } else {
+ p->primitives.rho *= expf(w);
+ }
p->primitives.v[0] += p->a_hydro[0] * dt;
p->primitives.v[1] += p->a_hydro[1] * dt;
p->primitives.v[2] += p->a_hydro[2] * dt;
@@ -299,6 +310,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part* p) {
+ p->force.h_dt *= p->h * hydro_dimension_inv;
+
/* Set the hydro acceleration, based on the new momentum and mass */
/* NOTE: the momentum and mass are only correct for active particles, since
only active particles have received flux contributions from all their
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index f19082547d..a615c080c4 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -265,6 +265,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
+ /* Compute h_dt */
+ // float dvdr = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
+ // (Wi[3]-Wj[3])*dx[2];
+ // float ri = 1.0f/r;
+ // float hidp1 = pow_dimension_plus_one(hi_inv);
+ // float hjdp1 = pow_dimension_plus_one(hj_inv);
+ // float wi_dr = hidp1 * wi_dx;
+ // float wj_dr = hjdp1 * wj_dx;
+ // dvdr *= ri;
+ // pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
+ // if(mode == 1){
+ // pj->force.h_dt -= pi->conserved.mass * dvdr / pi->primitives.rho * wj_dr;
+ //}
+
/* Compute area */
/* eqn. (7) */
Anorm = 0.0f;
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index ae5be8852f..9c7b0dfe5b 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/Gizmo/hydro_part.h
@@ -167,7 +167,7 @@ struct part {
/* Quantities used during the force loop. */
struct {
- /* Needed to compile the code. */
+ /* Needed to drift the primitive variables. */
float h_dt;
/* Physical time step of the particle. */
--
GitLab