diff --git a/src/engine.c b/src/engine.c
index a8ec2d1ddaf9220c7b2553fb223c968c53768eec..ac4ac50d49f244646c762e8efa6534cf9cb407f5 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 3dc8b62918c9df34cbfbbd0c7f4f34a7c524b185..a4d9e58087879645d5c98912ec2039015bcc145a 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -202,6 +202,13 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
struct part* p) {
hydro_gradients_finalize(p);
+
+ /* Initialize fluxes */
+ p->conserved.flux.mass = 0.0f;
+ p->conserved.flux.momentum[0] = 0.0f;
+ p->conserved.flux.momentum[1] = 0.0f;
+ p->conserved.flux.momentum[2] = 0.0f;
+ p->conserved.flux.energy = 0.0f;
}
/**
@@ -269,7 +276,20 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
* @param timeBase Conversion factor between integer and physical time.
*/
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
- struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {}
+ struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {
+
+ return;
+ float dt = (t1 - t0) * timeBase;
+
+ p->primitives.rho =
+ (p->conserved.mass + p->conserved.flux.mass * dt / p->force.dt) /
+ p->geometry.volume;
+ 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;
+ float u = p->conserved.energy + p->conserved.flux.energy * dt / p->force.dt;
+ p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u;
+}
/**
* @brief Set the particle acceleration after the flux loop
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index 1d5fdf15322767432c87db8d35f63b99695961bb..d4e7bcb6bdefebdf279ddbdba526b238a7c6cfa3 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -250,6 +250,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* The flux will be exchanged using the smallest time step of the two
* particles */
mindt = fminf(dti, dtj);
+ // dti = mindt;
+ // dtj = mindt;
/* Compute kernel of pi. */
hi_inv = 1.0 / hi;
@@ -363,44 +365,69 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* Update conserved variables */
/* eqn. (16) */
- pi->conserved.mass -= pi->force.dt * Anorm * totflux[0];
- pi->conserved.momentum[0] -= pi->force.dt * Anorm * totflux[1];
- pi->conserved.momentum[1] -= pi->force.dt * Anorm * totflux[2];
- pi->conserved.momentum[2] -= pi->force.dt * Anorm * totflux[3];
- pi->conserved.energy -= pi->force.dt * Anorm * totflux[4];
+ pi->conserved.mass -= dti * Anorm * totflux[0];
+ pi->conserved.momentum[0] -= dti * Anorm * totflux[1];
+ pi->conserved.momentum[1] -= dti * Anorm * totflux[2];
+ pi->conserved.momentum[2] -= dti * Anorm * totflux[3];
+ pi->conserved.energy -= dti * Anorm * totflux[4];
+
+ pi->conserved.flux.mass -= dti * Anorm * totflux[0];
+ pi->conserved.flux.momentum[0] -= dti * Anorm * totflux[1];
+ pi->conserved.flux.momentum[1] -= dti * Anorm * totflux[2];
+ pi->conserved.flux.momentum[2] -= dti * Anorm * totflux[3];
+ pi->conserved.flux.energy -= dti * Anorm * totflux[4];
float ekin = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
pi->primitives.v[1] * pi->primitives.v[1] +
pi->primitives.v[2] * pi->primitives.v[2]);
- pi->conserved.energy +=
- pi->force.dt * Anorm * totflux[1] * pi->primitives.v[0];
- pi->conserved.energy +=
- pi->force.dt * Anorm * totflux[2] * pi->primitives.v[1];
- pi->conserved.energy +=
- pi->force.dt * Anorm * totflux[3] * pi->primitives.v[2];
- pi->conserved.energy -= pi->force.dt * Anorm * totflux[0] * ekin;
-
- /* the non symmetric version is never called when using mindt, whether this
- * piece of code
- * should always be executed or only in the symmetric case is currently
- * unclear */
+ pi->conserved.energy += dti * Anorm * totflux[1] * pi->primitives.v[0];
+ pi->conserved.energy += dti * Anorm * totflux[2] * pi->primitives.v[1];
+ pi->conserved.energy += dti * Anorm * totflux[3] * pi->primitives.v[2];
+ pi->conserved.energy -= dti * Anorm * totflux[0] * ekin;
+
+ pi->conserved.flux.energy += dti * Anorm * totflux[1] * pi->primitives.v[0];
+ pi->conserved.flux.energy += dti * Anorm * totflux[2] * pi->primitives.v[1];
+ pi->conserved.flux.energy += dti * Anorm * totflux[3] * pi->primitives.v[2];
+ pi->conserved.flux.energy -= dti * Anorm * totflux[0] * ekin;
+
+ /* here is how it works:
+ Mode will only be 1 if both particles are ACTIVE and they are in the same
+ cell. In this case, this method IS the flux calculation for particle j, and
+ we HAVE TO UPDATE it.
+ Mode 0 can mean several things: it can mean that particle j is INACTIVE, in
+ which case we NEED TO UPDATE it, since otherwise the flux is lost from the
+ system and the conserved variable is not conserved.
+ It can also mean that particle j sits in another cell and is ACTIVE. In
+ this case, the flux exchange for particle j is done TWICE and we SHOULD NOT
+ UPDATE particle j.
+ ==> we update particle j if (MODE IS 1) OR (j IS INACTIVE)
+ */
+ // if (mode == 1 || pj->ti_end > pi->ti_end) {
if (mode == 1) {
- pj->conserved.mass += pj->force.dt * Anorm * totflux[0];
- pj->conserved.momentum[0] += pj->force.dt * Anorm * totflux[1];
- pj->conserved.momentum[1] += pj->force.dt * Anorm * totflux[2];
- pj->conserved.momentum[2] += pj->force.dt * Anorm * totflux[3];
- pj->conserved.energy += pj->force.dt * Anorm * totflux[4];
+ pj->conserved.mass += dtj * Anorm * totflux[0];
+ pj->conserved.momentum[0] += dtj * Anorm * totflux[1];
+ pj->conserved.momentum[1] += dtj * Anorm * totflux[2];
+ pj->conserved.momentum[2] += dtj * Anorm * totflux[3];
+ pj->conserved.energy += dtj * Anorm * totflux[4];
+
+ pj->conserved.flux.mass += dtj * Anorm * totflux[0];
+ pj->conserved.flux.momentum[0] += dtj * Anorm * totflux[1];
+ pj->conserved.flux.momentum[1] += dtj * Anorm * totflux[2];
+ pj->conserved.flux.momentum[2] += dtj * Anorm * totflux[3];
+ pj->conserved.flux.energy += dtj * Anorm * totflux[4];
ekin = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
pj->primitives.v[1] * pj->primitives.v[1] +
pj->primitives.v[2] * pj->primitives.v[2]);
- pj->conserved.energy -=
- pj->force.dt * Anorm * totflux[1] * pj->primitives.v[0];
- pj->conserved.energy -=
- pj->force.dt * Anorm * totflux[2] * pj->primitives.v[1];
- pj->conserved.energy -=
- pj->force.dt * Anorm * totflux[3] * pj->primitives.v[2];
- pj->conserved.energy += pj->force.dt * Anorm * totflux[0] * ekin;
+ pj->conserved.energy -= dtj * Anorm * totflux[1] * pj->primitives.v[0];
+ pj->conserved.energy -= dtj * Anorm * totflux[2] * pj->primitives.v[1];
+ pj->conserved.energy -= dtj * Anorm * totflux[3] * pj->primitives.v[2];
+ pj->conserved.energy += dtj * Anorm * totflux[0] * ekin;
+
+ pj->conserved.flux.energy -= dtj * Anorm * totflux[1] * pj->primitives.v[0];
+ pj->conserved.flux.energy -= dtj * Anorm * totflux[2] * pj->primitives.v[1];
+ pj->conserved.flux.energy -= dtj * Anorm * totflux[3] * pj->primitives.v[2];
+ pj->conserved.flux.energy += dtj * Anorm * totflux[0] * ekin;
}
}
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index ce1d7a5e00568c69acd1d53ac429dd35d1c56900..ae5be8852fa1965372e808c111f459d0b939ead6 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/Gizmo/hydro_part.h
@@ -111,6 +111,20 @@ struct part {
/* Fluid total energy (= internal energy + fluid kinetic energy). */
float energy;
+ /* Fluxes. */
+ struct {
+
+ /* Mass flux. */
+ float mass;
+
+ /* Momentum flux. */
+ float momentum[3];
+
+ /* Energy flux. */
+ float energy;
+
+ } flux;
+
} conserved;
/* Geometrical quantities used for hydro. */