diff --git a/src/Makefile.am b/src/Makefile.am
index e229d7f11e76aae4b6898a2f86145ec2bd025592..b820e8060778fdf8e7313d1b47596e370927d274 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -86,7 +86,6 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.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/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 6273f5009a051e17d2a2954ac994c3f189110464..99cbfc206a603fda23478f247dd20ba1790df3a7 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -1,4 +1,3 @@
-
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
@@ -19,7 +18,6 @@
*
******************************************************************************/
-#include <float.h>
#include "adiabatic_index.h"
#include "approx_math.h"
#include "equation_of_state.h"
@@ -31,6 +29,8 @@
#include "minmax.h"
#include "riemann.h"
+#include <float.h>
+
//#define GIZMO_LLOYD_ITERATION
/**
@@ -50,6 +50,9 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return CFL_condition;
#endif
+ /* v_full is the actual velocity of the particle, primitives.v is its
+ hydrodynamical velocity. The time step depends on the relative difference
+ of the two. */
float vrel[3];
vrel[0] = p->primitives.v[0] - xp->v_full[0];
vrel[1] = p->primitives.v[1] - xp->v_full[1];
@@ -71,12 +74,8 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* @brief Does some extra hydro operations once the actual physical time step
* for the particle is known.
*
- * We use this to store the physical time step, since it is used for the flux
- * exchange during the force loop.
- *
- * We also set the active flag of the particle to inactive. It will be set to
- * active in hydro_init_part, which is called the next time the particle becomes
- * active.
+ * This method is no longer used, as Gizmo is now unaware of the actual particle
+ * time step.
*
* @param p The particle to act upon.
* @param dt Physical time step of the particle during the next step.
@@ -93,9 +92,6 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
error("NaN time step assigned to particle!");
}
#endif
-
- p->force.dt = dt;
- p->force.active = 0;
}
/**
@@ -167,6 +163,11 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
/* initialize the particle velocity based on the primitive fluid velocity */
hydro_velocities_init(p, xp);
+ /* ignore accelerations present in the initial condition */
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+
/* we cannot initialize wcorr in init_part, as init_part gets called every
time the density loop is repeated, and the whole point of storing wcorr
is to have a way of remembering that we need more neighbours for this
@@ -200,10 +201,6 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->geometry.centroid[0] = 0.0f;
p->geometry.centroid[1] = 0.0f;
p->geometry.centroid[2] = 0.0f;
- p->geometry.Atot = 0.0f;
-
- /* Set the active flag to active. */
- p->force.active = 1;
}
/**
@@ -406,7 +403,6 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
p->geometry.centroid[0] = 0.0f;
p->geometry.centroid[1] = 0.0f;
p->geometry.centroid[2] = 0.0f;
- p->geometry.Atot = 1.0f;
}
/**
@@ -452,6 +448,12 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
p->gravity.mflux[1] = 0.0f;
p->gravity.mflux[2] = 0.0f;
+ 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;
+
#ifdef GIZMO_LLOYD_ITERATION
/* reset the gradients to zero, as we don't want them */
hydro_gradients_init(p);
@@ -534,27 +536,20 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->h *= h_corr;
}
-/* we temporarily disabled the primitive variable drift.
- This should be reenabled once gravity works, and we have time to check that
- drifting works properly. */
-// const float w2 = -hydro_dimension * w1;
-// if (fabsf(w2) < 0.2f) {
-// p->primitives.rho *= approx_expf(w2);
-// } else {
-// p->primitives.rho *= expf(w2);
-// }
-
-// p->primitives.v[0] += (p->a_hydro[0] + p->gravity.old_a[0]) * dt;
-// p->primitives.v[1] += (p->a_hydro[1] + p->gravity.old_a[1]) * dt;
-// p->primitives.v[2] += (p->a_hydro[2] + p->gravity.old_a[2]) * dt;
-
-//#if !defined(EOS_ISOTHERMAL_GAS)
-// if (p->conserved.mass > 0.) {
-// const float u = p->conserved.energy + p->du_dt * dt;
-// p->primitives.P =
-// hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass;
-// }
-//#endif
+ /* drift the primitive variables based on the old fluxes */
+ p->primitives.rho += p->conserved.flux.mass * dt / p->geometry.volume;
+
+ p->primitives.v[0] += p->conserved.flux.momentum[0] * dt / p->conserved.mass;
+ p->primitives.v[1] += p->conserved.flux.momentum[1] * dt / p->conserved.mass;
+ p->primitives.v[2] += p->conserved.flux.momentum[2] * dt / p->conserved.mass;
+
+#if !defined(EOS_ISOTHERMAL_GAS)
+ if (p->conserved.mass > 0.) {
+ const float u = p->conserved.energy + p->conserved.flux.energy * dt;
+ p->primitives.P =
+ hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass;
+ }
+#endif
#ifdef SWIFT_DEBUG_CHECKS
if (p->h <= 0.) {
@@ -580,12 +575,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
/* set the variables that are used to drift the primitive variables */
- if (p->force.dt > 0.) {
- p->du_dt = p->conserved.flux.energy / p->force.dt;
- } else {
- p->du_dt = 0.0f;
- }
-
hydro_velocities_end_force(p);
}
@@ -605,16 +594,16 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
float a_grav[3];
/* Update conserved variables. */
- p->conserved.mass += p->conserved.flux.mass;
- p->conserved.momentum[0] += p->conserved.flux.momentum[0];
- p->conserved.momentum[1] += p->conserved.flux.momentum[1];
- p->conserved.momentum[2] += p->conserved.flux.momentum[2];
+ p->conserved.mass += p->conserved.flux.mass * dt;
+ p->conserved.momentum[0] += p->conserved.flux.momentum[0] * dt;
+ p->conserved.momentum[1] += p->conserved.flux.momentum[1] * dt;
+ p->conserved.momentum[2] += p->conserved.flux.momentum[2] * dt;
#if defined(EOS_ISOTHERMAL_GAS)
/* We use the EoS equation in a sneaky way here just to get the constant u */
p->conserved.energy =
p->conserved.mass * gas_internal_energy_from_entropy(0.f, 0.f);
#else
- p->conserved.energy += p->conserved.flux.energy;
+ p->conserved.energy += p->conserved.flux.energy * dt;
#endif
gizmo_check_physical_quantity("mass", p->conserved.mass);
@@ -669,15 +658,6 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
}
- /* reset fluxes */
- /* we can only do this here, since we need to keep the fluxes for inactive
- particles */
- 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;
-
hydro_velocities_set(p, xp);
#ifdef GIZMO_LLOYD_ITERATION
diff --git a/src/hydro/Gizmo/hydro_debug.h b/src/hydro/Gizmo/hydro_debug.h
index 17e7f8a08570e355a701f8e165ee8af745fa34ab..f11104f8a6847300e47ffd6618acc4498bb5e648 100644
--- a/src/hydro/Gizmo/hydro_debug.h
+++ b/src/hydro/Gizmo/hydro_debug.h
@@ -48,9 +48,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"timestepvars={"
"vmax=%.3e},"
"density={"
- "div_v=%.3e, "
"wcount_dh=%.3e, "
- "curl_v=[%.3e,%.3e,%.3e], "
"wcount=%.3e}\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
p->a_hydro[1], p->a_hydro[2], p->h, p->time_bin, p->primitives.v[0],
@@ -75,7 +73,5 @@ __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.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);
+ p->timestepvars.vmax, p->density.wcount_dh, p->density.wcount);
}
diff --git a/src/hydro/Gizmo/hydro_flux_limiters.h b/src/hydro/Gizmo/hydro_flux_limiters.h
deleted file mode 100644
index dc91cf2808e02d903ff97efddc20c164db9c954e..0000000000000000000000000000000000000000
--- a/src/hydro/Gizmo/hydro_flux_limiters.h
+++ /dev/null
@@ -1,81 +0,0 @@
-
-/*******************************************************************************
- * 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 896128bd45d7964c1f4c8d63564f6fced38db770..b01b4d90b845520c2bd20b4d33279b5d533999d2 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -93,7 +93,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
*/
__attribute__((always_inline)) INLINE static void hydro_gradients_predict(
struct part* pi, struct part* pj, float hi, float hj, float* dx, float r,
- float* xij_i, float* Wi, float* Wj, float mindt) {
+ float* xij_i, float* Wi, float* Wj) {
float dWi[5], dWj[5];
float xij_j[3];
@@ -142,63 +142,6 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
- /* time */
- if (Wi[0] > 0.0f) {
- dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
- Wi[2] * pi->primitives.gradients.rho[1] +
- Wi[3] * pi->primitives.gradients.rho[2] +
- Wi[0] * (pi->primitives.gradients.v[0][0] +
- pi->primitives.gradients.v[1][1] +
- pi->primitives.gradients.v[2][2]));
- dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
- Wi[2] * pi->primitives.gradients.v[0][1] +
- Wi[3] * pi->primitives.gradients.v[0][2] +
- pi->primitives.gradients.P[0] / Wi[0]);
- dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
- Wi[2] * pi->primitives.gradients.v[1][1] +
- Wi[3] * pi->primitives.gradients.v[1][2] +
- pi->primitives.gradients.P[1] / Wi[0]);
- dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
- Wi[2] * pi->primitives.gradients.v[2][1] +
- Wi[3] * pi->primitives.gradients.v[2][2] +
- pi->primitives.gradients.P[2] / Wi[0]);
- dWi[4] -= 0.5 * mindt *
- (Wi[1] * pi->primitives.gradients.P[0] +
- Wi[2] * pi->primitives.gradients.P[1] +
- Wi[3] * pi->primitives.gradients.P[2] +
- hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
- pi->primitives.gradients.v[1][1] +
- pi->primitives.gradients.v[2][2]));
- }
-
- if (Wj[0] > 0.0f) {
- dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
- Wj[2] * pj->primitives.gradients.rho[1] +
- Wj[3] * pj->primitives.gradients.rho[2] +
- Wj[0] * (pj->primitives.gradients.v[0][0] +
- pj->primitives.gradients.v[1][1] +
- pj->primitives.gradients.v[2][2]));
- dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
- Wj[2] * pj->primitives.gradients.v[0][1] +
- Wj[3] * pj->primitives.gradients.v[0][2] +
- pj->primitives.gradients.P[0] / Wj[0]);
- dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
- Wj[2] * pj->primitives.gradients.v[1][1] +
- Wj[3] * pj->primitives.gradients.v[1][2] +
- pj->primitives.gradients.P[1] / Wj[0]);
- dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
- Wj[2] * pj->primitives.gradients.v[2][1] +
- Wj[3] * pj->primitives.gradients.v[2][2] +
- pj->primitives.gradients.P[2] / Wj[0]);
- dWj[4] -= 0.5 * mindt *
- (Wj[1] * pj->primitives.gradients.P[0] +
- Wj[2] * pj->primitives.gradients.P[1] +
- Wj[3] * pj->primitives.gradients.P[2] +
- hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
- pj->primitives.gradients.v[1][1] +
- pj->primitives.gradients.v[2][2]));
- }
-
Wi[0] += dWi[0];
Wi[1] += dWi[1];
Wi[2] += dWi[2];
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index 3dd66eecb408213fef5c03fae56c33b7d4de34cb..c2ed5dfd6fd4f137f44b35887b299d3c3ba8337a 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -20,7 +20,6 @@
******************************************************************************/
#include "adiabatic_index.h"
-#include "hydro_flux_limiters.h"
#include "hydro_gradients.h"
#include "riemann.h"
@@ -150,53 +149,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
__attribute__((always_inline)) INLINE static void runner_iact_gradient(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
- float hi_inv, hi_inv_dim, xi, wi, wi_dx;
- float hj_inv, hj_inv_dim, xj, wj, wj_dx;
- float Bi[3][3], Bj[3][3];
- float Vi, Vj;
- float A, Anorm;
- int k, l;
- float r;
-
- r = sqrtf(r2);
-
- hi_inv = 1.0 / hi;
- hi_inv_dim = pow_dimension(hi_inv);
- xi = r * hi_inv;
- kernel_deval(xi, &wi, &wi_dx);
-
- /* Compute kernel of pj. */
- hj_inv = 1.0 / hj;
- hj_inv_dim = pow_dimension(hj_inv);
- xj = r * hj_inv;
- kernel_deval(xj, &wj, &wj_dx);
-
- for (k = 0; k < 3; k++) {
- for (l = 0; l < 3; l++) {
- Bi[k][l] = pi->geometry.matrix_E[k][l];
- Bj[k][l] = pj->geometry.matrix_E[k][l];
- }
- }
- Vi = pi->geometry.volume;
- Vj = pj->geometry.volume;
-
- /* Compute area */
- /* eqn. (7) */
- Anorm = 0.0f;
- for (k = 0; k < 3; k++) {
- /* we add a minus sign since dx is pi->x - pj->x */
- A = -Vi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) * wi *
- hi_inv_dim -
- Vj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) * wj *
- hj_inv_dim;
- Anorm += A * A;
- }
-
- Anorm = sqrtf(Anorm);
-
- pi->geometry.Atot += Anorm;
- pj->geometry.Atot += Anorm;
-
hydro_gradients_collect(r2, dx, hi, hj, pi, pj);
}
@@ -217,52 +169,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_gradient(
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
- float hi_inv, hi_inv_dim, xi, wi, wi_dx;
- float hj_inv, hj_inv_dim, xj, wj, wj_dx;
- float Bi[3][3], Bj[3][3];
- float Vi, Vj;
- float A, Anorm;
- int k, l;
- float r;
-
- r = sqrtf(r2);
-
- hi_inv = 1.0 / hi;
- hi_inv_dim = pow_dimension(hi_inv);
- xi = r * hi_inv;
- kernel_deval(xi, &wi, &wi_dx);
-
- /* Compute kernel of pj. */
- hj_inv = 1.0 / hj;
- hj_inv_dim = pow_dimension(hj_inv);
- xj = r * hj_inv;
- kernel_deval(xj, &wj, &wj_dx);
-
- for (k = 0; k < 3; k++) {
- for (l = 0; l < 3; l++) {
- Bi[k][l] = pi->geometry.matrix_E[k][l];
- Bj[k][l] = pj->geometry.matrix_E[k][l];
- }
- }
- Vi = pi->geometry.volume;
- Vj = pj->geometry.volume;
-
- /* Compute area */
- /* eqn. (7) */
- Anorm = 0.0f;
- for (k = 0; k < 3; k++) {
- /* we add a minus sign since dx is pi->x - pj->x */
- A = -Vi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) * wi *
- hi_inv_dim -
- Vj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) * wj *
- hj_inv_dim;
- Anorm += A * A;
- }
-
- Anorm = sqrtf(Anorm);
-
- pi->geometry.Atot += Anorm;
-
hydro_gradients_nonsym_collect(r2, dx, hi, hj, pi, pj);
}
@@ -271,9 +177,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
*
* Since the only difference between the symmetric and non-symmetric version
* of the flux calculation is in the update of the conserved variables at the
- * very end (which is not done for particle j if mode is 0 and particle j is
- * active), both runner_iact_force and runner_iact_nonsym_force call this
- * method, with an appropriate mode.
+ * very end (which is not done for particle j if mode is 0), both
+ * runner_iact_force and runner_iact_nonsym_force call this method, with an
+ * appropriate mode.
*
* This method calculates the surface area of the interface between particle i
* and particle j, as well as the interface position and velocity. These are
@@ -310,7 +216,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float vmax, dvdotdx;
float vi[3], vj[3], vij[3];
float Wi[5], Wj[5];
- float dti, dtj, mindt;
float n_unit[3];
/* Initialize local variables */
@@ -319,8 +224,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Bi[k][l] = pi->geometry.matrix_E[k][l];
Bj[k][l] = pj->geometry.matrix_E[k][l];
}
- vi[k] = pi->force.v_full[k]; /* particle velocities */
- vj[k] = pj->force.v_full[k];
+ vi[k] = pi->v[k]; /* particle velocities */
+ vj[k] = pj->v[k];
}
Vi = pi->geometry.volume;
Vj = pj->geometry.volume;
@@ -335,9 +240,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Wj[3] = pj->primitives.v[2];
Wj[4] = pj->primitives.P;
- dti = pi->force.dt;
- dtj = pj->force.dt;
-
/* calculate the maximal signal velocity */
if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
vmax =
@@ -358,10 +260,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
pj->timestepvars.vmax = max(pj->timestepvars.vmax, vmax);
}
- /* The flux will be exchanged using the smallest time step of the two
- * particles */
- mindt = min(dti, dtj);
-
/* Compute kernel of pi. */
hi_inv = 1.0 / hi;
hi_inv_dim = pow_dimension(hi_inv);
@@ -411,9 +309,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
for (k = 0; k < 3; k++) {
/* we add a minus sign since dx is pi->x - pj->x */
A[k] = -Xi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) *
- wj * hj_inv_dim -
+ wi * hi_inv_dim -
Xj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) *
- wi * hi_inv_dim;
+ wj * hj_inv_dim;
Anorm += A[k] * A[k];
}
} else {
@@ -483,7 +381,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
Wj[2] -= vij[1];
Wj[3] -= vij[2];
- hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj, mindt);
+ hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj);
/* we don't need to rotate, we can use the unit vector in the Riemann problem
* itself (see GIZMO) */
@@ -491,69 +389,63 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float totflux[5];
riemann_solve_for_flux(Wi, Wj, n_unit, vij, totflux);
- hydro_flux_limiters_apply(totflux, pi, pj);
+ /* Multiply with the interface surface area */
+ totflux[0] *= Anorm;
+ totflux[1] *= Anorm;
+ totflux[2] *= Anorm;
+ totflux[3] *= Anorm;
+ totflux[4] *= Anorm;
/* Store mass flux */
- float mflux = Anorm * totflux[0];
+ float mflux = totflux[0];
pi->gravity.mflux[0] += mflux * dx[0];
pi->gravity.mflux[1] += mflux * dx[1];
pi->gravity.mflux[2] += mflux * dx[2];
/* Update conserved variables */
/* eqn. (16) */
- pi->conserved.flux.mass -= mindt * Anorm * totflux[0];
- pi->conserved.flux.momentum[0] -= mindt * Anorm * totflux[1];
- pi->conserved.flux.momentum[1] -= mindt * Anorm * totflux[2];
- pi->conserved.flux.momentum[2] -= mindt * Anorm * totflux[3];
- pi->conserved.flux.energy -= mindt * Anorm * totflux[4];
+ pi->conserved.flux.mass -= totflux[0];
+ pi->conserved.flux.momentum[0] -= totflux[1];
+ pi->conserved.flux.momentum[1] -= totflux[2];
+ pi->conserved.flux.momentum[2] -= totflux[3];
+ pi->conserved.flux.energy -= totflux[4];
#ifndef GIZMO_TOTAL_ENERGY
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.flux.energy += mindt * Anorm * totflux[1] * pi->primitives.v[0];
- pi->conserved.flux.energy += mindt * Anorm * totflux[2] * pi->primitives.v[1];
- pi->conserved.flux.energy += mindt * Anorm * totflux[3] * pi->primitives.v[2];
- pi->conserved.flux.energy -= mindt * Anorm * totflux[0] * ekin;
+ pi->conserved.flux.energy += totflux[1] * pi->primitives.v[0];
+ pi->conserved.flux.energy += totflux[2] * pi->primitives.v[1];
+ pi->conserved.flux.energy += totflux[3] * pi->primitives.v[2];
+ pi->conserved.flux.energy -= totflux[0] * ekin;
#endif
- /* 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->force.active == 0) {
+ /* Note that this used to be much more complicated in early implementations of
+ * the GIZMO scheme, as we wanted manifest conservation of conserved variables
+ * and had to do symmetric flux exchanges. Now we don't care about manifest
+ * conservation anymore and just assume the current fluxes are representative
+ * for the flux over the entire time step. */
+ if (mode == 1) {
/* Store mass flux */
- mflux = Anorm * totflux[0];
+ mflux = totflux[0];
pj->gravity.mflux[0] -= mflux * dx[0];
pj->gravity.mflux[1] -= mflux * dx[1];
pj->gravity.mflux[2] -= mflux * dx[2];
- pj->conserved.flux.mass += mindt * Anorm * totflux[0];
- pj->conserved.flux.momentum[0] += mindt * Anorm * totflux[1];
- pj->conserved.flux.momentum[1] += mindt * Anorm * totflux[2];
- pj->conserved.flux.momentum[2] += mindt * Anorm * totflux[3];
- pj->conserved.flux.energy += mindt * Anorm * totflux[4];
+ pj->conserved.flux.mass += totflux[0];
+ pj->conserved.flux.momentum[0] += totflux[1];
+ pj->conserved.flux.momentum[1] += totflux[2];
+ pj->conserved.flux.momentum[2] += totflux[3];
+ pj->conserved.flux.energy += totflux[4];
#ifndef GIZMO_TOTAL_ENERGY
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.flux.energy -=
- mindt * Anorm * totflux[1] * pj->primitives.v[0];
- pj->conserved.flux.energy -=
- mindt * Anorm * totflux[2] * pj->primitives.v[1];
- pj->conserved.flux.energy -=
- mindt * Anorm * totflux[3] * pj->primitives.v[2];
- pj->conserved.flux.energy += mindt * Anorm * totflux[0] * ekin;
+ pj->conserved.flux.energy -= totflux[1] * pj->primitives.v[0];
+ pj->conserved.flux.energy -= totflux[2] * pj->primitives.v[1];
+ pj->conserved.flux.energy -= totflux[3] * pj->primitives.v[2];
+ pj->conserved.flux.energy += totflux[0] * ekin;
#endif
}
}
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index 458ea6151e4b00b9ccf85ad6eff4ccf58277ff7d..9d6acf4155f02dc8caaa92be5366389ebb5ecaef 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -19,7 +19,6 @@
#include "adiabatic_index.h"
#include "hydro.h"
-#include "hydro_flux_limiters.h"
#include "hydro_gradients.h"
#include "hydro_slope_limiters.h"
#include "io_properties.h"
@@ -173,10 +172,6 @@ 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 47f722c5a2dcce2f3ce603ade3029821d6686067..f13da395bb8dfe0ccc04e7fb823c47c1f2db4b02 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/Gizmo/hydro_part.h
@@ -59,18 +59,15 @@ struct part {
/* Particle smoothing length. */
float h;
- /* Old internal energy flux */
- float du_dt;
-
/* The primitive hydrodynamical variables. */
struct {
- /* Fluid velocity. */
- float v[3];
-
/* Density. */
float rho;
+ /* Fluid velocity. */
+ float v[3];
+
/* Pressure. */
float P;
@@ -110,12 +107,12 @@ struct part {
/* The conserved hydrodynamical variables. */
struct {
- /* Fluid momentum. */
- float momentum[3];
-
/* Fluid mass */
float mass;
+ /* Fluid momentum. */
+ float momentum[3];
+
/* Fluid thermal energy (not per unit mass!). */
float energy;
@@ -145,9 +142,6 @@ struct part {
gradients */
float matrix_E[3][3];
- /* Total surface area of the particle. */
- float Atot;
-
/* Centroid of the "cell". */
float centroid[3];
@@ -167,15 +161,9 @@ struct part {
/* Quantities used during the volume (=density) loop. */
struct {
- /* Particle velocity divergence. */
- float div_v;
-
/* Derivative of particle number density. */
float wcount_dh;
- /* Particle velocity curl. */
- float curl_v[3];
-
/* Particle number density. */
float wcount;
@@ -190,15 +178,6 @@ struct part {
/* Needed to drift the primitive variables. */
float h_dt;
- /* Physical time step of the particle. */
- float dt;
-
- /* Flag keeping track of whether this is an active or inactive particle. */
- char active;
-
- /* Actual velocity of the particle. */
- float v_full[3];
-
} force;
/* Specific stuff for the gravity-hydro coupling. */
diff --git a/src/hydro/Gizmo/hydro_velocities.h b/src/hydro/Gizmo/hydro_velocities.h
index 08ba1f972b2f7a7b8a01ac4750c50a36f69784d0..eef5d408e66ea6b4bb1be0cf53452a81ea0b4e89 100644
--- a/src/hydro/Gizmo/hydro_velocities.h
+++ b/src/hydro/Gizmo/hydro_velocities.h
@@ -49,18 +49,11 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_init(
* velocities during the force loop.
*
* @param p The particle to act upon.
- * @param xp The extended particel data to act upon.
+ * @param xp The extended particle data to act upon.
*/
__attribute__((always_inline)) INLINE static void
hydro_velocities_prepare_force(struct part* restrict p,
- const struct xpart* restrict xp) {
-
-#ifndef GIZMO_FIX_PARTICLES
- p->force.v_full[0] = xp->v_full[0];
- p->force.v_full[1] = xp->v_full[1];
- p->force.v_full[2] = xp->v_full[2];
-#endif
-}
+ const struct xpart* restrict xp) {}
/**
* @brief Set the variables that will be used to update the smoothing length