diff --git a/examples/SineWavePotential_1D/makeIC.py b/examples/SineWavePotential_1D/makeIC.py
index 321af0714219dfa0c7cbb3d80845881dcbb8416d..afbf1bc0fa47a27677cb9c5645d439432bd9fd9a 100644
--- a/examples/SineWavePotential_1D/makeIC.py
+++ b/examples/SineWavePotential_1D/makeIC.py
@@ -31,12 +31,12 @@ cs2 = 2.*uconst/3.
A = 10.
fileName = "sineWavePotential.hdf5"
-numPart = 100
+numPart = 1000
boxSize = 1.
coords = np.zeros((numPart, 3))
v = np.zeros((numPart, 3))
-m = np.zeros(numPart) + 1.
+m = np.zeros(numPart) + 1000. / numPart
h = np.zeros(numPart) + 2./numPart
u = np.zeros(numPart) + uconst
ids = np.arange(numPart, dtype = 'L')
diff --git a/examples/SineWavePotential_1D/plotSolution.py b/examples/SineWavePotential_1D/plotSolution.py
index 65e981e4648fe0fe5d1da6cf3e753fb8a34f0fb4..3bb889aaabd3cdac0274afb09647d0e3aebb02cc 100644
--- a/examples/SineWavePotential_1D/plotSolution.py
+++ b/examples/SineWavePotential_1D/plotSolution.py
@@ -23,8 +23,13 @@
import numpy as np
import h5py
import sys
+import matplotlib
+matplotlib.use("Agg")
import pylab as pl
+pl.rcParams["figure.figsize"] = (12, 10)
+pl.rcParams["text.usetex"] = True
+
# these should be the same as in makeIC.py
uconst = 20.2615290634
cs2 = 2.*uconst/3.
@@ -39,15 +44,20 @@ fileName = sys.argv[1]
file = h5py.File(fileName, 'r')
coords = np.array(file["/PartType0/Coordinates"])
rho = np.array(file["/PartType0/Density"])
+P = np.array(file["/PartType0/Pressure"])
u = np.array(file["/PartType0/InternalEnergy"])
-agrav = np.array(file["/PartType0/GravAcceleration"])
m = np.array(file["/PartType0/Masses"])
+vs = np.array(file["/PartType0/Velocities"])
ids = np.array(file["/PartType0/ParticleIDs"])
x = np.linspace(0., 1., 1000)
rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
-P = cs2*rho
+a = A * np.sin(2. * np.pi * x)
+
+apart = A * np.sin(2. * np.pi * coords[:,0])
+tkin = -0.5 * np.dot(apart, coords[:,0])
+print tkin, 0.5 * np.dot(m, vs[:,0]**2)
ids_reverse = np.argsort(ids)
@@ -65,13 +75,38 @@ for i in range(len(P)):
corr = 1.
idxp1 = ids_reverse[ip1]
idxm1 = ids_reverse[im1]
- gradP[i] = (P[idxp1]-P[idxm1])/(coords[idxp1,0]-coords[idxm1,0])
+ gradP[i] = (P[idxp1] - P[idxm1])/(coords[idxp1,0] - coords[idxm1,0])
+
+fig, ax = pl.subplots(2, 2, sharex = True)
+
+ax[0][0].plot(coords[:,0], rho, ".", label = "gizmo-mfm")
+ax[0][0].plot(x, rho_x, "-", label = "stable solution")
+ax[0][0].set_ylabel("$\\rho{}$ (code units)")
+ax[0][0].legend(loc = "best")
+
+ax[0][1].plot(x, a, "-", label = "$\\nabla{}\\Phi{}$ external")
+ax[0][1].plot(coords[:,0], gradP/rho, ".",
+ label = "$\\nabla{}P/\\rho{}$ gizmo-mfm")
+ax[0][1].set_ylabel("force (code units)")
+ax[0][1].legend(loc = "best")
+
+ax[1][0].axhline(y = uconst, label = "isothermal $u$", color = "k",
+ linestyle = "--")
+ax[1][0].plot(coords[:,0], u, ".", label = "gizmo-mfm")
+ax[1][0].set_ylabel("$u$ (code units)")
+ax[1][0].set_xlabel("$x$ (code units)")
+ax[1][0].legend(loc = "best")
+
+#ax[1][1].plot(coords[:,0], m, "y.")
+#ax[1][1].set_ylabel("$m_i$ (code units)")
+#ax[1][1].set_xlabel("$x$ (code units)")
-fig, ax = pl.subplots(2, 2)
+ax[1][1].axhline(y = 0., label = "target", color = "k",
+ linestyle = "--")
+ax[1][1].plot(coords[:,0], vs[:,0], ".", label = "gizmo-mfm")
+ax[1][1].set_ylabel("$v_x$ (code units)")
+ax[1][1].set_xlabel("$x$ (code units)")
+ax[1][1].legend(loc = "best")
-ax[0][0].plot(coords[:,0], rho, "r.", markersize = 0.5)
-ax[0][0].plot(x, rho_x, "g-")
-ax[0][1].plot(coords[:,0], gradP/rho, "b.")
-ax[1][0].plot(coords[:,0], agrav[:,0], "g.", markersize = 0.5)
-ax[1][1].plot(coords[:,0], m, "y.")
+pl.tight_layout()
pl.savefig("{fileName}.png".format(fileName = fileName[:-5]))
diff --git a/examples/SineWavePotential_1D/sineWavePotential.yml b/examples/SineWavePotential_1D/sineWavePotential.yml
index 9662841032a12d48870f12de8c1bcfcd579a6b42..e6285785099f10902ea60b21334a0ad26c0593de 100644
--- a/examples/SineWavePotential_1D/sineWavePotential.yml
+++ b/examples/SineWavePotential_1D/sineWavePotential.yml
@@ -10,7 +10,7 @@ InternalUnitSystem:
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 10. # The end time of the simulation (in internal units).
- dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
+ dt_min: 1e-8 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the conserved quantities statistics
@@ -36,3 +36,6 @@ InitialConditions:
SineWavePotential:
amplitude: 10.
timestep_limit: 1.
+
+EoS:
+ isothermal_internal_energy: 20.2615290634
diff --git a/src/Makefile.am b/src/Makefile.am
index fbc1d9fdc025dca2a498208ddce98e716fa2fd03..0f61fb108d8d8acbd420c994266f4803a3f69d3e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -106,7 +106,6 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/GizmoMFM/hydro_slope_limiters_face.h \
hydro/GizmoMFM/hydro_slope_limiters.h \
hydro/GizmoMFM/hydro_unphysical.h \
- hydro/GizmoMFM/hydro_velocities.h \
hydro/Shadowswift/hydro_debug.h \
hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
hydro/Shadowswift/hydro_iact.h \
diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h
index 1ab142740b641bdc9a0dff5a02b19479bae8257e..7c38896a4bfb0c44a79c509954e11128244e350e 100644
--- a/src/hydro/GizmoMFM/hydro.h
+++ b/src/hydro/GizmoMFM/hydro.h
@@ -28,14 +28,11 @@
#include "hydro_properties.h"
#include "hydro_space.h"
#include "hydro_unphysical.h"
-#include "hydro_velocities.h"
#include "minmax.h"
#include "riemann.h"
#include <float.h>
-//#define GIZMO_LLOYD_ITERATION
-
/**
* @brief Computes the hydro time-step of a given particle
*
@@ -51,20 +48,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const float CFL_condition = hydro_properties->CFL_condition;
-#ifdef GIZMO_LLOYD_ITERATION
- return CFL_condition;
-#endif
-
- /* v_full is the actual velocity of the particle, primitives.v is its
+ /* v_full is the actual velocity of the particle, 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];
- vrel[2] = p->primitives.v[2] - xp->v_full[2];
+ vrel[0] = p->v[0] - xp->v_full[0];
+ vrel[1] = p->v[1] - xp->v_full[1];
+ vrel[2] = p->v[2] - xp->v_full[2];
float vmax =
sqrtf(vrel[0] * vrel[0] + vrel[1] * vrel[1] + vrel[2] * vrel[2]) +
- sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+ sqrtf(hydro_gamma * p->P / p->rho);
vmax = max(vmax, p->timestepvars.vmax);
// MATTHIEU: Bert is this correct? Do we need more cosmology terms here?
@@ -72,7 +65,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
cosmo->a * powf(p->geometry.volume / hydro_dimension_unit_sphere,
hydro_dimension_inv);
float dt = FLT_MAX;
- if (vmax > 0.) {
+ if (vmax > 0.0f) {
dt = psize / vmax;
}
return CFL_condition * dt;
@@ -92,7 +85,7 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
struct part* p, float dt) {
#ifdef SWIFT_DEBUG_CHECKS
- if (dt == 0.) {
+ if (dt == 0.0f) {
error("Zero time step assigned to particle!");
}
@@ -122,14 +115,10 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
const float mass = p->conserved.mass;
- p->primitives.v[0] = p->v[0];
- p->primitives.v[1] = p->v[1];
- p->primitives.v[2] = p->v[2];
-
/* we can already initialize the momentum */
- p->conserved.momentum[0] = mass * p->primitives.v[0];
- p->conserved.momentum[1] = mass * p->primitives.v[1];
- p->conserved.momentum[2] = mass * p->primitives.v[2];
+ p->conserved.momentum[0] = mass * p->v[0];
+ p->conserved.momentum[1] = mass * p->v[1];
+ p->conserved.momentum[2] = mass * p->v[2];
/* and the thermal energy */
/* remember that we store the total thermal energy, not the specific thermal
@@ -137,39 +126,22 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
#if defined(EOS_ISOTHERMAL_GAS)
/* this overwrites the internal energy from the initial condition file
* Note that we call the EoS function just to get the constant u here. */
- p->conserved.energy = mass * gas_internal_energy_from_entropy(0.f, 0.f);
+ p->conserved.energy = mass * gas_internal_energy_from_entropy(0.0f, 0.0f);
#else
p->conserved.energy *= mass;
#endif
#ifdef GIZMO_TOTAL_ENERGY
/* add the total kinetic energy */
- p->conserved.energy += 0.5f * (p->conserved.momentum[0] * p->primitives.v[0] +
- p->conserved.momentum[1] * p->primitives.v[1] +
- p->conserved.momentum[2] * p->primitives.v[2]);
-#endif
-
-#ifdef GIZMO_LLOYD_ITERATION
- /* overwrite all variables to make sure they have safe values */
- p->primitives.rho = 1.;
- p->primitives.v[0] = 0.;
- p->primitives.v[1] = 0.;
- p->primitives.v[2] = 0.;
- p->primitives.P = 1.;
-
- p->conserved.mass = 1.;
- p->conserved.momentum[0] = 0.;
- p->conserved.momentum[1] = 0.;
- p->conserved.momentum[2] = 0.;
- p->conserved.energy = 1.;
-
- p->v[0] = 0.;
- p->v[1] = 0.;
- p->v[2] = 0.;
+ p->conserved.energy += 0.5f * (p->conserved.momentum[0] * p->v[0] +
+ p->conserved.momentum[1] * p->v[1] +
+ p->conserved.momentum[2] * p->v[2]);
#endif
/* initialize the particle velocity based on the primitive fluid velocity */
- hydro_velocities_init(p, xp);
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
/* ignore accelerations present in the initial condition */
p->a_hydro[0] = 0.0f;
@@ -180,7 +152,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
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
particle */
- p->density.wcorr = 1.0f;
+ p->geometry.wcorr = 1.0f;
}
/**
@@ -250,7 +222,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Final operation on the geometry. */
/* we multiply with the smoothing kernel normalization ih3 and calculate the
* volume */
- const float volume = 1.f / (ihdim * (p->geometry.volume + kernel_root));
+ const float volume_inv = ihdim * (p->geometry.volume + kernel_root);
+ const float volume = 1.0f / volume_inv;
p->geometry.volume = volume;
/* we multiply with the smoothing kernel normalization */
@@ -268,9 +241,10 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->geometry.centroid[1] *= kernel_norm;
p->geometry.centroid[2] *= kernel_norm;
- p->geometry.centroid[0] /= p->density.wcount;
- p->geometry.centroid[1] /= p->density.wcount;
- p->geometry.centroid[2] /= p->density.wcount;
+ const float wcount_inv = 1.0f / p->density.wcount;
+ p->geometry.centroid[0] *= wcount_inv;
+ p->geometry.centroid[1] *= wcount_inv;
+ p->geometry.centroid[2] *= wcount_inv;
/* Check the condition number to see if we have a stable geometry. */
float condition_number_E = 0.0f;
@@ -296,7 +270,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
hydro_dimension_inv * sqrtf(condition_number_E * condition_number_Einv);
if (condition_number > const_gizmo_max_condition_number &&
- p->density.wcorr > const_gizmo_min_wcorr) {
+ p->geometry.wcorr > const_gizmo_min_wcorr) {
#ifdef GIZMO_PATHOLOGICAL_ERROR
error("Condition number larger than %g (%g)!",
const_gizmo_max_condition_number, condition_number);
@@ -306,21 +280,19 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
condition_number, const_gizmo_max_condition_number, p->id);
#endif
/* add a correction to the number of neighbours for this particle */
- p->density.wcorr *= const_gizmo_w_correction_factor;
+ p->geometry.wcorr *= const_gizmo_w_correction_factor;
}
- hydro_gradients_init(p);
-
/* compute primitive variables */
/* eqns (3)-(5) */
const float m = p->conserved.mass;
#ifdef SWIFT_DEBUG_CHECKS
- if (m < 0.) {
+ if (m < 0.0f) {
error("Mass is negative!");
}
- if (volume == 0.) {
+ if (volume == 0.0f) {
error("Volume is 0!");
}
#endif
@@ -330,55 +302,45 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
momentum[0] = p->conserved.momentum[0];
momentum[1] = p->conserved.momentum[1];
momentum[2] = p->conserved.momentum[2];
- p->primitives.rho = m / volume;
- if (m == 0.) {
- p->primitives.v[0] = 0.;
- p->primitives.v[1] = 0.;
- p->primitives.v[2] = 0.;
+ p->rho = m * volume_inv;
+ if (m == 0.0f) {
+ p->v[0] = 0.0f;
+ p->v[1] = 0.0f;
+ p->v[2] = 0.0f;
} else {
- p->primitives.v[0] = momentum[0] / m;
- p->primitives.v[1] = momentum[1] / m;
- p->primitives.v[2] = momentum[2] / m;
+ const float m_inv = 1.0f / m;
+ p->v[0] = momentum[0] * m_inv;
+ p->v[1] = momentum[1] * m_inv;
+ p->v[2] = momentum[2] * m_inv;
}
#ifdef EOS_ISOTHERMAL_GAS
/* although the pressure is not formally used anywhere if an isothermal eos
has been selected, we still make sure it is set to the correct value */
- p->primitives.P = gas_pressure_from_internal_energy(p->primitives.rho, 0.);
+ p->P = gas_pressure_from_internal_energy(p->rho, 0.0f);
#else
float energy = p->conserved.energy;
#ifdef GIZMO_TOTAL_ENERGY
/* subtract the kinetic energy; we want the thermal energy */
- energy -= 0.5f * (momentum[0] * p->primitives.v[0] +
- momentum[1] * p->primitives.v[1] +
- momentum[2] * p->primitives.v[2]);
+ energy -= 0.5f * (momentum[0] * p->v[0] + momentum[1] * p->v[1] +
+ momentum[2] * p->v[2]);
#endif
/* energy contains the total thermal energy, we want the specific energy.
this is why we divide by the volume, and not by the density */
- p->primitives.P = hydro_gamma_minus_one * energy / volume;
+ p->P = hydro_gamma_minus_one * energy * volume_inv;
#endif
/* sanity checks */
- gizmo_check_physical_quantities("density", "pressure", p->primitives.rho,
- p->primitives.v[0], p->primitives.v[1],
- p->primitives.v[2], p->primitives.P);
-
-#ifdef GIZMO_LLOYD_ITERATION
- /* overwrite primitive variables to make sure they still have safe values */
- p->primitives.rho = 1.;
- p->primitives.v[0] = 0.;
- p->primitives.v[1] = 0.;
- p->primitives.v[2] = 0.;
- p->primitives.P = 1.;
-#endif
+ gizmo_check_physical_quantities("density", "pressure", p->rho, p->v[0],
+ p->v[1], p->v[2], p->P);
/* Add a correction factor to wcount (to force a neighbour number increase if
the geometry matrix is close to singular) */
- p->density.wcount *= p->density.wcorr;
- p->density.wcount_dh *= p->density.wcorr;
+ p->density.wcount *= p->geometry.wcorr;
+ p->density.wcount_dh *= p->geometry.wcorr;
}
/**
@@ -399,7 +361,7 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
/* Re-set problematic values */
p->density.wcount = kernel_root * h_inv_dim;
- p->density.wcount_dh = 0.f;
+ p->density.wcount_dh = 0.0f;
p->geometry.volume = 1.0f;
p->geometry.matrix_E[0][0] = 1.0f;
p->geometry.matrix_E[0][1] = 0.0f;
@@ -437,12 +399,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
const struct cosmology* cosmo) {
/* Initialize time step criterion variables */
- p->timestepvars.vmax = 0.;
+ p->timestepvars.vmax = 0.0f;
- // MATTHIEU: Bert is this correct? Do we need cosmology terms here?
+ hydro_gradients_init(p);
- /* Set the actual velocity of the particle */
- hydro_velocities_prepare_force(p, xp);
+ // MATTHIEU: Bert is this correct? Do we need cosmology terms here?
}
/**
@@ -494,15 +455,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const struct cosmology* cosmo) {
/* Initialise values that are used in the force loop */
- p->gravity.mflux[0] = 0.0f;
- 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;
+ p->flux.momentum[0] = 0.0f;
+ p->flux.momentum[1] = 0.0f;
+ p->flux.momentum[2] = 0.0f;
+ p->flux.energy = 0.0f;
}
/**
@@ -559,10 +515,6 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part* p, struct xpart* xp, float dt_drift, float dt_therm) {
-#ifdef GIZMO_LLOYD_ITERATION
- return;
-#endif
-
const float h_inv = 1.0f / p->h;
/* Predict smoothing length */
@@ -575,39 +527,45 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Limit the smoothing length correction (and make sure it is always
positive). */
- if (h_corr < 2.0f && h_corr > 0.) {
+ if (h_corr < 2.0f && h_corr > 0.0f) {
p->h *= h_corr;
}
/* drift the primitive variables based on the old fluxes */
- if (p->geometry.volume > 0.) {
- p->primitives.rho += p->conserved.flux.mass * dt_drift / p->geometry.volume;
- }
+ if (p->conserved.mass > 0.0f) {
+ const float m_inv = 1.0f / p->conserved.mass;
- if (p->conserved.mass > 0.) {
- p->primitives.v[0] +=
- p->conserved.flux.momentum[0] * dt_drift / p->conserved.mass;
- p->primitives.v[1] +=
- p->conserved.flux.momentum[1] * dt_drift / p->conserved.mass;
- p->primitives.v[2] +=
- p->conserved.flux.momentum[2] * dt_drift / p->conserved.mass;
+ p->v[0] += p->flux.momentum[0] * dt_drift * m_inv;
+ p->v[1] += p->flux.momentum[1] * dt_drift * m_inv;
+ p->v[2] += p->flux.momentum[2] * dt_drift * m_inv;
#if !defined(EOS_ISOTHERMAL_GAS)
- const float u = p->conserved.energy + p->conserved.flux.energy * dt_therm;
- p->primitives.P =
- hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+ const float Etot = p->conserved.energy + p->flux.energy * dt_therm;
+ const float v2 =
+ (p->v[0] * p->v[0] + p->v[1] * p->v[1] + p->v[2] * p->v[2]);
+ const float u = (Etot * m_inv - 0.5f * v2);
+#else
+ const float u = (p->conserved.energy + p->flux.energy * dt_therm) * m_inv;
+#endif
+ p->P = hydro_gamma_minus_one * u * p->rho;
#endif
}
+ /* we use a sneaky way to get the gravitational contribtuion to the
+ velocity update */
+ p->v[0] += p->v[0] - xp->v_full[0];
+ p->v[1] += p->v[1] - xp->v_full[1];
+ p->v[2] += p->v[2] - xp->v_full[2];
+
#ifdef SWIFT_DEBUG_CHECKS
- if (p->h <= 0.) {
+ if (p->h <= 0.0f) {
error("Zero or negative smoothing length (%g)!", p->h);
}
#endif
- gizmo_check_physical_quantities("density", "pressure", p->primitives.rho,
- p->primitives.v[0], p->primitives.v[1],
- p->primitives.v[2], p->primitives.P);
+ gizmo_check_physical_quantities("density", "pressure", p->rho, p->v[0],
+ p->v[1], p->v[2], p->P);
}
/**
@@ -629,7 +587,9 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
/* set the variables that are used to drift the primitive variables */
// MATTHIEU: Bert is this correct? Do we need cosmology terms here?
- hydro_velocities_end_force(p);
+
+ /* Add normalization to h_dt. */
+ p->force.h_dt *= p->h * hydro_dimension_inv;
}
/**
@@ -648,17 +608,16 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
float a_grav[3];
- /* Update conserved variables. */
- 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;
+ /* Update conserved variables (note: the mass does not change). */
+ p->conserved.momentum[0] += p->flux.momentum[0] * dt;
+ p->conserved.momentum[1] += p->flux.momentum[1] * dt;
+ p->conserved.momentum[2] += p->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);
+ p->conserved.mass * gas_internal_energy_from_entropy(0.0f, 0.0f);
#else
- p->conserved.energy += p->conserved.flux.energy * dt;
+ p->conserved.energy += p->flux.energy * dt;
#endif
/* Apply the minimal energy limit */
@@ -666,7 +625,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
hydro_props->minimal_internal_energy * cosmo->a_factor_internal_energy;
if (p->conserved.energy < min_energy * p->conserved.mass) {
p->conserved.energy = min_energy * p->conserved.mass;
- p->conserved.flux.energy = 0.f;
+ p->flux.energy = 0.0f;
}
gizmo_check_physical_quantities(
@@ -676,17 +635,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
#ifdef SWIFT_DEBUG_CHECKS
/* Note that this check will only have effect if no GIZMO_UNPHYSICAL option
was selected. */
- if (p->conserved.mass < 0.) {
- error(
- "Negative mass after conserved variables update (mass: %g, dmass: %g)!",
- p->conserved.mass, p->conserved.flux.mass);
- }
-
- if (p->conserved.energy < 0.) {
+ if (p->conserved.energy < 0.0f) {
error(
"Negative energy after conserved variables update (energy: %g, "
"denergy: %g)!",
- p->conserved.energy, p->conserved.flux.energy);
+ p->conserved.energy, p->flux.energy);
}
#endif
@@ -699,44 +652,40 @@ __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];
- /* Make sure the gpart knows the mass has changed. */
- p->gpart->mass = p->conserved.mass;
-
/* Kick the momentum for half a time step */
/* Note that this also affects the particle movement, as the velocity for
the particles is set after this. */
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];
-
- p->conserved.energy += dt * (p->gravity.mflux[0] * a_grav[0] +
- p->gravity.mflux[1] * a_grav[1] +
- p->gravity.mflux[2] * a_grav[2]);
}
- hydro_velocities_set(p, xp);
+ /* Set the velocities: */
+ /* We first set the particle velocity */
+ if (p->conserved.mass > 0.0f && p->rho > 0.0f) {
-#ifdef GIZMO_LLOYD_ITERATION
- /* reset conserved variables to safe values */
- p->conserved.mass = 1.;
- p->conserved.momentum[0] = 0.;
- p->conserved.momentum[1] = 0.;
- p->conserved.momentum[2] = 0.;
- p->conserved.energy = 1.;
-
- /* set the particle velocities to the Lloyd velocities */
- /* note that centroid is the relative position of the centroid w.r.t. the
- particle position (position - centroid) */
- xp->v_full[0] = -p->geometry.centroid[0] / p->force.dt;
- xp->v_full[1] = -p->geometry.centroid[1] / p->force.dt;
- xp->v_full[2] = -p->geometry.centroid[2] / p->force.dt;
- p->v[0] = xp->v_full[0];
- p->v[1] = xp->v_full[1];
- p->v[2] = xp->v_full[2];
-#endif
+ const float inverse_mass = 1.0f / p->conserved.mass;
+
+ /* Normal case: set particle velocity to fluid velocity. */
+ xp->v_full[0] = p->conserved.momentum[0] * inverse_mass;
+ xp->v_full[1] = p->conserved.momentum[1] * inverse_mass;
+ xp->v_full[2] = p->conserved.momentum[2] * inverse_mass;
+
+ } else {
+ /* Vacuum particles have no fluid velocity. */
+ xp->v_full[0] = 0.0f;
+ xp->v_full[1] = 0.0f;
+ xp->v_full[2] = 0.0f;
+ }
+
+ if (p->gpart) {
+ p->gpart->v_full[0] = xp->v_full[0];
+ p->gpart->v_full[1] = xp->v_full[1];
+ p->gpart->v_full[2] = xp->v_full[2];
+ }
/* reset wcorr */
- p->density.wcorr = 1.0f;
+ p->geometry.wcorr = 1.0f;
}
/**
@@ -747,11 +696,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_internal_energy(const struct part* restrict p) {
- if (p->primitives.rho > 0.)
- return gas_internal_energy_from_pressure(p->primitives.rho,
- p->primitives.P);
- else
- return 0.;
+ if (p->rho > 0.0f) {
+ return gas_internal_energy_from_pressure(p->rho, p->P);
+ } else {
+ return 0.0f;
+ }
}
/**
@@ -776,10 +725,10 @@ hydro_get_physical_internal_energy(const struct part* restrict p,
__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
const struct part* restrict p) {
- if (p->primitives.rho > 0.) {
- return gas_entropy_from_pressure(p->primitives.rho, p->primitives.P);
+ if (p->rho > 0.0f) {
+ return gas_entropy_from_pressure(p->rho, p->P);
} else {
- return 0.;
+ return 0.0f;
}
}
@@ -805,10 +754,11 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_soundspeed(const struct part* restrict p) {
- if (p->primitives.rho > 0.)
- return gas_soundspeed_from_pressure(p->primitives.rho, p->primitives.P);
- else
- return 0.;
+ if (p->rho > 0.0f) {
+ return gas_soundspeed_from_pressure(p->rho, p->P);
+ } else {
+ return 0.0f;
+ }
}
/**
@@ -832,7 +782,7 @@ hydro_get_physical_soundspeed(const struct part* restrict p,
__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
const struct part* restrict p) {
- return p->primitives.P;
+ return p->P;
}
/**
@@ -844,7 +794,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
const struct part* restrict p, const struct cosmology* cosmo) {
- return cosmo->a_factor_pressure * p->primitives.P;
+ return cosmo->a_factor_pressure * p->P;
}
/**
@@ -883,17 +833,15 @@ __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
const struct part* restrict p, const struct xpart* xp, float dt_kick_hydro,
float dt_kick_grav, float v[3]) {
- if (p->conserved.mass > 0.) {
- v[0] = p->primitives.v[0] +
- p->conserved.flux.momentum[0] * dt_kick_hydro / p->conserved.mass;
- v[1] = p->primitives.v[1] +
- p->conserved.flux.momentum[1] * dt_kick_hydro / p->conserved.mass;
- v[2] = p->primitives.v[2] +
- p->conserved.flux.momentum[2] * dt_kick_hydro / p->conserved.mass;
+ if (p->conserved.mass > 0.0f) {
+ const float inverse_mass = 1.0f / p->conserved.mass;
+ v[0] = p->v[0] + p->flux.momentum[0] * dt_kick_hydro * inverse_mass;
+ v[1] = p->v[1] + p->flux.momentum[1] * dt_kick_hydro * inverse_mass;
+ v[2] = p->v[2] + p->flux.momentum[2] * dt_kick_hydro * inverse_mass;
} else {
- v[0] = p->primitives.v[0];
- v[1] = p->primitives.v[1];
- v[2] = p->primitives.v[2];
+ v[0] = p->v[0];
+ v[1] = p->v[1];
+ v[2] = p->v[2];
}
// MATTHIEU: Bert is this correct?
@@ -910,7 +858,7 @@ __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
const struct part* restrict p) {
- return p->primitives.rho;
+ return p->rho;
}
/**
@@ -922,7 +870,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
__attribute__((always_inline)) INLINE static float hydro_get_physical_density(
const struct part* restrict p, const struct cosmology* cosmo) {
- return cosmo->a3_inv * p->primitives.rho;
+ return cosmo->a3_inv * p->rho;
}
/**
@@ -943,12 +891,12 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
p->conserved.energy = u * p->conserved.mass;
#ifdef GIZMO_TOTAL_ENERGY
/* add the kinetic energy */
- p->conserved.energy += 0.5f * p->conserved.mass *
- (p->conserved.momentum[0] * p->primitives.v[0] +
- p->conserved.momentum[1] * p->primitives.v[1] +
- p->conserved.momentum[2] * p->primitives.v[2]);
+ p->conserved.energy +=
+ 0.5f * p->conserved.mass *
+ (p->conserved.momentum[0] * p->v[0] + p->conserved.momentum[1] * p->v[1] +
+ p->conserved.momentum[2] * p->v[2]);
#endif
- p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u;
+ p->P = hydro_gamma_minus_one * p->rho * u;
}
/**
@@ -963,16 +911,16 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
__attribute__((always_inline)) INLINE static void hydro_set_entropy(
struct part* restrict p, float S) {
- p->conserved.energy = S * pow_gamma_minus_one(p->primitives.rho) *
+ p->conserved.energy = S * pow_gamma_minus_one(p->rho) *
hydro_one_over_gamma_minus_one * p->conserved.mass;
#ifdef GIZMO_TOTAL_ENERGY
/* add the kinetic energy */
- p->conserved.energy += 0.5f * p->conserved.mass *
- (p->conserved.momentum[0] * p->primitives.v[0] +
- p->conserved.momentum[1] * p->primitives.v[1] +
- p->conserved.momentum[2] * p->primitives.v[2]);
+ p->conserved.energy +=
+ 0.5f * p->conserved.mass *
+ (p->conserved.momentum[0] * p->v[0] + p->conserved.momentum[1] * p->v[1] +
+ p->conserved.momentum[2] * p->v[2]);
#endif
- p->primitives.P = S * pow_gamma(p->primitives.rho);
+ p->P = S * pow_gamma(p->rho);
}
/**
@@ -992,12 +940,12 @@ hydro_set_init_internal_energy(struct part* p, float u_init) {
p->conserved.energy = u_init * p->conserved.mass;
#ifdef GIZMO_TOTAL_ENERGY
/* add the kinetic energy */
- p->conserved.energy += 0.5f * p->conserved.mass *
- (p->conserved.momentum[0] * p->primitives.v[0] +
- p->conserved.momentum[1] * p->primitives.v[1] +
- p->conserved.momentum[2] * p->primitives.v[2]);
+ p->conserved.energy +=
+ 0.5f * p->conserved.mass *
+ (p->conserved.momentum[0] * p->v[0] + p->conserved.momentum[1] * p->v[1] +
+ p->conserved.momentum[2] * p->v[2]);
#endif
- p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u_init;
+ p->P = hydro_gamma_minus_one * p->rho * u_init;
}
#endif /* SWIFT_GIZMO_MFM_HYDRO_H */
diff --git a/src/hydro/GizmoMFM/hydro_debug.h b/src/hydro/GizmoMFM/hydro_debug.h
index 6603bc216b986b40513383120587d3caec1adc87..e8b0914bd3cf6a99210399c6fc654e526319009f 100644
--- a/src/hydro/GizmoMFM/hydro_debug.h
+++ b/src/hydro/GizmoMFM/hydro_debug.h
@@ -27,8 +27,6 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"a=[%.3e,%.3e,%.3e], "
"h=%.3e, "
"time_bin=%d, "
- "primitives={"
- "v=[%.3e,%.3e,%.3e], "
"rho=%.3e, "
"P=%.3e, "
"gradients={"
@@ -39,7 +37,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"rho=[%.3e,%.3e], "
"v=[[%.3e,%.3e],[%.3e,%.3e],[%.3e,%.3e]], "
"P=[%.3e,%.3e], "
- "maxr=%.3e}}, "
+ "maxr=%.3e}, "
"conserved={"
"momentum=[%.3e,%.3e,%.3e], "
"mass=%.3e, "
@@ -53,24 +51,18 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"wcount_dh=%.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],
- p->primitives.v[1], p->primitives.v[2], p->primitives.rho,
- p->primitives.P, p->primitives.gradients.rho[0],
- p->primitives.gradients.rho[1], p->primitives.gradients.rho[2],
- p->primitives.gradients.v[0][0], p->primitives.gradients.v[0][1],
- p->primitives.gradients.v[0][2], p->primitives.gradients.v[1][0],
- p->primitives.gradients.v[1][1], p->primitives.gradients.v[1][2],
- p->primitives.gradients.v[2][0], p->primitives.gradients.v[2][1],
- p->primitives.gradients.v[2][2], p->primitives.gradients.P[0],
- p->primitives.gradients.P[1], p->primitives.gradients.P[2],
- p->primitives.limiter.rho[0], p->primitives.limiter.rho[1],
- p->primitives.limiter.v[0][0], p->primitives.limiter.v[0][1],
- p->primitives.limiter.v[1][0], p->primitives.limiter.v[1][1],
- p->primitives.limiter.v[2][0], p->primitives.limiter.v[2][1],
- p->primitives.limiter.P[0], p->primitives.limiter.P[1],
- p->primitives.limiter.maxr, p->conserved.momentum[0],
- p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.mass,
- p->conserved.energy, p->geometry.volume, p->geometry.matrix_E[0][0],
+ p->a_hydro[1], p->a_hydro[2], p->h, p->time_bin, p->rho, p->P,
+ p->gradients.rho[0], p->gradients.rho[1], p->gradients.rho[2],
+ p->gradients.v[0][0], p->gradients.v[0][1], p->gradients.v[0][2],
+ p->gradients.v[1][0], p->gradients.v[1][1], p->gradients.v[1][2],
+ p->gradients.v[2][0], p->gradients.v[2][1], p->gradients.v[2][2],
+ p->gradients.P[0], p->gradients.P[1], p->gradients.P[2],
+ p->limiter.rho[0], p->limiter.rho[1], p->limiter.v[0][0],
+ p->limiter.v[0][1], p->limiter.v[1][0], p->limiter.v[1][1],
+ p->limiter.v[2][0], p->limiter.v[2][1], p->limiter.P[0], p->limiter.P[1],
+ p->limiter.maxr, p->conserved.momentum[0], p->conserved.momentum[1],
+ p->conserved.momentum[2], p->conserved.mass, p->conserved.energy,
+ p->geometry.volume, p->geometry.matrix_E[0][0],
p->geometry.matrix_E[0][1], p->geometry.matrix_E[0][2],
p->geometry.matrix_E[1][0], p->geometry.matrix_E[1][1],
p->geometry.matrix_E[1][2], p->geometry.matrix_E[2][0],
diff --git a/src/hydro/GizmoMFM/hydro_gradients.h b/src/hydro/GizmoMFM/hydro_gradients.h
index 85617f05b6a5350e7b1b60de53f2c999b6ba7453..6f751d970287ca7ba137c1138ca6f3bf00e6c4cb 100644
--- a/src/hydro/GizmoMFM/hydro_gradients.h
+++ b/src/hydro/GizmoMFM/hydro_gradients.h
@@ -101,38 +101,28 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
const float xij_j[3] = {xij_i[0] + dx[0], xij_i[1] + dx[1], xij_i[2] + dx[2]};
float dWi[5];
- dWi[0] = pi->primitives.gradients.rho[0] * xij_i[0] +
- pi->primitives.gradients.rho[1] * xij_i[1] +
- pi->primitives.gradients.rho[2] * xij_i[2];
- dWi[1] = pi->primitives.gradients.v[0][0] * xij_i[0] +
- pi->primitives.gradients.v[0][1] * xij_i[1] +
- pi->primitives.gradients.v[0][2] * xij_i[2];
- dWi[2] = pi->primitives.gradients.v[1][0] * xij_i[0] +
- pi->primitives.gradients.v[1][1] * xij_i[1] +
- pi->primitives.gradients.v[1][2] * xij_i[2];
- dWi[3] = pi->primitives.gradients.v[2][0] * xij_i[0] +
- pi->primitives.gradients.v[2][1] * xij_i[1] +
- pi->primitives.gradients.v[2][2] * xij_i[2];
- dWi[4] = pi->primitives.gradients.P[0] * xij_i[0] +
- pi->primitives.gradients.P[1] * xij_i[1] +
- pi->primitives.gradients.P[2] * xij_i[2];
+ dWi[0] = pi->gradients.rho[0] * xij_i[0] + pi->gradients.rho[1] * xij_i[1] +
+ pi->gradients.rho[2] * xij_i[2];
+ dWi[1] = pi->gradients.v[0][0] * xij_i[0] + pi->gradients.v[0][1] * xij_i[1] +
+ pi->gradients.v[0][2] * xij_i[2];
+ dWi[2] = pi->gradients.v[1][0] * xij_i[0] + pi->gradients.v[1][1] * xij_i[1] +
+ pi->gradients.v[1][2] * xij_i[2];
+ dWi[3] = pi->gradients.v[2][0] * xij_i[0] + pi->gradients.v[2][1] * xij_i[1] +
+ pi->gradients.v[2][2] * xij_i[2];
+ dWi[4] = pi->gradients.P[0] * xij_i[0] + pi->gradients.P[1] * xij_i[1] +
+ pi->gradients.P[2] * xij_i[2];
float dWj[5];
- dWj[0] = pj->primitives.gradients.rho[0] * xij_j[0] +
- pj->primitives.gradients.rho[1] * xij_j[1] +
- pj->primitives.gradients.rho[2] * xij_j[2];
- dWj[1] = pj->primitives.gradients.v[0][0] * xij_j[0] +
- pj->primitives.gradients.v[0][1] * xij_j[1] +
- pj->primitives.gradients.v[0][2] * xij_j[2];
- dWj[2] = pj->primitives.gradients.v[1][0] * xij_j[0] +
- pj->primitives.gradients.v[1][1] * xij_j[1] +
- pj->primitives.gradients.v[1][2] * xij_j[2];
- dWj[3] = pj->primitives.gradients.v[2][0] * xij_j[0] +
- pj->primitives.gradients.v[2][1] * xij_j[1] +
- pj->primitives.gradients.v[2][2] * xij_j[2];
- dWj[4] = pj->primitives.gradients.P[0] * xij_j[0] +
- pj->primitives.gradients.P[1] * xij_j[1] +
- pj->primitives.gradients.P[2] * xij_j[2];
+ dWj[0] = pj->gradients.rho[0] * xij_j[0] + pj->gradients.rho[1] * xij_j[1] +
+ pj->gradients.rho[2] * xij_j[2];
+ dWj[1] = pj->gradients.v[0][0] * xij_j[0] + pj->gradients.v[0][1] * xij_j[1] +
+ pj->gradients.v[0][2] * xij_j[2];
+ dWj[2] = pj->gradients.v[1][0] * xij_j[0] + pj->gradients.v[1][1] * xij_j[1] +
+ pj->gradients.v[1][2] * xij_j[2];
+ dWj[3] = pj->gradients.v[2][0] * xij_j[0] + pj->gradients.v[2][1] * xij_j[1] +
+ pj->gradients.v[2][2] * xij_j[2];
+ dWj[4] = pj->gradients.P[0] * xij_j[0] + pj->gradients.P[1] * xij_j[1] +
+ pj->gradients.P[2] * xij_j[2];
/* Apply the slope limiter at this interface */
hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
diff --git a/src/hydro/GizmoMFM/hydro_gradients_gizmo.h b/src/hydro/GizmoMFM/hydro_gradients_gizmo.h
index 1c3b68bb28375259628e09f16730710fbbd80149..90c8096a85b3418c8ee4c5382d6f087c14f8907e 100644
--- a/src/hydro/GizmoMFM/hydro_gradients_gizmo.h
+++ b/src/hydro/GizmoMFM/hydro_gradients_gizmo.h
@@ -28,25 +28,25 @@
__attribute__((always_inline)) INLINE static void hydro_gradients_init(
struct part *p) {
- p->primitives.gradients.rho[0] = 0.0f;
- p->primitives.gradients.rho[1] = 0.0f;
- p->primitives.gradients.rho[2] = 0.0f;
+ p->gradients.rho[0] = 0.0f;
+ p->gradients.rho[1] = 0.0f;
+ p->gradients.rho[2] = 0.0f;
- p->primitives.gradients.v[0][0] = 0.0f;
- p->primitives.gradients.v[0][1] = 0.0f;
- p->primitives.gradients.v[0][2] = 0.0f;
+ p->gradients.v[0][0] = 0.0f;
+ p->gradients.v[0][1] = 0.0f;
+ p->gradients.v[0][2] = 0.0f;
- p->primitives.gradients.v[1][0] = 0.0f;
- p->primitives.gradients.v[1][1] = 0.0f;
- p->primitives.gradients.v[1][2] = 0.0f;
+ p->gradients.v[1][0] = 0.0f;
+ p->gradients.v[1][1] = 0.0f;
+ p->gradients.v[1][2] = 0.0f;
- p->primitives.gradients.v[2][0] = 0.0f;
- p->primitives.gradients.v[2][1] = 0.0f;
- p->primitives.gradients.v[2][2] = 0.0f;
+ p->gradients.v[2][0] = 0.0f;
+ p->gradients.v[2][1] = 0.0f;
+ p->gradients.v[2][2] = 0.0f;
- p->primitives.gradients.P[0] = 0.0f;
- p->primitives.gradients.P[1] = 0.0f;
- p->primitives.gradients.P[2] = 0.0f;
+ p->gradients.P[0] = 0.0f;
+ p->gradients.P[1] = 0.0f;
+ p->gradients.P[2] = 0.0f;
hydro_slope_limit_cell_init(p);
}
@@ -65,7 +65,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj) {
- const float r_inv = 1.f / sqrtf(r2);
+ const float r_inv = 1.0f / sqrtf(r2);
const float r = r2 * r_inv;
float wi, wj, wi_dx, wj_dx;
@@ -80,211 +80,100 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
Bj[k][l] = pj->geometry.matrix_E[k][l];
}
}
- Wi[0] = pi->primitives.rho;
- Wi[1] = pi->primitives.v[0];
- Wi[2] = pi->primitives.v[1];
- Wi[3] = pi->primitives.v[2];
- Wi[4] = pi->primitives.P;
- Wj[0] = pj->primitives.rho;
- Wj[1] = pj->primitives.v[0];
- Wj[2] = pj->primitives.v[1];
- Wj[3] = pj->primitives.v[2];
- Wj[4] = pj->primitives.P;
+ Wi[0] = pi->rho;
+ Wi[1] = pi->v[0];
+ Wi[2] = pi->v[1];
+ Wi[3] = pi->v[2];
+ Wi[4] = pi->P;
+ Wj[0] = pj->rho;
+ Wj[1] = pj->v[0];
+ Wj[2] = pj->v[1];
+ Wj[3] = pj->v[2];
+ Wj[4] = pj->P;
/* Compute kernel of pi. */
- const float hi_inv = 1.f / hi;
+ const float hi_inv = 1.0f / hi;
const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
- if (pi->density.wcorr > const_gizmo_min_wcorr) {
- /* Compute gradients for pi */
- /* there is a sign difference w.r.t. eqn. (6) because of the inverse
- * definition of dx */
- pi->primitives.gradients.rho[0] +=
- (Wi[0] - Wj[0]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.rho[1] +=
- (Wi[0] - Wj[0]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.rho[2] +=
- (Wi[0] - Wj[0]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
- pi->primitives.gradients.v[0][0] +=
- (Wi[1] - Wj[1]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.v[0][1] +=
- (Wi[1] - Wj[1]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.v[0][2] +=
- (Wi[1] - Wj[1]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->primitives.gradients.v[1][0] +=
- (Wi[2] - Wj[2]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.v[1][1] +=
- (Wi[2] - Wj[2]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.v[1][2] +=
- (Wi[2] - Wj[2]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->primitives.gradients.v[2][0] +=
- (Wi[3] - Wj[3]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.v[2][1] +=
- (Wi[3] - Wj[3]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.v[2][2] +=
- (Wi[3] - Wj[3]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
- pi->primitives.gradients.P[0] +=
- (Wi[4] - Wj[4]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.P[1] +=
- (Wi[4] - Wj[4]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.P[2] +=
- (Wi[4] - Wj[4]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ const float dW[5] = {Wi[0] - Wj[0], Wi[1] - Wj[1], Wi[2] - Wj[2],
+ Wi[3] - Wj[3], Wi[4] - Wj[4]};
+ float wiBidx[3];
+ if (pi->geometry.wcorr > const_gizmo_min_wcorr) {
+ wiBidx[0] = wi * (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ wiBidx[1] = wi * (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ wiBidx[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 */
-
- pi->primitives.gradients.rho[0] -=
- wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
- pi->primitives.gradients.rho[1] -=
- wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
- pi->primitives.gradients.rho[2] -=
- wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-
- pi->primitives.gradients.v[0][0] -=
- wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
- pi->primitives.gradients.v[0][1] -=
- wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
- pi->primitives.gradients.v[0][2] -=
- wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-
- pi->primitives.gradients.v[1][0] -=
- wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pi->primitives.gradients.v[1][1] -=
- wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pi->primitives.gradients.v[1][2] -=
- wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
-
- pi->primitives.gradients.v[2][0] -=
- wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
- pi->primitives.gradients.v[2][1] -=
- wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
- pi->primitives.gradients.v[2][2] -=
- wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-
- pi->primitives.gradients.P[0] -=
- wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
- pi->primitives.gradients.P[1] -=
- wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
- pi->primitives.gradients.P[2] -=
- wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ const float norm = -wi_dx * r_inv;
+ wiBidx[0] = norm * dx[0];
+ wiBidx[1] = norm * dx[1];
+ wiBidx[2] = norm * dx[2];
}
+ /* Compute gradients for pi */
+ /* there is a sign difference w.r.t. eqn. (6) because of the inverse
+ * definition of dx */
+ pi->gradients.rho[0] += dW[0] * wiBidx[0];
+ pi->gradients.rho[1] += dW[0] * wiBidx[1];
+ pi->gradients.rho[2] += dW[0] * wiBidx[2];
+
+ pi->gradients.v[0][0] += dW[1] * wiBidx[0];
+ pi->gradients.v[0][1] += dW[1] * wiBidx[1];
+ pi->gradients.v[0][2] += dW[1] * wiBidx[2];
+ pi->gradients.v[1][0] += dW[2] * wiBidx[0];
+ pi->gradients.v[1][1] += dW[2] * wiBidx[1];
+ pi->gradients.v[1][2] += dW[2] * wiBidx[2];
+ pi->gradients.v[2][0] += dW[3] * wiBidx[0];
+ pi->gradients.v[2][1] += dW[3] * wiBidx[1];
+ pi->gradients.v[2][2] += dW[3] * wiBidx[2];
+
+ pi->gradients.P[0] += dW[4] * wiBidx[0];
+ pi->gradients.P[1] += dW[4] * wiBidx[1];
+ pi->gradients.P[2] += dW[4] * wiBidx[2];
+
hydro_slope_limit_cell_collect(pi, pj, r);
/* Compute kernel of pj. */
- const float hj_inv = 1.f / hj;
+ const float hj_inv = 1.0f / hj;
const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
- if (pj->density.wcorr > const_gizmo_min_wcorr) {
- /* Compute gradients for pj */
- /* there is no sign difference w.r.t. eqn. (6) because dx is now what we
- * want
- * it to be */
- pj->primitives.gradients.rho[0] +=
- (Wi[0] - Wj[0]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.rho[1] +=
- (Wi[0] - Wj[0]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.rho[2] +=
- (Wi[0] - Wj[0]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-
- pj->primitives.gradients.v[0][0] +=
- (Wi[1] - Wj[1]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.v[0][1] +=
- (Wi[1] - Wj[1]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.v[0][2] +=
- (Wi[1] - Wj[1]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
- pj->primitives.gradients.v[1][0] +=
- (Wi[2] - Wj[2]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.v[1][1] +=
- (Wi[2] - Wj[2]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.v[1][2] +=
- (Wi[2] - Wj[2]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
- pj->primitives.gradients.v[2][0] +=
- (Wi[3] - Wj[3]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.v[2][1] +=
- (Wi[3] - Wj[3]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.v[2][2] +=
- (Wi[3] - Wj[3]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
-
- pj->primitives.gradients.P[0] +=
- (Wi[4] - Wj[4]) * wj *
- (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
- pj->primitives.gradients.P[1] +=
- (Wi[4] - Wj[4]) * wj *
- (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
- pj->primitives.gradients.P[2] +=
- (Wi[4] - Wj[4]) * wj *
- (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+ float wjBjdx[3];
+ if (pj->geometry.wcorr > const_gizmo_min_wcorr) {
+
+ wjBjdx[0] = wj * (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ wjBjdx[1] = wj * (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ wjBjdx[2] = wj * (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
} else {
- /* SPH gradients */
-
- pj->primitives.gradients.rho[0] -=
- wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
- pj->primitives.gradients.rho[1] -=
- wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
- pj->primitives.gradients.rho[2] -=
- wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-
- pj->primitives.gradients.v[0][0] -=
- wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
- pj->primitives.gradients.v[0][1] -=
- wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
- pj->primitives.gradients.v[0][2] -=
- wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-
- pj->primitives.gradients.v[1][0] -=
- wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pj->primitives.gradients.v[1][1] -=
- wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pj->primitives.gradients.v[1][2] -=
- wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pj->primitives.gradients.v[2][0] -=
- wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
- pj->primitives.gradients.v[2][1] -=
- wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
- pj->primitives.gradients.v[2][2] -=
- wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-
- pj->primitives.gradients.P[0] -=
- wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
- pj->primitives.gradients.P[1] -=
- wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
- pj->primitives.gradients.P[2] -=
- wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ const float norm = -wj_dx * r_inv;
+ wjBjdx[0] = norm * dx[0];
+ wjBjdx[1] = norm * dx[1];
+ wjBjdx[2] = norm * dx[2];
}
+ /* Compute gradients for pj */
+ /* there is no sign difference w.r.t. eqn. (6) because dx is now what we
+ * want it to be */
+ pj->gradients.rho[0] += dW[0] * wjBjdx[0];
+ pj->gradients.rho[1] += dW[0] * wjBjdx[1];
+ pj->gradients.rho[2] += dW[0] * wjBjdx[2];
+
+ pj->gradients.v[0][0] += dW[1] * wjBjdx[0];
+ pj->gradients.v[0][1] += dW[1] * wjBjdx[1];
+ pj->gradients.v[0][2] += dW[1] * wjBjdx[2];
+ pj->gradients.v[1][0] += dW[2] * wjBjdx[0];
+ pj->gradients.v[1][1] += dW[2] * wjBjdx[1];
+ pj->gradients.v[1][2] += dW[2] * wjBjdx[2];
+ pj->gradients.v[2][0] += dW[3] * wjBjdx[0];
+ pj->gradients.v[2][1] += dW[3] * wjBjdx[1];
+ pj->gradients.v[2][2] += dW[3] * wjBjdx[2];
+
+ pj->gradients.P[0] += dW[4] * wjBjdx[0];
+ pj->gradients.P[1] += dW[4] * wjBjdx[1];
+ pj->gradients.P[2] += dW[4] * wjBjdx[2];
+
hydro_slope_limit_cell_collect(pj, pi, r);
}
@@ -303,7 +192,7 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
struct part *restrict pi,
struct part *restrict pj) {
- const float r_inv = 1.f / sqrtf(r2);
+ const float r_inv = 1.0f / sqrtf(r2);
const float r = r2 * r_inv;
float Bi[3][3];
@@ -315,113 +204,59 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
Bi[k][l] = pi->geometry.matrix_E[k][l];
}
}
- Wi[0] = pi->primitives.rho;
- Wi[1] = pi->primitives.v[0];
- Wi[2] = pi->primitives.v[1];
- Wi[3] = pi->primitives.v[2];
- Wi[4] = pi->primitives.P;
- Wj[0] = pj->primitives.rho;
- Wj[1] = pj->primitives.v[0];
- Wj[2] = pj->primitives.v[1];
- Wj[3] = pj->primitives.v[2];
- Wj[4] = pj->primitives.P;
+ Wi[0] = pi->rho;
+ Wi[1] = pi->v[0];
+ Wi[2] = pi->v[1];
+ Wi[3] = pi->v[2];
+ Wi[4] = pi->P;
+ Wj[0] = pj->rho;
+ Wj[1] = pj->v[0];
+ Wj[2] = pj->v[1];
+ Wj[3] = pj->v[2];
+ Wj[4] = pj->P;
/* Compute kernel of pi. */
float wi, wi_dx;
- const float hi_inv = 1.f / hi;
+ const float hi_inv = 1.0f / hi;
const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
- if (pi->density.wcorr > const_gizmo_min_wcorr) {
- /* Compute gradients for pi */
- /* there is a sign difference w.r.t. eqn. (6) because of the inverse
- * definition of dx */
- pi->primitives.gradients.rho[0] +=
- (Wi[0] - Wj[0]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.rho[1] +=
- (Wi[0] - Wj[0]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.rho[2] +=
- (Wi[0] - Wj[0]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
- pi->primitives.gradients.v[0][0] +=
- (Wi[1] - Wj[1]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.v[0][1] +=
- (Wi[1] - Wj[1]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.v[0][2] +=
- (Wi[1] - Wj[1]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->primitives.gradients.v[1][0] +=
- (Wi[2] - Wj[2]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.v[1][1] +=
- (Wi[2] - Wj[2]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.v[1][2] +=
- (Wi[2] - Wj[2]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
- pi->primitives.gradients.v[2][0] +=
- (Wi[3] - Wj[3]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.v[2][1] +=
- (Wi[3] - Wj[3]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.v[2][2] +=
- (Wi[3] - Wj[3]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
-
- pi->primitives.gradients.P[0] +=
- (Wi[4] - Wj[4]) * wi *
- (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
- pi->primitives.gradients.P[1] +=
- (Wi[4] - Wj[4]) * wi *
- (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
- pi->primitives.gradients.P[2] +=
- (Wi[4] - Wj[4]) * wi *
- (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ const float dW[5] = {Wi[0] - Wj[0], Wi[1] - Wj[1], Wi[2] - Wj[2],
+ Wi[3] - Wj[3], Wi[4] - Wj[4]};
+ float wiBidx[3];
+ if (pi->geometry.wcorr > const_gizmo_min_wcorr) {
+ wiBidx[0] = wi * (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ wiBidx[1] = wi * (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ wiBidx[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 */
-
- pi->primitives.gradients.rho[0] -=
- wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
- pi->primitives.gradients.rho[1] -=
- wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
- pi->primitives.gradients.rho[2] -=
- wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-
- pi->primitives.gradients.v[0][0] -=
- wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
- pi->primitives.gradients.v[0][1] -=
- wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
- pi->primitives.gradients.v[0][2] -=
- wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
- pi->primitives.gradients.v[1][0] -=
- wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pi->primitives.gradients.v[1][1] -=
- wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pi->primitives.gradients.v[1][2] -=
- wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
-
- pi->primitives.gradients.v[2][0] -=
- wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
- pi->primitives.gradients.v[2][1] -=
- wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
- pi->primitives.gradients.v[2][2] -=
- wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-
- pi->primitives.gradients.P[0] -=
- wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
- pi->primitives.gradients.P[1] -=
- wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
- pi->primitives.gradients.P[2] -=
- wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ const float norm = -wi_dx * r_inv;
+ wiBidx[0] = norm * dx[0];
+ wiBidx[1] = norm * dx[1];
+ wiBidx[2] = norm * dx[2];
}
+ /* Compute gradients for pi */
+ /* there is a sign difference w.r.t. eqn. (6) because of the inverse
+ * definition of dx */
+ pi->gradients.rho[0] += dW[0] * wiBidx[0];
+ pi->gradients.rho[1] += dW[0] * wiBidx[1];
+ pi->gradients.rho[2] += dW[0] * wiBidx[2];
+
+ pi->gradients.v[0][0] += dW[1] * wiBidx[0];
+ pi->gradients.v[0][1] += dW[1] * wiBidx[1];
+ pi->gradients.v[0][2] += dW[1] * wiBidx[2];
+ pi->gradients.v[1][0] += dW[2] * wiBidx[0];
+ pi->gradients.v[1][1] += dW[2] * wiBidx[1];
+ pi->gradients.v[1][2] += dW[2] * wiBidx[2];
+ pi->gradients.v[2][0] += dW[3] * wiBidx[0];
+ pi->gradients.v[2][1] += dW[3] * wiBidx[1];
+ pi->gradients.v[2][2] += dW[3] * wiBidx[2];
+
+ pi->gradients.P[0] += dW[4] * wiBidx[0];
+ pi->gradients.P[1] += dW[4] * wiBidx[1];
+ pi->gradients.P[2] += dW[4] * wiBidx[2];
+
hydro_slope_limit_cell_collect(pi, pj, r);
}
@@ -438,50 +273,33 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
const float h = p->h;
const float h_inv = 1.0f / h;
const float ihdim = pow_dimension(h_inv);
- const float ihdimp1 = pow_dimension_plus_one(h_inv);
-
- if (p->density.wcorr > const_gizmo_min_wcorr) {
- p->primitives.gradients.rho[0] *= ihdim;
- p->primitives.gradients.rho[1] *= ihdim;
- p->primitives.gradients.rho[2] *= ihdim;
-
- p->primitives.gradients.v[0][0] *= ihdim;
- p->primitives.gradients.v[0][1] *= ihdim;
- p->primitives.gradients.v[0][2] *= ihdim;
- p->primitives.gradients.v[1][0] *= ihdim;
- p->primitives.gradients.v[1][1] *= ihdim;
- p->primitives.gradients.v[1][2] *= ihdim;
- p->primitives.gradients.v[2][0] *= ihdim;
- p->primitives.gradients.v[2][1] *= ihdim;
- p->primitives.gradients.v[2][2] *= ihdim;
-
- p->primitives.gradients.P[0] *= ihdim;
- p->primitives.gradients.P[1] *= ihdim;
- p->primitives.gradients.P[2] *= ihdim;
+ float norm;
+ if (p->geometry.wcorr > const_gizmo_min_wcorr) {
+ norm = ihdim;
} else {
-
- /* finalize gradients by multiplying with volume */
- p->primitives.gradients.rho[0] *= ihdimp1 * volume;
- p->primitives.gradients.rho[1] *= ihdimp1 * volume;
- p->primitives.gradients.rho[2] *= ihdimp1 * volume;
-
- p->primitives.gradients.v[0][0] *= ihdimp1 * volume;
- p->primitives.gradients.v[0][1] *= ihdimp1 * volume;
- p->primitives.gradients.v[0][2] *= ihdimp1 * volume;
-
- p->primitives.gradients.v[1][0] *= ihdimp1 * volume;
- p->primitives.gradients.v[1][1] *= ihdimp1 * volume;
- p->primitives.gradients.v[1][2] *= ihdimp1 * volume;
- p->primitives.gradients.v[2][0] *= ihdimp1 * volume;
- p->primitives.gradients.v[2][1] *= ihdimp1 * volume;
- p->primitives.gradients.v[2][2] *= ihdimp1 * volume;
-
- p->primitives.gradients.P[0] *= ihdimp1 * volume;
- p->primitives.gradients.P[1] *= ihdimp1 * volume;
- p->primitives.gradients.P[2] *= ihdimp1 * volume;
+ const float ihdimp1 = pow_dimension_plus_one(h_inv);
+ norm = ihdimp1 * volume;
}
+ p->gradients.rho[0] *= norm;
+ p->gradients.rho[1] *= norm;
+ p->gradients.rho[2] *= norm;
+
+ p->gradients.v[0][0] *= norm;
+ p->gradients.v[0][1] *= norm;
+ p->gradients.v[0][2] *= norm;
+ p->gradients.v[1][0] *= norm;
+ p->gradients.v[1][1] *= norm;
+ p->gradients.v[1][2] *= norm;
+ p->gradients.v[2][0] *= norm;
+ p->gradients.v[2][1] *= norm;
+ p->gradients.v[2][2] *= norm;
+
+ p->gradients.P[0] *= norm;
+ p->gradients.P[1] *= norm;
+ p->gradients.P[2] *= norm;
+
hydro_slope_limit_cell(p);
}
diff --git a/src/hydro/GizmoMFM/hydro_gradients_sph.h b/src/hydro/GizmoMFM/hydro_gradients_sph.h
index 169bed74f0b1b7e966f9880248f811d100bec13b..58233c3b75b22c4ba6347bb6e5db5e4e0e2ebe71 100644
--- a/src/hydro/GizmoMFM/hydro_gradients_sph.h
+++ b/src/hydro/GizmoMFM/hydro_gradients_sph.h
@@ -28,24 +28,24 @@
__attribute__((always_inline)) INLINE static void hydro_gradients_init(
struct part *p) {
- p->primitives.gradients.rho[0] = 0.0f;
- p->primitives.gradients.rho[1] = 0.0f;
- p->primitives.gradients.rho[2] = 0.0f;
+ p->gradients.rho[0] = 0.0f;
+ p->gradients.rho[1] = 0.0f;
+ p->gradients.rho[2] = 0.0f;
- p->primitives.gradients.v[0][0] = 0.0f;
- p->primitives.gradients.v[0][1] = 0.0f;
- p->primitives.gradients.v[0][2] = 0.0f;
+ p->gradients.v[0][0] = 0.0f;
+ p->gradients.v[0][1] = 0.0f;
+ p->gradients.v[0][2] = 0.0f;
- p->primitives.gradients.v[1][0] = 0.0f;
- p->primitives.gradients.v[1][1] = 0.0f;
- p->primitives.gradients.v[1][2] = 0.0f;
- p->primitives.gradients.v[2][0] = 0.0f;
- p->primitives.gradients.v[2][1] = 0.0f;
- p->primitives.gradients.v[2][2] = 0.0f;
+ p->gradients.v[1][0] = 0.0f;
+ p->gradients.v[1][1] = 0.0f;
+ p->gradients.v[1][2] = 0.0f;
+ p->gradients.v[2][0] = 0.0f;
+ p->gradients.v[2][1] = 0.0f;
+ p->gradients.v[2][2] = 0.0f;
- p->primitives.gradients.P[0] = 0.0f;
- p->primitives.gradients.P[1] = 0.0f;
- p->primitives.gradients.P[2] = 0.0f;
+ p->gradients.P[0] = 0.0f;
+ p->gradients.P[1] = 0.0f;
+ p->gradients.P[2] = 0.0f;
hydro_slope_limit_cell_init(p);
}
@@ -64,7 +64,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj) {
- const float r_inv = 1.f / sqrtf(r2);
+ const float r_inv = 1.0f / sqrtf(r2);
const float r = r2 * r_inv;
float wi, wi_dx;
@@ -72,41 +72,29 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
- /* very basic gradient estimate */
- pi->primitives.gradients.rho[0] -=
- wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
- pi->primitives.gradients.rho[1] -=
- wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
- pi->primitives.gradients.rho[2] -=
- wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-
- pi->primitives.gradients.v[0][0] -=
- wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
- pi->primitives.gradients.v[0][1] -=
- wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
- pi->primitives.gradients.v[0][2] -=
- wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-
- pi->primitives.gradients.v[1][0] -=
- wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pi->primitives.gradients.v[1][1] -=
- wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pi->primitives.gradients.v[1][2] -=
- wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
-
- pi->primitives.gradients.v[2][0] -=
- wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
- pi->primitives.gradients.v[2][1] -=
- wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
- pi->primitives.gradients.v[2][2] -=
- wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-
- pi->primitives.gradients.P[0] -=
- wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
- pi->primitives.gradients.P[1] -=
- wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
- pi->primitives.gradients.P[2] -=
- wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ const float dW[5] = {pi->rho - pj->rho, pi->v[0] - pj->v[0],
+ pi->v[1] - pj->v[1], pi->v[2] - pj->v[2], pi->P - pj->P};
+
+ const float normi = wi_dx * r_inv;
+ const float nidx[3] = {normi * dx[0], normi * dx[1], normi * dx[2]};
+
+ pi->gradients.rho[0] -= dW[0] * nidx[0];
+ pi->gradients.rho[1] -= dW[0] * nidx[1];
+ pi->gradients.rho[2] -= dW[0] * nidx[2];
+
+ pi->gradients.v[0][0] -= dW[1] * nidx[0];
+ pi->gradients.v[0][1] -= dW[1] * nidx[1];
+ pi->gradients.v[0][2] -= dW[1] * nidx[2];
+ pi->gradients.v[1][0] -= dW[2] * nidx[0];
+ pi->gradients.v[1][1] -= dW[2] * nidx[1];
+ pi->gradients.v[1][2] -= dW[2] * nidx[2];
+ pi->gradients.v[2][0] -= dW[3] * nidx[0];
+ pi->gradients.v[2][1] -= dW[3] * nidx[1];
+ pi->gradients.v[2][2] -= dW[3] * nidx[2];
+
+ pi->gradients.P[0] -= dW[4] * nidx[0];
+ pi->gradients.P[1] -= dW[4] * nidx[1];
+ pi->gradients.P[2] -= dW[4] * nidx[2];
hydro_slope_limit_cell_collect(pi, pj, r);
@@ -115,40 +103,27 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
+ const float normj = wj_dx * r_inv;
+ const float njdx[3] = {normj * dx[0], normj * dx[1], normj * dx[2]};
+
/* signs are the same as before, since we swap i and j twice */
- pj->primitives.gradients.rho[0] -=
- wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
- pj->primitives.gradients.rho[1] -=
- wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
- pj->primitives.gradients.rho[2] -=
- wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-
- pj->primitives.gradients.v[0][0] -=
- wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
- pj->primitives.gradients.v[0][1] -=
- wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
- pj->primitives.gradients.v[0][2] -=
- wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-
- pj->primitives.gradients.v[1][0] -=
- wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pj->primitives.gradients.v[1][1] -=
- wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pj->primitives.gradients.v[1][2] -=
- wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pj->primitives.gradients.v[2][0] -=
- wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
- pj->primitives.gradients.v[2][1] -=
- wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
- pj->primitives.gradients.v[2][2] -=
- wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-
- pj->primitives.gradients.P[0] -=
- wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
- pj->primitives.gradients.P[1] -=
- wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
- pj->primitives.gradients.P[2] -=
- wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ pj->gradients.rho[0] -= dW[0] * njdx[0];
+ pj->gradients.rho[1] -= dW[0] * njdx[1];
+ pj->gradients.rho[2] -= dW[0] * njdx[2];
+
+ pj->gradients.v[0][0] -= dW[1] * njdx[0];
+ pj->gradients.v[0][1] -= dW[1] * njdx[1];
+ pj->gradients.v[0][2] -= dW[1] * njdx[2];
+ pj->gradients.v[1][0] -= dW[2] * njdx[0];
+ pj->gradients.v[1][1] -= dW[2] * njdx[1];
+ pj->gradients.v[1][2] -= dW[2] * njdx[2];
+ pj->gradients.v[2][0] -= dW[3] * njdx[0];
+ pj->gradients.v[2][1] -= dW[3] * njdx[1];
+ pj->gradients.v[2][2] -= dW[3] * njdx[2];
+
+ pj->gradients.P[0] -= dW[4] * njdx[0];
+ pj->gradients.P[1] -= dW[4] * njdx[1];
+ pj->gradients.P[2] -= dW[4] * njdx[2];
hydro_slope_limit_cell_collect(pj, pi, r);
}
@@ -169,7 +144,7 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
struct part *restrict pi,
struct part *restrict pj) {
- const float r_inv = 1.f / sqrtf(r2);
+ const float r_inv = 1.0f / sqrtf(r2);
const float r = r2 * r_inv;
float wi, wi_dx;
@@ -177,41 +152,29 @@ hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
- /* very basic gradient estimate */
- pi->primitives.gradients.rho[0] -=
- wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
- pi->primitives.gradients.rho[1] -=
- wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
- pi->primitives.gradients.rho[2] -=
- wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
-
- pi->primitives.gradients.v[0][0] -=
- wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
- pi->primitives.gradients.v[0][1] -=
- wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
- pi->primitives.gradients.v[0][2] -=
- wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
-
- pi->primitives.gradients.v[1][0] -=
- wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pi->primitives.gradients.v[1][1] -=
- wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
- pi->primitives.gradients.v[1][2] -=
- wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
-
- pi->primitives.gradients.v[2][0] -=
- wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
- pi->primitives.gradients.v[2][1] -=
- wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
- pi->primitives.gradients.v[2][2] -=
- wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
-
- pi->primitives.gradients.P[0] -=
- wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
- pi->primitives.gradients.P[1] -=
- wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
- pi->primitives.gradients.P[2] -=
- wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ const float dW[5] = {pi->rho - pj->rho, pi->v[0] - pj->v[0],
+ pi->v[1] - pj->v[1], pi->v[2] - pj->v[2], pi->P - pj->P};
+
+ const float normi = wi_dx * r_inv;
+ const float nidx[3] = {normi * dx[0], normi * dx[1], normi * dx[2]};
+
+ pi->gradients.rho[0] -= dW[0] * nidx[0];
+ pi->gradients.rho[1] -= dW[0] * nidx[1];
+ pi->gradients.rho[2] -= dW[0] * nidx[2];
+
+ pi->gradients.v[0][0] -= dW[1] * nidx[0];
+ pi->gradients.v[0][1] -= dW[1] * nidx[1];
+ pi->gradients.v[0][2] -= dW[1] * nidx[2];
+ pi->gradients.v[1][0] -= dW[2] * nidx[0];
+ pi->gradients.v[1][1] -= dW[2] * nidx[1];
+ pi->gradients.v[1][2] -= dW[2] * nidx[2];
+ pi->gradients.v[2][0] -= dW[3] * nidx[0];
+ pi->gradients.v[2][1] -= dW[3] * nidx[1];
+ pi->gradients.v[2][2] -= dW[3] * nidx[2];
+
+ pi->gradients.P[0] -= dW[4] * nidx[0];
+ pi->gradients.P[1] -= dW[4] * nidx[1];
+ pi->gradients.P[2] -= dW[4] * nidx[2];
hydro_slope_limit_cell_collect(pi, pj, r);
}
@@ -229,26 +192,28 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
const float ihdimp1 = pow_dimension_plus_one(ih);
const float volume = p->geometry.volume;
+ const float norm = ihdimp1 * volume;
+
/* finalize gradients by multiplying with volume */
- p->primitives.gradients.rho[0] *= ihdimp1 * volume;
- p->primitives.gradients.rho[1] *= ihdimp1 * volume;
- p->primitives.gradients.rho[2] *= ihdimp1 * volume;
+ p->gradients.rho[0] *= norm;
+ p->gradients.rho[1] *= norm;
+ p->gradients.rho[2] *= norm;
- p->primitives.gradients.v[0][0] *= ihdimp1 * volume;
- p->primitives.gradients.v[0][1] *= ihdimp1 * volume;
- p->primitives.gradients.v[0][2] *= ihdimp1 * volume;
+ p->gradients.v[0][0] *= norm;
+ p->gradients.v[0][1] *= norm;
+ p->gradients.v[0][2] *= norm;
- p->primitives.gradients.v[1][0] *= ihdimp1 * volume;
- p->primitives.gradients.v[1][1] *= ihdimp1 * volume;
- p->primitives.gradients.v[1][2] *= ihdimp1 * volume;
+ p->gradients.v[1][0] *= norm;
+ p->gradients.v[1][1] *= norm;
+ p->gradients.v[1][2] *= norm;
- p->primitives.gradients.v[2][0] *= ihdimp1 * volume;
- p->primitives.gradients.v[2][1] *= ihdimp1 * volume;
- p->primitives.gradients.v[2][2] *= ihdimp1 * volume;
+ p->gradients.v[2][0] *= norm;
+ p->gradients.v[2][1] *= norm;
+ p->gradients.v[2][2] *= norm;
- p->primitives.gradients.P[0] *= ihdimp1 * volume;
- p->primitives.gradients.P[1] *= ihdimp1 * volume;
- p->primitives.gradients.P[2] *= ihdimp1 * volume;
+ p->gradients.P[0] *= norm;
+ p->gradients.P[1] *= norm;
+ p->gradients.P[2] *= norm;
hydro_slope_limit_cell(p);
}
diff --git a/src/hydro/GizmoMFM/hydro_iact.h b/src/hydro/GizmoMFM/hydro_iact.h
index fb54202a86b3e83aad7b2759c18dd3ce548cda13..5bed20d7f894a76d5fe3642c7438dc03195e43d6 100644
--- a/src/hydro/GizmoMFM/hydro_iact.h
+++ b/src/hydro/GizmoMFM/hydro_iact.h
@@ -57,7 +57,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
const float r = sqrtf(r2);
/* Compute density of pi. */
- const float hi_inv = 1.f / hi;
+ const float hi_inv = 1.0f / hi;
const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
@@ -75,7 +75,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pi->geometry.centroid[2] -= dx[2] * wi;
/* Compute density of pj. */
- const float hj_inv = 1.f / hj;
+ const float hj_inv = 1.0f / hj;
const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
@@ -123,7 +123,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* Get r and h inverse. */
const float r = sqrtf(r2);
- const float hi_inv = 1.f / hi;
+ const float hi_inv = 1.0f / hi;
const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
@@ -220,7 +220,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj, int mode, float a, float H) {
- const float r_inv = 1.f / sqrtf(r2);
+ const float r_inv = 1.0f / sqrtf(r2);
const float r = r2 * r_inv;
/* Initialize local variables */
@@ -238,16 +238,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
const float Vi = pi->geometry.volume;
const float Vj = pj->geometry.volume;
float Wi[5], Wj[5];
- Wi[0] = pi->primitives.rho;
- Wi[1] = pi->primitives.v[0];
- Wi[2] = pi->primitives.v[1];
- Wi[3] = pi->primitives.v[2];
- Wi[4] = pi->primitives.P;
- Wj[0] = pj->primitives.rho;
- Wj[1] = pj->primitives.v[0];
- Wj[2] = pj->primitives.v[1];
- Wj[3] = pj->primitives.v[2];
- Wj[4] = pj->primitives.P;
+ Wi[0] = pi->rho;
+ Wi[1] = pi->v[0];
+ Wi[2] = pi->v[1];
+ Wi[3] = pi->v[2];
+ Wi[4] = pi->P;
+ Wj[0] = pj->rho;
+ Wj[1] = pj->v[0];
+ Wj[2] = pj->v[1];
+ Wj[3] = pj->v[2];
+ Wj[4] = pj->P;
/* calculate the maximal signal velocity */
float vmax;
@@ -258,21 +258,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
} else
vmax = 0.0f;
- float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2];
-
/* Velocity on the axis linking the particles */
- float dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
- (Wi[3] - Wj[3]) * dx[2];
- dvdotdx = min(dvdotdx, dvdr);
+ float dvdr = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
+ (Wi[3] - Wj[3]) * dx[2];
- /* We only care about this velocity for particles moving towards each others
+ /* We only care about this velocity for particles moving towards each other
*/
- dvdotdx = min(dvdotdx, 0.f);
+ const float dvdotdx = min(dvdr, 0.0f);
/* Get the signal velocity */
/* the magical factor 3 also appears in Gadget2 */
- vmax -= 3.f * dvdotdx * r_inv;
+ vmax -= 3.0f * dvdotdx * r_inv;
/* Store the signal velocity */
pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
@@ -280,14 +276,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* Compute kernel of pi. */
float wi, wi_dx;
- const float hi_inv = 1.f / hi;
+ const float hi_inv = 1.0f / hi;
const float hi_inv_dim = pow_dimension(hi_inv);
const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
/* Compute kernel of pj. */
float wj, wj_dx;
- const float hj_inv = 1.f / hj;
+ const float hj_inv = 1.0f / hj;
const float hj_inv_dim = pow_dimension(hj_inv);
const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
@@ -298,17 +294,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
const float wi_dr = hidp1 * wi_dx;
const float wj_dr = hjdp1 * wj_dx;
dvdr *= r_inv;
- if (pj->primitives.rho > 0.)
- pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
- if (mode == 1 && pi->primitives.rho > 0.)
- pj->force.h_dt -= pi->conserved.mass * dvdr / pi->primitives.rho * wj_dr;
+ if (pj->rho > 0.0f)
+ pi->force.h_dt -= pj->conserved.mass * dvdr / pj->rho * wi_dr;
+ if (mode == 1 && pi->rho > 0.0f)
+ pj->force.h_dt -= pi->conserved.mass * dvdr / pi->rho * wj_dr;
/* Compute (square of) area */
/* eqn. (7) */
float Anorm2 = 0.0f;
float A[3];
- if (pi->density.wcorr > const_gizmo_min_wcorr &&
- pj->density.wcorr > const_gizmo_min_wcorr) {
+ if (pi->geometry.wcorr > const_gizmo_min_wcorr &&
+ pj->geometry.wcorr > const_gizmo_min_wcorr) {
/* in principle, we use Vi and Vj as weights for the left and right
contributions to the generalized surface vector.
However, if Vi and Vj are very different (because they have very
@@ -342,10 +338,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* if the interface has no area, nothing happens and we return */
/* continuing results in dividing by zero and NaN's... */
- if (Anorm2 == 0.f) return;
+ if (Anorm2 == 0.0f) return;
/* Compute the area */
- const float Anorm_inv = 1. / sqrtf(Anorm2);
+ const float Anorm_inv = 1.0f / sqrtf(Anorm2);
const float Anorm = Anorm2 * Anorm_inv;
#ifdef SWIFT_DEBUG_CHECKS
@@ -407,34 +403,24 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
riemann_solve_for_middle_state_flux(Wi, Wj, n_unit, vij, totflux);
/* Multiply with the interface surface area */
- totflux[0] *= Anorm;
totflux[1] *= Anorm;
totflux[2] *= Anorm;
totflux[3] *= Anorm;
totflux[4] *= Anorm;
- /* Store mass flux */
- const float mflux_i = totflux[0];
- pi->gravity.mflux[0] += mflux_i * dx[0];
- pi->gravity.mflux[1] += mflux_i * dx[1];
- pi->gravity.mflux[2] += mflux_i * dx[2];
-
/* Update conserved variables */
+ /* We shamelessly exploit the fact that the mass flux is zero and omit all
+ terms involving it */
/* eqn. (16) */
- 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];
+ pi->flux.momentum[0] -= totflux[1];
+ pi->flux.momentum[1] -= totflux[2];
+ pi->flux.momentum[2] -= totflux[3];
+ pi->flux.energy -= totflux[4];
#ifndef GIZMO_TOTAL_ENERGY
- const float ekin_i = 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 += 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_i;
+ pi->flux.energy += totflux[1] * pi->v[0];
+ pi->flux.energy += totflux[2] * pi->v[1];
+ pi->flux.energy += totflux[3] * pi->v[2];
#endif
/* Note that this used to be much more complicated in early implementations of
@@ -443,26 +429,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
* conservation anymore and just assume the current fluxes are representative
* for the flux over the entire time step. */
if (mode == 1) {
- /* Store mass flux */
- const float mflux_j = totflux[0];
- pj->gravity.mflux[0] -= mflux_j * dx[0];
- pj->gravity.mflux[1] -= mflux_j * dx[1];
- pj->gravity.mflux[2] -= mflux_j * dx[2];
-
- 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];
+ pj->flux.momentum[0] += totflux[1];
+ pj->flux.momentum[1] += totflux[2];
+ pj->flux.momentum[2] += totflux[3];
+ pj->flux.energy += totflux[4];
#ifndef GIZMO_TOTAL_ENERGY
- const float ekin_j = 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 -= 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_j;
+ pj->flux.energy -= totflux[1] * pj->v[0];
+ pj->flux.energy -= totflux[2] * pj->v[1];
+ pj->flux.energy -= totflux[3] * pj->v[2];
#endif
}
}
diff --git a/src/hydro/GizmoMFM/hydro_io.h b/src/hydro/GizmoMFM/hydro_io.h
index 59d579f70cd4aedc728dbf42038eff78d4c507d5..1f956edf3fdc31990c6aba254603ea69a98238eb 100644
--- a/src/hydro/GizmoMFM/hydro_io.h
+++ b/src/hydro/GizmoMFM/hydro_io.h
@@ -63,7 +63,7 @@ INLINE static void hydro_read_particles(struct part* parts,
list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
UNIT_CONV_ACCELERATION, parts, a_hydro);
list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
- UNIT_CONV_DENSITY, parts, primitives.rho);
+ UNIT_CONV_DENSITY, parts, rho);
}
/**
@@ -203,12 +203,12 @@ INLINE static void hydro_write_particles(const struct part* parts,
parts, xparts, convert_u);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
- list[6] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
- primitives.rho);
+ list[6] =
+ io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
list[7] = io_make_output_field_convert_part(
"Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, xparts, convert_A);
- list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
- parts, primitives.P);
+ list[8] =
+ io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, P);
list[9] = io_make_output_field_convert_part(
"TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, xparts, convert_Etot);
diff --git a/src/hydro/GizmoMFM/hydro_part.h b/src/hydro/GizmoMFM/hydro_part.h
index 857c429ec933ad3eea730e8f0b6f830782cdf77b..0055d7d86a35746a8ba90015b3a6986f8ddb5f9f 100644
--- a/src/hydro/GizmoMFM/hydro_part.h
+++ b/src/hydro/GizmoMFM/hydro_part.h
@@ -63,31 +63,23 @@ struct part {
/* Particle smoothing length. */
float h;
- /* The primitive hydrodynamical variables. */
- struct {
-
- /* Density. */
- float rho;
-
- /* Fluid velocity. */
- float v[3];
+ /* Density. */
+ float rho;
- /* Pressure. */
- float P;
+ /* Pressure. */
+ float P;
- /* Gradients of the primitive variables. */
+ union {
+ /* Quantities used during the volume (=density) loop. */
struct {
- /* Density gradients. */
- float rho[3];
+ /* Derivative of particle number density. */
+ float wcount_dh;
- /* Fluid velocity gradients. */
- float v[3][3];
+ /* Particle number density. */
+ float wcount;
- /* Pressure gradients. */
- float P[3];
-
- } gradients;
+ } density;
/* Quantities needed by the slope limiter. */
struct {
@@ -106,7 +98,56 @@ struct part {
} limiter;
- } primitives;
+ struct {
+ /* Fluxes. */
+ struct {
+
+ /* No mass flux, since it is always zero. */
+
+ /* Momentum flux. */
+ float momentum[3];
+
+ /* Energy flux. */
+ float energy;
+
+ } flux;
+
+ /* Variables used for timestep calculation. */
+ struct {
+
+ /* 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;
+
+ /* Quantities used during the force loop. */
+ struct {
+
+ /* Needed to drift the primitive variables. */
+ float h_dt;
+
+ } force;
+ };
+ };
+
+ /* Gradients of the primitive variables. */
+ struct {
+
+ /* Density gradients. */
+ float rho[3];
+
+ /* Fluid velocity gradients. */
+ float v[3][3];
+
+ /* Pressure gradients. */
+ float P[3];
+
+ } gradients;
/* The conserved hydrodynamical variables. */
struct {
@@ -120,20 +161,6 @@ struct part {
/* Fluid thermal energy (not per unit mass!). */
float energy;
- /* Fluxes. */
- struct {
-
- /* Mass flux. */
- float mass;
-
- /* Momentum flux. */
- float momentum[3];
-
- /* Energy flux. */
- float energy;
-
- } flux;
-
} conserved;
/* Geometrical quantities used for hydro. */
@@ -149,48 +176,10 @@ struct part {
/* Centroid of the "cell". */
float centroid[3];
- } geometry;
-
- /* Variables used for timestep calculation. */
- struct {
-
- /* 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;
-
- /* Quantities used during the volume (=density) loop. */
- struct {
-
- /* Derivative of particle number density. */
- float wcount_dh;
-
- /* Particle number density. */
- float wcount;
-
/* Correction factor for wcount. */
float wcorr;
- } density;
-
- /* Quantities used during the force loop. */
- struct {
-
- /* Needed to drift the primitive variables. */
- float h_dt;
-
- } force;
-
- /* Specific stuff for the gravity-hydro coupling. */
- struct {
-
- /* Current value of the mass flux vector. */
- float mflux[3];
-
- } gravity;
+ } geometry;
/* Chemistry information */
struct chemistry_part_data chemistry_data;
diff --git a/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h b/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h
index 7dec6f499da31de1f10652a31781a788166957cc..329ca5fda4510a3672bfe698bed52c8f0bebc895 100644
--- a/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h
+++ b/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h
@@ -29,18 +29,18 @@
__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
struct part* p) {
- p->primitives.limiter.rho[0] = FLT_MAX;
- p->primitives.limiter.rho[1] = -FLT_MAX;
- p->primitives.limiter.v[0][0] = FLT_MAX;
- p->primitives.limiter.v[0][1] = -FLT_MAX;
- p->primitives.limiter.v[1][0] = FLT_MAX;
- p->primitives.limiter.v[1][1] = -FLT_MAX;
- p->primitives.limiter.v[2][0] = FLT_MAX;
- p->primitives.limiter.v[2][1] = -FLT_MAX;
- p->primitives.limiter.P[0] = FLT_MAX;
- p->primitives.limiter.P[1] = -FLT_MAX;
-
- p->primitives.limiter.maxr = -FLT_MAX;
+ p->limiter.rho[0] = FLT_MAX;
+ p->limiter.rho[1] = -FLT_MAX;
+ p->limiter.v[0][0] = FLT_MAX;
+ p->limiter.v[0][1] = -FLT_MAX;
+ p->limiter.v[1][0] = FLT_MAX;
+ p->limiter.v[1][1] = -FLT_MAX;
+ p->limiter.v[2][0] = FLT_MAX;
+ p->limiter.v[2][1] = -FLT_MAX;
+ p->limiter.P[0] = FLT_MAX;
+ p->limiter.P[1] = -FLT_MAX;
+
+ p->limiter.maxr = -FLT_MAX;
}
/**
@@ -56,30 +56,20 @@ hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
/* basic slope limiter: collect the maximal and the minimal value for the
* primitive variables among the ngbs */
- pi->primitives.limiter.rho[0] =
- min(pj->primitives.rho, pi->primitives.limiter.rho[0]);
- pi->primitives.limiter.rho[1] =
- max(pj->primitives.rho, pi->primitives.limiter.rho[1]);
-
- pi->primitives.limiter.v[0][0] =
- min(pj->primitives.v[0], pi->primitives.limiter.v[0][0]);
- pi->primitives.limiter.v[0][1] =
- max(pj->primitives.v[0], pi->primitives.limiter.v[0][1]);
- pi->primitives.limiter.v[1][0] =
- min(pj->primitives.v[1], pi->primitives.limiter.v[1][0]);
- pi->primitives.limiter.v[1][1] =
- max(pj->primitives.v[1], pi->primitives.limiter.v[1][1]);
- pi->primitives.limiter.v[2][0] =
- min(pj->primitives.v[2], pi->primitives.limiter.v[2][0]);
- pi->primitives.limiter.v[2][1] =
- max(pj->primitives.v[2], pi->primitives.limiter.v[2][1]);
-
- pi->primitives.limiter.P[0] =
- min(pj->primitives.P, pi->primitives.limiter.P[0]);
- pi->primitives.limiter.P[1] =
- max(pj->primitives.P, pi->primitives.limiter.P[1]);
-
- pi->primitives.limiter.maxr = max(r, pi->primitives.limiter.maxr);
+ pi->limiter.rho[0] = min(pj->rho, pi->limiter.rho[0]);
+ pi->limiter.rho[1] = max(pj->rho, pi->limiter.rho[1]);
+
+ pi->limiter.v[0][0] = min(pj->v[0], pi->limiter.v[0][0]);
+ pi->limiter.v[0][1] = max(pj->v[0], pi->limiter.v[0][1]);
+ pi->limiter.v[1][0] = min(pj->v[1], pi->limiter.v[1][0]);
+ pi->limiter.v[1][1] = max(pj->v[1], pi->limiter.v[1][1]);
+ pi->limiter.v[2][0] = min(pj->v[2], pi->limiter.v[2][0]);
+ pi->limiter.v[2][1] = max(pj->v[2], pi->limiter.v[2][1]);
+
+ pi->limiter.P[0] = min(pj->P, pi->limiter.P[0]);
+ pi->limiter.P[1] = max(pj->P, pi->limiter.P[1]);
+
+ pi->limiter.maxr = max(r, pi->limiter.maxr);
}
/**
@@ -92,94 +82,94 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
float gradtrue, gradrho[3], gradv[3][3], gradP[3];
- gradrho[0] = p->primitives.gradients.rho[0];
- gradrho[1] = p->primitives.gradients.rho[1];
- gradrho[2] = p->primitives.gradients.rho[2];
+ gradrho[0] = p->gradients.rho[0];
+ gradrho[1] = p->gradients.rho[1];
+ gradrho[2] = p->gradients.rho[2];
- gradv[0][0] = p->primitives.gradients.v[0][0];
- gradv[0][1] = p->primitives.gradients.v[0][1];
- gradv[0][2] = p->primitives.gradients.v[0][2];
+ gradv[0][0] = p->gradients.v[0][0];
+ gradv[0][1] = p->gradients.v[0][1];
+ gradv[0][2] = p->gradients.v[0][2];
- gradv[1][0] = p->primitives.gradients.v[1][0];
- gradv[1][1] = p->primitives.gradients.v[1][1];
- gradv[1][2] = p->primitives.gradients.v[1][2];
+ gradv[1][0] = p->gradients.v[1][0];
+ gradv[1][1] = p->gradients.v[1][1];
+ gradv[1][2] = p->gradients.v[1][2];
- gradv[2][0] = p->primitives.gradients.v[2][0];
- gradv[2][1] = p->primitives.gradients.v[2][1];
- gradv[2][2] = p->primitives.gradients.v[2][2];
+ gradv[2][0] = p->gradients.v[2][0];
+ gradv[2][1] = p->gradients.v[2][1];
+ gradv[2][2] = p->gradients.v[2][2];
- gradP[0] = p->primitives.gradients.P[0];
- gradP[1] = p->primitives.gradients.P[1];
- gradP[2] = p->primitives.gradients.P[2];
+ gradP[0] = p->gradients.P[0];
+ gradP[1] = p->gradients.P[1];
+ gradP[2] = p->gradients.P[2];
gradtrue = sqrtf(gradrho[0] * gradrho[0] + gradrho[1] * gradrho[1] +
gradrho[2] * gradrho[2]);
if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- const float gradmax = p->primitives.limiter.rho[1] - p->primitives.rho;
- const float gradmin = p->primitives.rho - p->primitives.limiter.rho[0];
- const float gradtrue_inv = 1.f / gradtrue;
+ gradtrue *= p->limiter.maxr;
+ const float gradmax = p->limiter.rho[1] - p->rho;
+ const float gradmin = p->rho - p->limiter.rho[0];
+ const float gradtrue_inv = 1.0f / gradtrue;
const float alpha =
min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
- p->primitives.gradients.rho[0] *= alpha;
- p->primitives.gradients.rho[1] *= alpha;
- p->primitives.gradients.rho[2] *= alpha;
+ p->gradients.rho[0] *= alpha;
+ p->gradients.rho[1] *= alpha;
+ p->gradients.rho[2] *= alpha;
}
gradtrue = sqrtf(gradv[0][0] * gradv[0][0] + gradv[0][1] * gradv[0][1] +
gradv[0][2] * gradv[0][2]);
if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- const float gradmax = p->primitives.limiter.v[0][1] - p->primitives.v[0];
- const float gradmin = p->primitives.v[0] - p->primitives.limiter.v[0][0];
- const float gradtrue_inv = 1.f / gradtrue;
+ gradtrue *= p->limiter.maxr;
+ const float gradmax = p->limiter.v[0][1] - p->v[0];
+ const float gradmin = p->v[0] - p->limiter.v[0][0];
+ const float gradtrue_inv = 1.0f / gradtrue;
const float alpha =
min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
- p->primitives.gradients.v[0][0] *= alpha;
- p->primitives.gradients.v[0][1] *= alpha;
- p->primitives.gradients.v[0][2] *= alpha;
+ p->gradients.v[0][0] *= alpha;
+ p->gradients.v[0][1] *= alpha;
+ p->gradients.v[0][2] *= alpha;
}
gradtrue = sqrtf(gradv[1][0] * gradv[1][0] + gradv[1][1] * gradv[1][1] +
gradv[1][2] * gradv[1][2]);
if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- const float gradmax = p->primitives.limiter.v[1][1] - p->primitives.v[1];
- const float gradmin = p->primitives.v[1] - p->primitives.limiter.v[1][0];
- const float gradtrue_inv = 1.f / gradtrue;
+ gradtrue *= p->limiter.maxr;
+ const float gradmax = p->limiter.v[1][1] - p->v[1];
+ const float gradmin = p->v[1] - p->limiter.v[1][0];
+ const float gradtrue_inv = 1.0f / gradtrue;
const float alpha =
min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
- p->primitives.gradients.v[1][0] *= alpha;
- p->primitives.gradients.v[1][1] *= alpha;
- p->primitives.gradients.v[1][2] *= alpha;
+ p->gradients.v[1][0] *= alpha;
+ p->gradients.v[1][1] *= alpha;
+ p->gradients.v[1][2] *= alpha;
}
gradtrue = sqrtf(gradv[2][0] * gradv[2][0] + gradv[2][1] * gradv[2][1] +
gradv[2][2] * gradv[2][2]);
if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- const float gradmax = p->primitives.limiter.v[2][1] - p->primitives.v[2];
- const float gradmin = p->primitives.v[2] - p->primitives.limiter.v[2][0];
- const float gradtrue_inv = 1.f / gradtrue;
+ gradtrue *= p->limiter.maxr;
+ const float gradmax = p->limiter.v[2][1] - p->v[2];
+ const float gradmin = p->v[2] - p->limiter.v[2][0];
+ const float gradtrue_inv = 1.0f / gradtrue;
const float alpha =
min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
- p->primitives.gradients.v[2][0] *= alpha;
- p->primitives.gradients.v[2][1] *= alpha;
- p->primitives.gradients.v[2][2] *= alpha;
+ p->gradients.v[2][0] *= alpha;
+ p->gradients.v[2][1] *= alpha;
+ p->gradients.v[2][2] *= alpha;
}
gradtrue =
sqrtf(gradP[0] * gradP[0] + gradP[1] * gradP[1] + gradP[2] * gradP[2]);
if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- const float gradmax = p->primitives.limiter.P[1] - p->primitives.P;
- const float gradmin = p->primitives.P - p->primitives.limiter.P[0];
- const float gradtrue_inv = 1.f / gradtrue;
+ gradtrue *= p->limiter.maxr;
+ const float gradmax = p->limiter.P[1] - p->P;
+ const float gradmin = p->P - p->limiter.P[0];
+ const float gradtrue_inv = 1.0f / gradtrue;
const float alpha =
min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
- p->primitives.gradients.P[0] *= alpha;
- p->primitives.gradients.P[1] *= alpha;
- p->primitives.gradients.P[2] *= alpha;
+ p->gradients.P[0] *= alpha;
+ p->gradients.P[1] *= alpha;
+ p->gradients.P[2] *= alpha;
}
}
diff --git a/src/hydro/GizmoMFM/hydro_slope_limiters_face.h b/src/hydro/GizmoMFM/hydro_slope_limiters_face.h
index 8f0636c992ccbd615cc62ca09c1b0ca9b08ea56b..a7c0e37d3953d1d385128b4a1c72dd2311ac00ae 100644
--- a/src/hydro/GizmoMFM/hydro_slope_limiters_face.h
+++ b/src/hydro/GizmoMFM/hydro_slope_limiters_face.h
@@ -56,15 +56,17 @@ hydro_slope_limit_face_quantity(float phi_i, float phi_j, float phi_mid0,
float phiplus, phiminus, phi_mid;
- if (same_signf(phimax + delta1, phimax))
+ if (same_signf(phimax + delta1, phimax)) {
phiplus = phimax + delta1;
- else
+ } else {
phiplus = phimax / (1.0f + delta1 / (fabsf(phimax) + FLT_MIN));
+ }
- if (same_signf(phimin - delta1, phimin))
+ if (same_signf(phimin - delta1, phimin)) {
phiminus = phimin - delta1;
- else
+ } else {
phiminus = phimin / (1.0f + delta1 / (fabsf(phimin) + FLT_MIN));
+ }
if (phi_i < phi_j) {
const float temp = min(phibar + delta2, phi_mid0);
diff --git a/src/hydro/GizmoMFM/hydro_velocities.h b/src/hydro/GizmoMFM/hydro_velocities.h
deleted file mode 100644
index ea30d5a6c9c74d34ffd73fa6ab941640f37b02e4..0000000000000000000000000000000000000000
--- a/src/hydro/GizmoMFM/hydro_velocities.h
+++ /dev/null
@@ -1,157 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Coypright (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_VELOCITIES_H
-#define SWIFT_HYDRO_VELOCITIES_H
-
-/**
- * @brief Initialize the GIZMO particle velocities before the start of the
- * actual run based on the initial value of the primitive velocity.
- *
- * @param p The particle to act upon.
- * @param xp The extended particle data to act upon.
- */
-__attribute__((always_inline)) INLINE static void hydro_velocities_init(
- struct part* restrict p, struct xpart* restrict xp) {
-
-#ifdef GIZMO_FIX_PARTICLES
- p->v[0] = 0.f;
- p->v[1] = 0.f;
- p->v[2] = 0.f;
-#else
- p->v[0] = p->primitives.v[0];
- p->v[1] = p->primitives.v[1];
- p->v[2] = p->primitives.v[2];
-#endif
-
- xp->v_full[0] = p->v[0];
- xp->v_full[1] = p->v[1];
- xp->v_full[2] = p->v[2];
-}
-
-/**
- * @brief Set the particle velocity field that will be used to deboost fluid
- * velocities during the force loop.
- *
- * @param p The particle 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) {}
-
-/**
- * @brief Set the variables that will be used to update the smoothing length
- * during the drift (these will depend on the movement of the particles).
- *
- * @param p The particle to act upon.
- */
-__attribute__((always_inline)) INLINE static void hydro_velocities_end_force(
- struct part* restrict p) {
-
-#ifdef GIZMO_FIX_PARTICLES
- /* disable the smoothing length update, since the smoothing lengths should
- stay the same for all steps (particles don't move) */
- p->force.h_dt = 0.0f;
-#else
- /* Add normalization to h_dt. */
- p->force.h_dt *= p->h * hydro_dimension_inv;
-#endif
-}
-
-/**
- * @brief Set the velocity of a GIZMO particle, based on the values of its
- * primitive variables and the geometry of its mesh-free "cell".
- *
- * @param p The particle to act upon.
- * @param xp The extended particle data to act upon.
- */
-__attribute__((always_inline)) INLINE static void hydro_velocities_set(
- struct part* restrict p, struct xpart* restrict xp) {
-
-/* We first set the particle velocity. */
-#ifdef GIZMO_FIX_PARTICLES
-
- p->v[0] = 0.f;
- p->v[1] = 0.f;
- p->v[2] = 0.f;
-
-#else // GIZMO_FIX_PARTICLES
-
- if (p->conserved.mass > 0.f && p->primitives.rho > 0.f) {
-
- const float inverse_mass = 1.f / p->conserved.mass;
-
- /* Normal case: set particle velocity to fluid velocity. */
- p->v[0] = p->conserved.momentum[0] * inverse_mass;
- p->v[1] = p->conserved.momentum[1] * inverse_mass;
- p->v[2] = p->conserved.momentum[2] * inverse_mass;
-
-#ifdef GIZMO_STEER_MOTION
-
- /* Add a correction to the velocity to keep particle positions close enough
- to
- the centroid of their mesh-free "cell". */
- /* The correction term below is the same one described in Springel (2010).
- */
- float ds[3];
- ds[0] = p->geometry.centroid[0];
- ds[1] = p->geometry.centroid[1];
- ds[2] = p->geometry.centroid[2];
- const float d = sqrtf(ds[0] * ds[0] + ds[1] * ds[1] + ds[2] * ds[2]);
- const float R = get_radius_dimension_sphere(p->geometry.volume);
- const float eta = 0.25f;
- const float etaR = eta * R;
- const float xi = 1.f;
- const float soundspeed =
- sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
- /* We only apply the correction if the offset between centroid and position
- is too large. */
- if (d > 0.9f * etaR) {
- float fac = xi * soundspeed / d;
- if (d < 1.1f * etaR) {
- fac *= 5.f * (d - 0.9f * etaR) / etaR;
- }
- p->v[0] -= ds[0] * fac;
- p->v[1] -= ds[1] * fac;
- p->v[2] -= ds[2] * fac;
- }
-
-#endif // GIZMO_STEER_MOTION
- } else {
- /* Vacuum particles have no fluid velocity. */
- p->v[0] = 0.f;
- p->v[1] = 0.f;
- p->v[2] = 0.f;
- }
-
-#endif // GIZMO_FIX_PARTICLES
-
- /* Now make sure all velocity variables are up to date. */
- xp->v_full[0] = p->v[0];
- xp->v_full[1] = p->v[1];
- xp->v_full[2] = p->v[2];
-
- if (p->gpart) {
- p->gpart->v_full[0] = p->v[0];
- p->gpart->v_full[1] = p->v[1];
- p->gpart->v_full[2] = p->v[2];
- }
-}
-
-#endif /* SWIFT_HYDRO_VELOCITIES_H */
diff --git a/src/riemann/riemann_hllc.h b/src/riemann/riemann_hllc.h
index d7e376fe667dd3cfeac58d5820053057ab928d3e..a34b907ba5c6142755f2c4107d7b456393c12b42 100644
--- a/src/riemann/riemann_hllc.h
+++ b/src/riemann/riemann_hllc.h
@@ -33,11 +33,6 @@
#include "riemann_checks.h"
#include "riemann_vacuum.h"
-#ifndef EOS_IDEAL_GAS
-#error \
- "The HLLC Riemann solver currently only supports and ideal gas equation of state. Either select this equation of state, or try using another Riemann solver!"
-#endif
-
__attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
const float *WL, const float *WR, const float *n, const float *vij,
float *totflux) {
@@ -48,19 +43,21 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
/* Handle pure vacuum */
if (!WL[0] && !WR[0]) {
- totflux[0] = 0.f;
- totflux[1] = 0.f;
- totflux[2] = 0.f;
- totflux[3] = 0.f;
- totflux[4] = 0.f;
+ totflux[0] = 0.0f;
+ totflux[1] = 0.0f;
+ totflux[2] = 0.0f;
+ totflux[3] = 0.0f;
+ totflux[4] = 0.0f;
return;
}
/* STEP 0: obtain velocity in interface frame */
const float uL = WL[1] * n[0] + WL[2] * n[1] + WL[3] * n[2];
const float uR = WR[1] * n[0] + WR[2] * n[1] + WR[3] * n[2];
- const float aL = sqrtf(hydro_gamma * WL[4] / WL[0]);
- const float aR = sqrtf(hydro_gamma * WR[4] / WR[0]);
+ const float rhoLinv = 1.0f / WL[0];
+ const float rhoRinv = 1.0f / WR[0];
+ const float aL = sqrtf(hydro_gamma * WL[4] * rhoLinv);
+ const float aR = sqrtf(hydro_gamma * WR[4] * rhoRinv);
/* Handle vacuum: vacuum does not require iteration and is always exact */
if (riemann_is_vacuum(WL, WR, uL, uR, aL, aR)) {
@@ -69,98 +66,90 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
}
/* STEP 1: pressure estimate */
- const float rhobar = 0.5f * (WL[0] + WR[0]);
- const float abar = 0.5f * (aL + aR);
- const float pPVRS = 0.5f * (WL[4] + WR[4]) - 0.5f * (uR - uL) * rhobar * abar;
- const float pstar = max(0.f, pPVRS);
+ const float rhobar = WL[0] + WR[0];
+ const float abar = aL + aR;
+ const float pPVRS =
+ 0.5f * ((WL[4] + WR[4]) - 0.25f * (uR - uL) * rhobar * abar);
+ const float pstar = max(0.0f, pPVRS);
/* STEP 2: wave speed estimates
all these speeds are along the interface normal, since uL and uR are */
- float qL = 1.f;
- if (pstar > WL[4] && WL[4] > 0.f) {
- qL = sqrtf(1.f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
- (pstar / WL[4] - 1.f));
+ float qL = 1.0f;
+ if (pstar > WL[4] && WL[4] > 0.0f) {
+ qL = sqrtf(1.0f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
+ (pstar / WL[4] - 1.0f));
}
- float qR = 1.f;
- if (pstar > WR[4] && WR[4] > 0.f) {
- qR = sqrtf(1.f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
- (pstar / WR[4] - 1.f));
+ float qR = 1.0f;
+ if (pstar > WR[4] && WR[4] > 0.0f) {
+ qR = sqrtf(1.0f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
+ (pstar / WR[4] - 1.0f));
}
- const float SL = uL - aL * qL;
- const float SR = uR + aR * qR;
+ const float SLmuL = -aL * qL;
+ const float SRmuR = aR * qR;
const float Sstar =
- (WR[4] - WL[4] + WL[0] * uL * (SL - uL) - WR[0] * uR * (SR - uR)) /
- (WL[0] * (SL - uL) - WR[0] * (SR - uR));
+ (WR[4] - WL[4] + WL[0] * uL * SLmuL - WR[0] * uR * SRmuR) /
+ (WL[0] * SLmuL - WR[0] * SRmuR);
/* STEP 3: HLLC flux in a frame moving with the interface velocity */
- if (Sstar >= 0.f) {
+ if (Sstar >= 0.0f) {
+ const float rhoLuL = WL[0] * uL;
+ const float v2 = WL[1] * WL[1] + WL[2] * WL[2] + WL[3] * WL[3];
+ const float eL =
+ WL[4] * rhoLinv * hydro_one_over_gamma_minus_one + 0.5f * v2;
+ const float SL = SLmuL + uL;
+
/* flux FL */
- totflux[0] = WL[0] * uL;
+ totflux[0] = rhoLuL;
/* these are the actual correct fluxes in the boosted lab frame
(not rotated to interface frame) */
- totflux[1] = WL[0] * WL[1] * uL + WL[4] * n[0];
- totflux[2] = WL[0] * WL[2] * uL + WL[4] * n[1];
- totflux[3] = WL[0] * WL[3] * uL + WL[4] * n[2];
- const float v2 = WL[1] * WL[1] + WL[2] * WL[2] + WL[3] * WL[3];
- const float eL = WL[4] * hydro_one_over_gamma_minus_one / WL[0] + 0.5f * v2;
- totflux[4] = WL[0] * eL * uL + WL[4] * uL;
- if (SL < 0.f) {
-
- float UstarL[5];
-
- /* add flux FstarL */
- UstarL[0] = 1.f;
- /* we need UstarL in the lab frame:
- * subtract the velocity in the interface frame from the lab frame
- * velocity and then add Sstar in interface frame */
- UstarL[1] = WL[1] + (Sstar - uL) * n[0];
- UstarL[2] = WL[2] + (Sstar - uL) * n[1];
- UstarL[3] = WL[3] + (Sstar - uL) * n[2];
- UstarL[4] = eL + (Sstar - uL) * (Sstar + WL[4] / (WL[0] * (SL - uL)));
- UstarL[0] *= WL[0] * (SL - uL) / (SL - Sstar);
- UstarL[1] *= WL[0] * (SL - uL) / (SL - Sstar);
- UstarL[2] *= WL[0] * (SL - uL) / (SL - Sstar);
- UstarL[3] *= WL[0] * (SL - uL) / (SL - Sstar);
- UstarL[4] *= WL[0] * (SL - uL) / (SL - Sstar);
- totflux[0] += SL * (UstarL[0] - WL[0]);
- totflux[1] += SL * (UstarL[1] - WL[0] * WL[1]);
- totflux[2] += SL * (UstarL[2] - WL[0] * WL[2]);
- totflux[3] += SL * (UstarL[3] - WL[0] * WL[3]);
- totflux[4] += SL * (UstarL[4] - WL[0] * eL);
+ totflux[1] = rhoLuL * WL[1] + WL[4] * n[0];
+ totflux[2] = rhoLuL * WL[2] + WL[4] * n[1];
+ totflux[3] = rhoLuL * WL[3] + WL[4] * n[2];
+ totflux[4] = rhoLuL * eL + WL[4] * uL;
+
+ if (SL < 0.0f) {
+
+ const float starfac = SLmuL / (SL - Sstar) - 1.0f;
+ const float rhoLSL = WL[0] * SL;
+ const float SstarmuL = Sstar - uL;
+ const float rhoLSLstarfac = rhoLSL * starfac;
+ const float rhoLSLSstarmuL = rhoLSL * SstarmuL;
+
+ totflux[0] += rhoLSLstarfac;
+ totflux[1] += rhoLSLstarfac * WL[1] + rhoLSLSstarmuL * n[0];
+ totflux[2] += rhoLSLstarfac * WL[2] + rhoLSLSstarmuL * n[1];
+ totflux[3] += rhoLSLstarfac * WL[3] + rhoLSLSstarmuL * n[2];
+ totflux[4] += rhoLSLstarfac * eL +
+ rhoLSLSstarmuL * (Sstar + WL[4] / (WL[0] * SLmuL));
}
} else {
- /* flux FR */
- totflux[0] = WR[0] * uR;
- totflux[1] = WR[0] * WR[1] * uR + WR[4] * n[0];
- totflux[2] = WR[0] * WR[2] * uR + WR[4] * n[1];
- totflux[3] = WR[0] * WR[3] * uR + WR[4] * n[2];
+ const float rhoRuR = WR[0] * uR;
const float v2 = WR[1] * WR[1] + WR[2] * WR[2] + WR[3] * WR[3];
- const float eR = WR[4] * hydro_one_over_gamma_minus_one / WR[0] + 0.5f * v2;
- totflux[4] = WR[0] * eR * uR + WR[4] * uR;
- if (SR > 0.f) {
+ const float eR =
+ WR[4] * rhoRinv * hydro_one_over_gamma_minus_one + 0.5f * v2;
+ const float SR = SRmuR + uR;
- float UstarR[5];
-
- /* add flux FstarR */
- UstarR[0] = 1.f;
-
- /* we need UstarR in the lab frame:
- * subtract the velocity in the interface frame from the lab frame
- * velocity and then add Sstar in interface frame */
- UstarR[1] = WR[1] + (Sstar - uR) * n[0];
- UstarR[2] = WR[2] + (Sstar - uR) * n[1];
- UstarR[3] = WR[3] + (Sstar - uR) * n[2];
- UstarR[4] = eR + (Sstar - uR) * (Sstar + WR[4] / (WR[0] * (SR - uR)));
- UstarR[0] *= WR[0] * (SR - uR) / (SR - Sstar);
- UstarR[1] *= WR[0] * (SR - uR) / (SR - Sstar);
- UstarR[2] *= WR[0] * (SR - uR) / (SR - Sstar);
- UstarR[3] *= WR[0] * (SR - uR) / (SR - Sstar);
- UstarR[4] *= WR[0] * (SR - uR) / (SR - Sstar);
- totflux[0] += SR * (UstarR[0] - WR[0]);
- totflux[1] += SR * (UstarR[1] - WR[0] * WR[1]);
- totflux[2] += SR * (UstarR[2] - WR[0] * WR[2]);
- totflux[3] += SR * (UstarR[3] - WR[0] * WR[3]);
- totflux[4] += SR * (UstarR[4] - WR[0] * eR);
+ /* flux FR */
+ totflux[0] = rhoRuR;
+ totflux[1] = rhoRuR * WR[1] + WR[4] * n[0];
+ totflux[2] = rhoRuR * WR[2] + WR[4] * n[1];
+ totflux[3] = rhoRuR * WR[3] + WR[4] * n[2];
+ totflux[4] = rhoRuR * eR + WR[4] * uR;
+
+ if (SR > 0.0f) {
+
+ const float starfac = SRmuR / (SR - Sstar) - 1.0f;
+ const float rhoRSR = WR[0] * SR;
+ const float SstarmuR = Sstar - uR;
+ const float rhoRSRstarfac = rhoRSR * starfac;
+ const float rhoRSRSstarmuR = rhoRSR * SstarmuR;
+
+ totflux[0] += rhoRSRstarfac;
+ totflux[1] += rhoRSRstarfac * WR[1] + rhoRSRSstarmuR * n[0];
+ totflux[2] += rhoRSRstarfac * WR[2] + rhoRSRSstarmuR * n[1];
+ totflux[3] += rhoRSRstarfac * WR[3] + rhoRSRSstarmuR * n[2];
+ totflux[4] += rhoRSRstarfac * eR +
+ rhoRSRSstarmuR * (Sstar + WR[4] / (WR[0] * SRmuR));
}
}
@@ -194,11 +183,11 @@ riemann_solve_for_middle_state_flux(const float *WL, const float *WR,
/* Handle pure vacuum */
if (!WL[0] && !WR[0]) {
- totflux[0] = 0.f;
- totflux[1] = 0.f;
- totflux[2] = 0.f;
- totflux[3] = 0.f;
- totflux[4] = 0.f;
+ totflux[0] = 0.0f;
+ totflux[1] = 0.0f;
+ totflux[2] = 0.0f;
+ totflux[3] = 0.0f;
+ totflux[4] = 0.0f;
return;
}
@@ -210,37 +199,38 @@ riemann_solve_for_middle_state_flux(const float *WL, const float *WR,
/* Handle vacuum: vacuum does not require iteration and is always exact */
if (riemann_is_vacuum(WL, WR, uL, uR, aL, aR)) {
- totflux[0] = 0.f;
- totflux[1] = 0.f;
- totflux[2] = 0.f;
- totflux[3] = 0.f;
- totflux[4] = 0.f;
+ totflux[0] = 0.0f;
+ totflux[1] = 0.0f;
+ totflux[2] = 0.0f;
+ totflux[3] = 0.0f;
+ totflux[4] = 0.0f;
return;
}
/* STEP 1: pressure estimate */
- const float rhobar = 0.5f * (WL[0] + WR[0]);
- const float abar = 0.5f * (aL + aR);
- const float pPVRS = 0.5f * (WL[4] + WR[4]) - 0.5f * (uR - uL) * rhobar * abar;
+ const float rhobar = WL[0] + WR[0];
+ const float abar = aL + aR;
+ const float pPVRS =
+ 0.5f * ((WL[4] + WR[4]) - 0.25f * (uR - uL) * rhobar * abar);
const float pstar = max(0.f, pPVRS);
/* STEP 2: wave speed estimates
all these speeds are along the interface normal, since uL and uR are */
- float qL = 1.f;
- if (pstar > WL[4] && WL[4] > 0.f) {
- qL = sqrtf(1.f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
- (pstar / WL[4] - 1.f));
+ float qL = 1.0f;
+ if (pstar > WL[4] && WL[4] > 0.0f) {
+ qL = sqrtf(1.0f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
+ (pstar / WL[4] - 1.0f));
}
- float qR = 1.f;
- if (pstar > WR[4] && WR[4] > 0.f) {
- qR = sqrtf(1.f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
- (pstar / WR[4] - 1.f));
+ float qR = 1.0f;
+ if (pstar > WR[4] && WR[4] > 0.0f) {
+ qR = sqrtf(1.0f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
+ (pstar / WR[4] - 1.0f));
}
- const float SL = uL - aL * qL;
- const float SR = uR + aR * qR;
+ const float SLmuL = -aL * qL;
+ const float SRmuR = aR * qR;
const float Sstar =
- (WR[4] - WL[4] + WL[0] * uL * (SL - uL) - WR[0] * uR * (SR - uR)) /
- (WL[0] * (SL - uL) - WR[0] * (SR - uR));
+ (WR[4] - WL[4] + WL[0] * uL * SLmuL - WR[0] * uR * SRmuR) /
+ (WL[0] * SLmuL - WR[0] * SRmuR);
totflux[0] = 0.0f;
totflux[1] = pstar * n[0];