From 7cc5943f6e16b0f8b01704afcbe65bb4e228c923 Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <bert.vandenbroucke@ugent.be>
Date: Fri, 28 Jul 2017 11:37:18 +0100
Subject: [PATCH] Partially restored old time step criterion. Isolated flux
limiter into a separate file and made it configurable. Removed some
superfluous gravity variables.
---
src/Makefile.am | 1 +
src/const.h | 3 +
src/hydro/Gizmo/hydro.h | 18 +--
src/hydro/Gizmo/hydro_debug.h | 4 +-
src/hydro/Gizmo/hydro_flux_limiters.h | 81 +++++++++++
src/hydro/Gizmo/hydro_gradients.h | 37 -----
src/hydro/Gizmo/hydro_gradients_gizmo.h | 185 ------------------------
src/hydro/Gizmo/hydro_iact.h | 42 +-----
src/hydro/Gizmo/hydro_io.h | 9 +-
src/hydro/Gizmo/hydro_part.h | 15 +-
10 files changed, 114 insertions(+), 281 deletions(-)
create mode 100644 src/hydro/Gizmo/hydro_flux_limiters.h
diff --git a/src/Makefile.am b/src/Makefile.am
index 2ddcdb0908..3784886278 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -86,6 +86,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
hydro/Gizmo/hydro_slope_limiters_cell.h \
hydro/Gizmo/hydro_slope_limiters_face.h \
hydro/Gizmo/hydro_slope_limiters.h \
+ hydro/Gizmo/hydro_flux_limiters.h \
hydro/Gizmo/hydro_unphysical.h \
hydro/Gizmo/hydro_velocities.h \
hydro/Shadowswift/hydro_debug.h \
diff --git a/src/const.h b/src/const.h
index f5232de0b2..c8060a2be5 100644
--- a/src/const.h
+++ b/src/const.h
@@ -49,6 +49,9 @@
#define SLOPE_LIMITER_PER_FACE
#define SLOPE_LIMITER_CELL_WIDE
+/* Types of flux limiter to use (GIZMO_SPH only) */
+#define GIZMO_FLUX_LIMITER
+
/* Options to control the movement of particles for GIZMO_SPH. */
/* This option disables particle movement */
//#define GIZMO_FIX_PARTICLES
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 301700379b..2c2f54699b 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -53,12 +53,17 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
vrel[0] = p->primitives.v[0] - xp->v_full[0];
vrel[1] = p->primitives.v[1] - xp->v_full[1];
vrel[2] = p->primitives.v[2] - xp->v_full[2];
- const float vmax =
+ float vmax =
sqrtf(vrel[0] * vrel[0] + vrel[1] * vrel[1] + vrel[2] * vrel[2]) +
sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+ vmax = max(vmax, p->timestepvars.vmax);
const float psize = powf(p->geometry.volume / hydro_dimension_unit_sphere,
hydro_dimension_inv);
- return CFL_condition * min(psize / vmax, p->timestepvars.dt_min);
+ float dt = FLT_MAX;
+ if (vmax > 0.) {
+ dt = psize / vmax;
+ }
+ return CFL_condition * dt;
}
/**
@@ -420,7 +425,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* restrict p, struct xpart* restrict xp) {
/* Initialize time step criterion variables */
- p->timestepvars.dt_min = FLT_MAX;
+ p->timestepvars.vmax = 0.;
/* Set the actual velocity of the particle */
hydro_velocities_prepare_force(p, xp);
@@ -637,12 +642,6 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
a_grav[1] = p->gpart->a_grav[1];
a_grav[2] = p->gpart->a_grav[2];
- /* Store the gravitational acceleration for later use. */
- /* This is used for the prediction step. */
- p->gravity.old_a[0] = a_grav[0];
- p->gravity.old_a[1] = a_grav[1];
- p->gravity.old_a[2] = a_grav[2];
-
/* Make sure the gpart knows the mass has changed. */
p->gpart->mass = p->conserved.mass;
@@ -664,7 +663,6 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0];
p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
-
}
/* reset fluxes */
diff --git a/src/hydro/Gizmo/hydro_debug.h b/src/hydro/Gizmo/hydro_debug.h
index d22ff03b63..17e7f8a085 100644
--- a/src/hydro/Gizmo/hydro_debug.h
+++ b/src/hydro/Gizmo/hydro_debug.h
@@ -46,7 +46,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"volume=%.3e, "
"matrix_E=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]]}, "
"timestepvars={"
- "dt_min=%.3e},"
+ "vmax=%.3e},"
"density={"
"div_v=%.3e, "
"wcount_dh=%.3e, "
@@ -75,7 +75,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
p->geometry.matrix_E[1][0], p->geometry.matrix_E[1][1],
p->geometry.matrix_E[1][2], p->geometry.matrix_E[2][0],
p->geometry.matrix_E[2][1], p->geometry.matrix_E[2][2],
- p->timestepvars.dt_min, p->density.div_v, p->density.wcount_dh,
+ p->timestepvars.vmax, p->density.div_v, p->density.wcount_dh,
p->density.curl_v[0], p->density.curl_v[1], p->density.curl_v[2],
p->density.wcount);
}
diff --git a/src/hydro/Gizmo/hydro_flux_limiters.h b/src/hydro/Gizmo/hydro_flux_limiters.h
new file mode 100644
index 0000000000..dc91cf2808
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_flux_limiters.h
@@ -0,0 +1,81 @@
+
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_HYDRO_FLUX_LIMITERS_H
+#define SWIFT_HYDRO_FLUX_LIMITERS_H
+
+#ifdef GIZMO_FLUX_LIMITER
+
+#define HYDRO_FLUX_LIMITER_IMPLEMENTATION "GIZMO flux limiter"
+
+/**
+ * @brief Limit the flux between two particles.
+ *
+ * @param flux Unlimited flux between the particles.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_flux_limiters_apply(
+ float* flux, struct part* pi, struct part* pj) {
+
+ float flux_limit_factor = 1.;
+ const float timefac = max(pi->force.dt, pj->force.dt);
+ const float areafac = max(pi->geometry.Atot, pj->geometry.Atot);
+ const float totfac = timefac * areafac;
+ if (flux[0] * totfac > pi->conserved.mass) {
+ flux_limit_factor = pi->conserved.mass / (flux[0] * totfac);
+ }
+ if (flux[0] * totfac > pj->conserved.mass) {
+ flux_limit_factor =
+ min(pj->conserved.mass / (flux[0] * totfac), flux_limit_factor);
+ }
+ if (flux[4] * totfac > pi->conserved.energy) {
+ flux_limit_factor =
+ min(pi->conserved.energy / (flux[4] * totfac), flux_limit_factor);
+ }
+ if (flux[4] * totfac > pj->conserved.energy) {
+ flux_limit_factor =
+ min(pj->conserved.energy / (flux[4] * totfac), flux_limit_factor);
+ }
+
+ flux[0] *= flux_limit_factor;
+ flux[1] *= flux_limit_factor;
+ flux[2] *= flux_limit_factor;
+ flux[3] *= flux_limit_factor;
+ flux[4] *= flux_limit_factor;
+}
+
+#else
+
+#define HYDRO_FLUX_LIMITER_IMPLEMENTATION "No flux limiter"
+
+/**
+ * @brief Limit the flux between two particles.
+ *
+ * @param flux Unlimited flux between the particles.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_flux_limiters_apply(
+ float* flux, struct part* pi, struct part* pj) {}
+
+#endif
+
+#endif // SWIFT_HYDRO_FLUX_LIMITERS_H
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index 5ad6d87619..896128bd45 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -99,7 +99,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
float xij_j[3];
int k;
float xfac;
- float a_grav_i[3], a_grav_j[3];
/* perform gradient reconstruction in space and time */
/* space */
@@ -141,34 +140,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
pj->primitives.gradients.P[1] * xij_j[1] +
pj->primitives.gradients.P[2] * xij_j[2];
- a_grav_i[0] = pi->gravity.old_a[0];
- a_grav_i[1] = pi->gravity.old_a[1];
- a_grav_i[2] = pi->gravity.old_a[2];
-
- a_grav_i[0] += pi->gravity.grad_a[0][0] * xij_i[0] +
- pi->gravity.grad_a[0][1] * xij_i[1] +
- pi->gravity.grad_a[0][2] * xij_i[2];
- a_grav_i[1] += pi->gravity.grad_a[1][0] * xij_i[0] +
- pi->gravity.grad_a[1][1] * xij_i[1] +
- pi->gravity.grad_a[1][2] * xij_i[2];
- a_grav_i[2] += pi->gravity.grad_a[2][0] * xij_i[0] +
- pi->gravity.grad_a[2][1] * xij_i[1] +
- pi->gravity.grad_a[2][2] * xij_i[2];
-
- a_grav_j[0] = pj->gravity.old_a[0];
- a_grav_j[1] = pj->gravity.old_a[1];
- a_grav_j[2] = pj->gravity.old_a[2];
-
- a_grav_j[0] += pj->gravity.grad_a[0][0] * xij_j[0] +
- pj->gravity.grad_a[0][1] * xij_j[1] +
- pj->gravity.grad_a[0][2] * xij_j[2];
- a_grav_j[1] += pj->gravity.grad_a[1][0] * xij_j[0] +
- pj->gravity.grad_a[1][1] * xij_j[1] +
- pj->gravity.grad_a[1][2] * xij_j[2];
- a_grav_j[2] += pj->gravity.grad_a[2][0] * xij_j[0] +
- pj->gravity.grad_a[2][1] * xij_j[1] +
- pj->gravity.grad_a[2][2] * xij_j[2];
-
hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
/* time */
@@ -198,10 +169,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
pi->primitives.gradients.v[1][1] +
pi->primitives.gradients.v[2][2]));
-
- dWi[1] += 0.5 * mindt * a_grav_i[0];
- dWi[2] += 0.5 * mindt * a_grav_i[1];
- dWi[3] += 0.5 * mindt * a_grav_i[2];
}
if (Wj[0] > 0.0f) {
@@ -230,10 +197,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
pj->primitives.gradients.v[1][1] +
pj->primitives.gradients.v[2][2]));
-
- dWj[1] += 0.5 * mindt * a_grav_j[0];
- dWj[2] += 0.5 * mindt * a_grav_j[1];
- dWj[3] += 0.5 * mindt * a_grav_j[2];
}
Wi[0] += dWi[0];
diff --git a/src/hydro/Gizmo/hydro_gradients_gizmo.h b/src/hydro/Gizmo/hydro_gradients_gizmo.h
index ee3ad6919f..bc50c10d84 100644
--- a/src/hydro/Gizmo/hydro_gradients_gizmo.h
+++ b/src/hydro/Gizmo/hydro_gradients_gizmo.h
@@ -45,18 +45,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_init(
p->primitives.gradients.P[1] = 0.0f;
p->primitives.gradients.P[2] = 0.0f;
- p->gravity.grad_a[0][0] = 0.0f;
- p->gravity.grad_a[0][1] = 0.0f;
- p->gravity.grad_a[0][2] = 0.0f;
-
- p->gravity.grad_a[1][0] = 0.0f;
- p->gravity.grad_a[1][1] = 0.0f;
- p->gravity.grad_a[1][2] = 0.0f;
-
- p->gravity.grad_a[2][0] = 0.0f;
- p->gravity.grad_a[2][1] = 0.0f;
- p->gravity.grad_a[2][2] = 0.0f;
-
hydro_slope_limit_cell_init(p);
}
@@ -157,35 +145,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
(Wi[4] - Wj[4]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->gravity.grad_a[0][0] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->gravity.grad_a[0][1] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->gravity.grad_a[0][2] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
- pi->gravity.grad_a[1][0] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->gravity.grad_a[1][1] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->gravity.grad_a[1][2] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
- pi->gravity.grad_a[2][0] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->gravity.grad_a[2][1] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->gravity.grad_a[2][2] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
} else {
/* The gradient matrix was not well-behaved, switch to SPH gradients */
@@ -223,27 +182,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
pi->primitives.gradients.P[2] -=
wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
-
- pi->gravity.grad_a[0][0] -=
- wi_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
- pi->gravity.grad_a[0][1] -=
- wi_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
- pi->gravity.grad_a[0][2] -=
- wi_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-
- pi->gravity.grad_a[1][0] -=
- wi_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
- pi->gravity.grad_a[1][1] -=
- wi_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
- pi->gravity.grad_a[1][2] -=
- wi_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-
- pi->gravity.grad_a[2][0] -=
- wi_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
- pi->gravity.grad_a[2][1] -=
- wi_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
- pi->gravity.grad_a[2][2] -=
- wi_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
}
hydro_slope_limit_cell_collect(pi, pj, r);
@@ -306,35 +244,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
(Wi[4] - Wj[4]) * wj *
(Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
- pj->gravity.grad_a[0][0] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->gravity.grad_a[0][1] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->gravity.grad_a[0][2] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-
- pj->gravity.grad_a[1][0] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->gravity.grad_a[1][1] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->gravity.grad_a[1][2] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-
- pj->gravity.grad_a[2][0] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->gravity.grad_a[2][1] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->gravity.grad_a[2][2] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
} else {
/* SPH gradients */
@@ -371,27 +280,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
pj->primitives.gradients.P[2] -=
wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
-
- pj->gravity.grad_a[0][0] -=
- wj_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
- pj->gravity.grad_a[0][1] -=
- wj_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
- pj->gravity.grad_a[0][2] -=
- wj_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-
- pj->gravity.grad_a[1][0] -=
- wj_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
- pj->gravity.grad_a[1][1] -=
- wj_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
- pj->gravity.grad_a[1][2] -=
- wj_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-
- pj->gravity.grad_a[2][0] -=
- wj_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
- pj->gravity.grad_a[2][1] -=
- wj_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
- pj->gravity.grad_a[2][2] -=
- wj_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
}
hydro_slope_limit_cell_collect(pj, pi, r);
@@ -493,35 +381,6 @@ hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
(Wi[4] - Wj[4]) * wi *
(Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->gravity.grad_a[0][0] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->gravity.grad_a[0][1] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->gravity.grad_a[0][2] +=
- (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
- pi->gravity.grad_a[1][0] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->gravity.grad_a[1][1] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->gravity.grad_a[1][2] +=
- (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
- pi->gravity.grad_a[2][0] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->gravity.grad_a[2][1] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->gravity.grad_a[2][2] +=
- (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
} else {
/* Gradient matrix is not well-behaved, switch to SPH gradients */
@@ -558,27 +417,6 @@ hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
pi->primitives.gradients.P[2] -=
wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
-
- pi->gravity.grad_a[0][0] -=
- wi_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
- pi->gravity.grad_a[0][1] -=
- wi_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
- pi->gravity.grad_a[0][2] -=
- wi_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
-
- pi->gravity.grad_a[1][0] -=
- wi_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
- pi->gravity.grad_a[1][1] -=
- wi_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
- pi->gravity.grad_a[1][2] -=
- wi_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
-
- pi->gravity.grad_a[2][0] -=
- wi_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
- pi->gravity.grad_a[2][1] -=
- wi_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
- pi->gravity.grad_a[2][2] -=
- wi_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
}
hydro_slope_limit_cell_collect(pi, pj, r);
@@ -618,17 +456,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
p->primitives.gradients.P[1] *= ihdim;
p->primitives.gradients.P[2] *= ihdim;
- p->gravity.grad_a[0][0] *= ihdim;
- p->gravity.grad_a[0][1] *= ihdim;
- p->gravity.grad_a[0][2] *= ihdim;
-
- p->gravity.grad_a[1][0] *= ihdim;
- p->gravity.grad_a[1][1] *= ihdim;
- p->gravity.grad_a[1][2] *= ihdim;
-
- p->gravity.grad_a[2][0] *= ihdim;
- p->gravity.grad_a[2][1] *= ihdim;
- p->gravity.grad_a[2][2] *= ihdim;
} else {
const float ihdimp1 = pow_dimension_plus_one(ih);
@@ -653,18 +480,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
p->primitives.gradients.P[0] *= ihdimp1 * volume;
p->primitives.gradients.P[1] *= ihdimp1 * volume;
p->primitives.gradients.P[2] *= ihdimp1 * volume;
-
- p->gravity.grad_a[0][0] *= ihdimp1 * volume;
- p->gravity.grad_a[0][1] *= ihdimp1 * volume;
- p->gravity.grad_a[0][2] *= ihdimp1 * volume;
-
- p->gravity.grad_a[1][0] *= ihdimp1 * volume;
- p->gravity.grad_a[1][1] *= ihdimp1 * volume;
- p->gravity.grad_a[1][2] *= ihdimp1 * volume;
-
- p->gravity.grad_a[2][0] *= ihdimp1 * volume;
- p->gravity.grad_a[2][1] *= ihdimp1 * volume;
- p->gravity.grad_a[2][2] *= ihdimp1 * volume;
}
hydro_slope_limit_cell(p);
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index b2eb24a582..0c7c8251b7 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -20,6 +20,7 @@
******************************************************************************/
#include "adiabatic_index.h"
+#include "hydro_flux_limiters.h"
#include "hydro_gradients.h"
#include "riemann.h"
@@ -306,7 +307,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float Bj[3][3];
float Vi, Vj;
float xij_i[3], xfac, xijdotdx;
- float vmax, dvdotdx, dt_min;
+ float vmax, dvdotdx;
float vi[3], vj[3], vij[3];
float Wi[5], Wj[5];
float dti, dtj, mindt;
@@ -346,16 +347,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
}
dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
(Wi[3] - Wj[3]) * dx[2];
- dvdotdx = min(dvdotdx,
- (vi[0] - vj[0])*dx[0] + (vi[1]-vj[1])*dx[1] + (vi[2]-vj[2])*dx[2]);
+ dvdotdx = min(dvdotdx, (vi[0] - vj[0]) * dx[0] + (vi[1] - vj[1]) * dx[1] +
+ (vi[2] - vj[2]) * dx[2]);
if (dvdotdx < 0.) {
/* the magical factor 3 also appears in Gadget2 */
- vmax -= 3.*dvdotdx / r;
+ vmax -= 3. * dvdotdx / r;
}
- dt_min = r / vmax;
- pi->timestepvars.dt_min = min(pi->timestepvars.dt_min, dt_min);
+ pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
if (mode == 1) {
- pj->timestepvars.dt_min = min(pj->timestepvars.dt_min, dt_min);
+ pj->timestepvars.vmax = max(pj->timestepvars.vmax, vmax);
}
/* The flux will be exchanged using the smallest time step of the two
@@ -491,33 +491,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float totflux[5];
riemann_solve_for_flux(Wi, Wj, n_unit, vij, totflux);
- /* Flux limiter */
- float flux_limit_factor = 1.;
- float timefac = max(dti, dtj);
- float areafac = max(pi->geometry.Atot, pj->geometry.Atot);
- if (totflux[0] * areafac * timefac > pi->conserved.mass) {
- flux_limit_factor = pi->conserved.mass / (totflux[0] * areafac * timefac);
- }
- if (totflux[0] * areafac * timefac > pj->conserved.mass) {
- flux_limit_factor =
- min(pj->conserved.mass / (totflux[0] * areafac * timefac),
- flux_limit_factor);
- }
- if (totflux[4] * areafac * timefac > pi->conserved.energy) {
- flux_limit_factor =
- min(pi->conserved.energy / (totflux[4] * areafac * timefac),
- flux_limit_factor);
- }
- if (totflux[4] * areafac * timefac > pj->conserved.energy) {
- flux_limit_factor =
- min(pj->conserved.energy / (totflux[4] * areafac * timefac),
- flux_limit_factor);
- }
- totflux[0] *= flux_limit_factor;
- totflux[1] *= flux_limit_factor;
- totflux[2] *= flux_limit_factor;
- totflux[3] *= flux_limit_factor;
- totflux[4] *= flux_limit_factor;
+ hydro_flux_limiters_apply(totflux, pi, pj);
/* Store mass flux */
float mflux = Anorm * totflux[0];
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index 3d58be2f47..d20f7e2eb1 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -18,6 +18,7 @@
******************************************************************************/
#include "adiabatic_index.h"
+#include "hydro_flux_limiters.h"
#include "hydro_gradients.h"
#include "hydro_slope_limiters.h"
#include "io_properties.h"
@@ -127,7 +128,7 @@ float convert_Etot(struct engine* e, struct part* p) {
void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) {
- *num_fields = 11;
+ *num_fields = 10;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -152,8 +153,6 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
list[9] =
io_make_output_field_convert_part("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
parts, conserved.energy, convert_Etot);
- list[10] = io_make_output_field("GravAcceleration", FLOAT, 3,
- UNIT_CONV_ACCELERATION, parts, gravity.old_a);
}
/**
@@ -171,6 +170,10 @@ void writeSPHflavour(hid_t h_grpsph) {
io_write_attribute_s(h_grpsph, "Piecewise slope limiter model",
HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION);
+ /* Flux limiter information */
+ io_write_attribute_s(h_grpsph, "Flux limiter model",
+ HYDRO_FLUX_LIMITER_IMPLEMENTATION);
+
/* Riemann solver information */
io_write_attribute_s(h_grpsph, "Riemann solver type",
RIEMANN_SOLVER_IMPLEMENTATION);
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index e6ed233f0f..47f722c5a2 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/Gizmo/hydro_part.h
@@ -156,8 +156,11 @@ struct part {
/* Variables used for timestep calculation. */
struct {
- /* Minimum non-normalized timestep based on the neighbours. */
- float dt_min;
+ /* Maximum signal velocity among all the neighbours of the particle. The
+ * signal velocity encodes information about the relative fluid velocities
+ * AND particle velocities of the neighbour and this particle, as well as
+ * the sound speed of both particles. */
+ float vmax;
} timestepvars;
@@ -201,14 +204,6 @@ struct part {
/* Specific stuff for the gravity-hydro coupling. */
struct {
- /* Previous value of the gravitational acceleration. */
- float old_a[3];
-
- float grad_a[3][3];
-
- /* Previous value of the mass flux vector. */
- float old_mflux[3];
-
/* Current value of the mass flux vector. */
float mflux[3];
--
GitLab