Skip to content
GitLab
Commit 5a540b01 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'master' into update_cooling_grackle

parents af67ece6 55e98e2f
Hide whitespace changes
Inline Side-by-side
src/Makefile.am
View file @ 5a540b01
......@@ -94,7 +94,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 \
......
src/hydro/Gizmo/hydro.h
View file @ 5a540b01
/*******************************************************************************
* 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
......
src/hydro/Gizmo/hydro_debug.h
View file @ 5a540b01
......@@ -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);
}
src/hydro/Gizmo/hydro_flux_limiters.h deleted 100644 → 0
View file @ af67ece6
/*******************************************************************************
* 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
src/hydro/Gizmo/hydro_gradients.h
View file @ 5a540b01
......@@ -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];
......
src/hydro/Gizmo/hydro_iact.h
View file @ 5a540b01
......@@ -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];