diff --git a/src/Makefile.am b/src/Makefile.am
index 809f5c34642392d001f1987c1c8926b8b97eae1e..2ddcdb0908201c65053d7cc5380a4217277b5c13 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -86,6 +86,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
hydro/Gizmo/hydro_slope_limiters_cell.h \
hydro/Gizmo/hydro_slope_limiters_face.h \
hydro/Gizmo/hydro_slope_limiters.h \
+ hydro/Gizmo/hydro_unphysical.h \
+ hydro/Gizmo/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/const.h b/src/const.h
index 6962ee8bca32e92664e3f20cdb23e7cb6fbc4abd..144049d6b6897235067d840a150f9a0a49411aed 100644
--- a/src/const.h
+++ b/src/const.h
@@ -52,8 +52,43 @@
/* Options to control the movement of particles for GIZMO_SPH. */
/* This option disables particle movement */
//#define GIZMO_FIX_PARTICLES
+/* Try to keep cells regular by adding a correction velocity. */
+#define GIZMO_STEER_MOTION
//#define GIZMO_TOTAL_ENERGY
+/* Options to control handling of unphysical values (GIZMO_SPH only). */
+/* In GIZMO, mass and energy (and hence density and pressure) can in principle
+ become negative, which will cause unwanted behaviour that can make the code
+ crash.
+ If no options are selected below, we assume (and pray) that this will not
+ happen, and add no restrictions to how these variables are treated. */
+/* Check for unphysical values and crash if they occur. */
+//#define GIZMO_UNPHYSICAL_ERROR
+/* Check for unphysical values and reset them to safe values. */
+#define GIZMO_UNPHYSICAL_RESCUE
+/* Show a warning message if an unphysical value was reset (only works if
+ GIZMO_UNPHYSICAL_RESCUE is also selected). */
+#define GIZMO_UNPHYSICAL_WARNING
+
+/* Parameters that control how GIZMO handles pathological particle
+ configurations. */
+/* Show a warning message if a pathological configuration has been detected. */
+#define GIZMO_PATHOLOGICAL_WARNING
+/* Crash if a pathological configuration has been detected. */
+//#define GIZMO_PATHOLOGICAL_ERROR
+/* Maximum allowed gradient matrix condition number. If the condition number of
+ the gradient matrix (defined in equation C1 in Hopkins, 2015) is larger than
+ this value, we artificially increase the number of neighbours to get a more
+ homogeneous sampling. */
+#define const_gizmo_max_condition_number 100.0f
+/* Correction factor applied to the particle wcount to force more neighbours if
+ the condition number is too large. */
+#define const_gizmo_w_correction_factor 0.9f
+/* Lower limit on the wcount correction factor. If the condition number is still
+ too high after this wcount correction has been applied, we give up on the
+ gradient matrix and use SPH gradients instead. */
+#define const_gizmo_min_wcorr 0.5f
+
/* Types of gradients to use for SHADOWFAX_SPH */
/* If no option is chosen, no gradients are used (first order scheme) */
#define SHADOWFAX_GRADIENTS
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 2e340a03b99ae51bc49a2e57456f4d6838d62f21..6d39c54d2ddc3571ac34c54fc9eede6f7dee6ac5 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -2,6 +2,7 @@
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2016, 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
@@ -24,9 +25,13 @@
#include "equation_of_state.h"
#include "hydro_gradients.h"
#include "hydro_space.h"
+#include "hydro_unphysical.h"
+#include "hydro_velocities.h"
#include "minmax.h"
#include "riemann.h"
+//#define GIZMO_LLOYD_ITERATION
+
/**
* @brief Computes the hydro time-step of a given particle
*
@@ -40,6 +45,10 @@ __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
+
if (p->timestepvars.vmax == 0.) {
/* vmax can be zero in vacuum cells that only have vacuum neighbours */
/* in this case, the time step should be limited by the maximally
@@ -47,7 +56,9 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
the time step to a very large value */
return FLT_MAX;
} else {
- return CFL_condition * p->h / fabsf(p->timestepvars.vmax);
+ const float psize = powf(p->geometry.volume / hydro_dimension_unit_sphere,
+ hydro_dimension_inv);
+ return 2. * CFL_condition * psize / fabsf(p->timestepvars.vmax);
}
}
@@ -128,16 +139,27 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->conserved.momentum[2] * p->primitives.v[2]);
#endif
-#if defined(GIZMO_FIX_PARTICLES)
- /* make sure the particles are initially at rest */
+#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.;
#endif
- xp->v_full[0] = p->v[0];
- xp->v_full[1] = p->v[1];
- xp->v_full[2] = p->v[2];
+ /* initialize the particle velocity based on the primitive fluid velocity */
+ hydro_velocities_init(p, xp);
/* we cannot initialize wcorr in init_part, as init_part gets called every
time the density loop is repeated, and the whole point of storing wcorr
@@ -169,6 +191,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->geometry.matrix_E[2][0] = 0.0f;
p->geometry.matrix_E[2][1] = 0.0f;
p->geometry.matrix_E[2][2] = 0.0f;
+ p->geometry.centroid[0] = 0.0f;
+ p->geometry.centroid[1] = 0.0f;
+ p->geometry.centroid[2] = 0.0f;
p->geometry.Atot = 0.0f;
/* Set the active flag to active. */
@@ -226,6 +251,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->geometry.matrix_E[2][1] = ihdim * p->geometry.matrix_E[2][1];
p->geometry.matrix_E[2][2] = ihdim * p->geometry.matrix_E[2][2];
+ p->geometry.centroid[0] *= kernel_norm;
+ 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;
+
/* Check the condition number to see if we have a stable geometry. */
float condition_number_E = 0.0f;
int i, j;
@@ -249,12 +282,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
float condition_number =
hydro_dimension_inv * sqrtf(condition_number_E * condition_number_Einv);
- if (condition_number > 100.0f) {
- // error("Condition number larger than 100!");
- // message("Condition number too large: %g (p->id: %llu)!",
- // condition_number, p->id);
+ if (condition_number > const_gizmo_max_condition_number &&
+ p->density.wcorr > const_gizmo_min_wcorr) {
+#ifdef GIZMO_PATHOLOGICAL_ERROR
+ error("Condition number larger than %g (%g)!",
+ const_gizmo_max_condition_number, condition_number);
+#endif
+#ifdef GIZMO_PATHOLOGICAL_WARNING
+ message("Condition number too large: %g (> %g, p->id: %llu)!",
+ condition_number, const_gizmo_max_condition_number, p->id);
+#endif
/* add a correction to the number of neighbours for this particle */
- p->density.wcorr *= 0.75;
+ p->density.wcorr *= const_gizmo_w_correction_factor;
}
hydro_gradients_init(p);
@@ -264,8 +303,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
const float m = p->conserved.mass;
#ifdef SWIFT_DEBUG_CHECKS
- if (m == 0.) {
- error("Mass is 0!");
+ if (m < 0.) {
+ error("Mass is negative!");
}
if (volume == 0.) {
@@ -278,15 +317,20 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
momentum[1] = p->conserved.momentum[1];
momentum[2] = p->conserved.momentum[2];
p->primitives.rho = m / volume;
- p->primitives.v[0] = momentum[0] / m;
- p->primitives.v[1] = momentum[1] / m;
- p->primitives.v[2] = momentum[2] / m;
+ if (m == 0.) {
+ p->primitives.v[0] = 0.;
+ p->primitives.v[1] = 0.;
+ p->primitives.v[2] = 0.;
+ } else {
+ p->primitives.v[0] = momentum[0] / m;
+ p->primitives.v[1] = momentum[1] / m;
+ p->primitives.v[2] = momentum[2] / m;
+ }
#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 = const_isothermal_soundspeed * const_isothermal_soundspeed *
- p->primitives.rho;
+ p->primitives.P = gas_pressure_from_internal_energy(p->primitives.rho, 0.);
#else
float energy = p->conserved.energy;
@@ -304,12 +348,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
#endif
/* sanity checks */
- /* it would probably be safer to throw a warning if netive densities or
- pressures occur */
- if (p->primitives.rho < 0.0f || p->primitives.P < 0.0f) {
- p->primitives.rho = 0.0f;
- p->primitives.P = 0.0f;
- }
+ gizmo_check_physical_quantity("density", p->primitives.rho);
+ gizmo_check_physical_quantity("pressure", 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
/* Add a correction factor to wcount (to force a neighbour number increase if
the geometry matrix is close to singular) */
@@ -330,8 +379,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
*
* @param p The particle to act upon.
* @param xp The extended particle data to act upon.
- * @param ti_current Current integer time.
- * @param timeBase Conversion factor between integer time and physical time.
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* restrict p, struct xpart* restrict xp) {
@@ -340,10 +387,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->timestepvars.vmax = 0.0f;
/* Set the actual velocity of the particle */
- /* if GIZMO_FIX_PARTICLES has been selected, v_full will always be zero */
- p->force.v_full[0] = xp->v_full[0];
- p->force.v_full[1] = xp->v_full[1];
- p->force.v_full[2] = xp->v_full[2];
+ hydro_velocities_prepare_force(p, xp);
}
/**
@@ -364,6 +408,11 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
p->gravity.mflux[0] = 0.0f;
p->gravity.mflux[1] = 0.0f;
p->gravity.mflux[2] = 0.0f;
+
+#ifdef GIZMO_LLOYD_ITERATION
+ /* reset the gradients to zero, as we don't want them */
+ hydro_gradients_init(p);
+#endif
}
/**
@@ -422,6 +471,10 @@ __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) {
+#ifdef GIZMO_LLOYD_ITERATION
+ return;
+#endif
+
const float h_inv = 1.0f / p->h;
/* Predict smoothing length */
@@ -432,8 +485,9 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
else
h_corr = expf(w1);
- /* Limit the smoothing length correction. */
- if (h_corr < 2.0f) {
+ /* Limit the smoothing length correction (and make sure it is always
+ positive). */
+ if (h_corr < 2.0f && h_corr > 0.) {
p->h *= h_corr;
}
@@ -483,22 +537,13 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
/* set the variables that are used to drift the primitive variables */
- /* Add normalization to h_dt. */
- p->force.h_dt *= p->h * hydro_dimension_inv;
-
- if (p->force.dt) {
+ if (p->force.dt > 0.) {
p->du_dt = p->conserved.flux.energy / p->force.dt;
} else {
p->du_dt = 0.0f;
}
-#if defined(GIZMO_FIX_PARTICLES)
- p->du_dt = 0.0f;
-
- /* 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;
-#endif
+ hydro_velocities_end_force(p);
}
/**
@@ -527,7 +572,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.energy += p->conserved.flux.energy;
#endif
+ gizmo_check_physical_quantity("mass", p->conserved.mass);
+ gizmo_check_physical_quantity("energy", p->conserved.energy);
+
#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)!",
@@ -535,7 +585,10 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
}
if (p->conserved.energy < 0.) {
- error("Negative energy after conserved variables update!");
+ error(
+ "Negative energy after conserved variables update (energy: %g, "
+ "denergy: %g)!",
+ p->conserved.energy, p->conserved.flux.energy);
}
#endif
@@ -549,7 +602,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
a_grav[2] = p->gpart->a_grav[2];
/* Store the gravitational acceleration for later use. */
- /* This is currently only used for output purposes. */
+ /* This is used for the prediction step. */
p->gravity.old_a[0] = a_grav[0];
p->gravity.old_a[1] = a_grav[1];
p->gravity.old_a[2] = a_grav[2];
@@ -564,7 +617,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
-#if !defined(EOS_ISOTHERMAL_GAS) && defined(GIZMO_TOTAL_ENERGY)
+#if !defined(EOS_ISOTHERMAL_GAS)
/* This part still needs to be tested! */
p->conserved.energy += dt * (p->conserved.momentum[0] * a_grav[0] +
p->conserved.momentum[1] * a_grav[1] +
@@ -585,45 +638,25 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.flux.momentum[2] = 0.0f;
p->conserved.flux.energy = 0.0f;
-#if defined(GIZMO_FIX_PARTICLES)
- xp->v_full[0] = 0.;
- xp->v_full[1] = 0.;
- xp->v_full[2] = 0.;
-
- p->v[0] = 0.;
- p->v[1] = 0.;
- p->v[2] = 0.;
-
- if (p->gpart) {
- p->gpart->v_full[0] = 0.;
- p->gpart->v_full[1] = 0.;
- p->gpart->v_full[2] = 0.;
- }
-#else
- /* Set particle movement */
- if (p->conserved.mass > 0.) {
- xp->v_full[0] = p->conserved.momentum[0] / p->conserved.mass;
- xp->v_full[1] = p->conserved.momentum[1] / p->conserved.mass;
- xp->v_full[2] = p->conserved.momentum[2] / p->conserved.mass;
- } else {
- /* vacuum particles don't move */
- xp->v_full[0] = 0.;
- xp->v_full[1] = 0.;
- xp->v_full[2] = 0.;
- }
+ hydro_velocities_set(p, xp);
+
+#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];
-
- /* Update gpart! */
- /* This is essential, as the gpart drift is done independently from the part
- drift, and we don't want the gpart and the part to have different
- positions! */
- 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];
- }
#endif
/* reset wcorr */
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index a5c1e9038d0d3de6896afe773e3193a2304a6b6b..5ad6d87619a7629a703a8b9c03d089e69ffbdf7d 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -22,6 +22,7 @@
#define SWIFT_HYDRO_GRADIENTS_H
#include "hydro_slope_limiters.h"
+#include "hydro_unphysical.h"
#include "riemann.h"
#if defined(GRADIENTS_SPH)
@@ -98,6 +99,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
float xij_j[3];
int k;
float xfac;
+ float a_grav_i[3], a_grav_j[3];
/* perform gradient reconstruction in space and time */
/* space */
@@ -139,37 +141,38 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
pj->primitives.gradients.P[1] * xij_j[1] +
pj->primitives.gradients.P[2] * xij_j[2];
+ a_grav_i[0] = pi->gravity.old_a[0];
+ a_grav_i[1] = pi->gravity.old_a[1];
+ a_grav_i[2] = pi->gravity.old_a[2];
+
+ a_grav_i[0] += pi->gravity.grad_a[0][0] * xij_i[0] +
+ pi->gravity.grad_a[0][1] * xij_i[1] +
+ pi->gravity.grad_a[0][2] * xij_i[2];
+ a_grav_i[1] += pi->gravity.grad_a[1][0] * xij_i[0] +
+ pi->gravity.grad_a[1][1] * xij_i[1] +
+ pi->gravity.grad_a[1][2] * xij_i[2];
+ a_grav_i[2] += pi->gravity.grad_a[2][0] * xij_i[0] +
+ pi->gravity.grad_a[2][1] * xij_i[1] +
+ pi->gravity.grad_a[2][2] * xij_i[2];
+
+ a_grav_j[0] = pj->gravity.old_a[0];
+ a_grav_j[1] = pj->gravity.old_a[1];
+ a_grav_j[2] = pj->gravity.old_a[2];
+
+ a_grav_j[0] += pj->gravity.grad_a[0][0] * xij_j[0] +
+ pj->gravity.grad_a[0][1] * xij_j[1] +
+ pj->gravity.grad_a[0][2] * xij_j[2];
+ a_grav_j[1] += pj->gravity.grad_a[1][0] * xij_j[0] +
+ pj->gravity.grad_a[1][1] * xij_j[1] +
+ pj->gravity.grad_a[1][2] * xij_j[2];
+ a_grav_j[2] += pj->gravity.grad_a[2][0] * xij_j[0] +
+ pj->gravity.grad_a[2][1] * xij_j[1] +
+ pj->gravity.grad_a[2][2] * xij_j[2];
+
hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
/* time */
if (Wi[0] > 0.0f) {
-#ifdef EOS_ISOTHERMAL_GAS
- dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
- Wi[2] * pi->primitives.gradients.rho[1] +
- Wi[3] * pi->primitives.gradients.rho[2] +
- Wi[0] * (pi->primitives.gradients.v[0][0] +
- pi->primitives.gradients.v[1][1] +
- pi->primitives.gradients.v[2][2]));
- dWi[1] -= 0.5 * mindt *
- (Wi[1] * pi->primitives.gradients.v[0][0] +
- Wi[2] * pi->primitives.gradients.v[0][1] +
- Wi[3] * pi->primitives.gradients.v[0][2] +
- const_isothermal_soundspeed * const_isothermal_soundspeed *
- pi->primitives.gradients.rho[0] / Wi[0]);
- dWi[2] -= 0.5 * mindt *
- (Wi[1] * pi->primitives.gradients.v[1][0] +
- Wi[2] * pi->primitives.gradients.v[1][1] +
- Wi[3] * pi->primitives.gradients.v[1][2] +
- const_isothermal_soundspeed * const_isothermal_soundspeed *
- pi->primitives.gradients.rho[1] / Wi[0]);
- dWi[3] -= 0.5 * mindt *
- (Wi[1] * pi->primitives.gradients.v[2][0] +
- Wi[2] * pi->primitives.gradients.v[2][1] +
- Wi[3] * pi->primitives.gradients.v[2][2] +
- const_isothermal_soundspeed * const_isothermal_soundspeed *
- pi->primitives.gradients.rho[2] / Wi[0]);
-/* we don't care about P in this case */
-#else
dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
Wi[2] * pi->primitives.gradients.rho[1] +
Wi[3] * pi->primitives.gradients.rho[2] +
@@ -195,36 +198,13 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
pi->primitives.gradients.v[1][1] +
pi->primitives.gradients.v[2][2]));
-#endif
+
+ dWi[1] += 0.5 * mindt * a_grav_i[0];
+ dWi[2] += 0.5 * mindt * a_grav_i[1];
+ dWi[3] += 0.5 * mindt * a_grav_i[2];
}
if (Wj[0] > 0.0f) {
-#ifdef EOS_ISOTHERMAL_GAS
- dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
- Wj[2] * pj->primitives.gradients.rho[1] +
- Wj[3] * pj->primitives.gradients.rho[2] +
- Wj[0] * (pj->primitives.gradients.v[0][0] +
- pj->primitives.gradients.v[1][1] +
- pj->primitives.gradients.v[2][2]));
- dWj[1] -= 0.5 * mindt *
- (Wj[1] * pj->primitives.gradients.v[0][0] +
- Wj[2] * pj->primitives.gradients.v[0][1] +
- Wj[3] * pj->primitives.gradients.v[0][2] +
- const_isothermal_soundspeed * const_isothermal_soundspeed *
- pj->primitives.gradients.rho[0] / Wj[0]);
- dWj[2] -= 0.5 * mindt *
- (Wj[1] * pj->primitives.gradients.v[1][0] +
- Wj[2] * pj->primitives.gradients.v[1][1] +
- Wj[3] * pj->primitives.gradients.v[1][2] +
- const_isothermal_soundspeed * const_isothermal_soundspeed *
- pj->primitives.gradients.rho[1] / Wj[0]);
- dWj[3] -= 0.5 * mindt *
- (Wj[1] * pj->primitives.gradients.v[2][0] +
- Wj[2] * pj->primitives.gradients.v[2][1] +
- Wj[3] * pj->primitives.gradients.v[2][2] +
- const_isothermal_soundspeed * const_isothermal_soundspeed *
- pj->primitives.gradients.rho[2] / Wj[0]);
-#else
dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
Wj[2] * pj->primitives.gradients.rho[1] +
Wj[3] * pj->primitives.gradients.rho[2] +
@@ -250,36 +230,28 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
pj->primitives.gradients.v[1][1] +
pj->primitives.gradients.v[2][2]));
-#endif
- }
- if (-dWi[0] > Wi[0]) {
- Wi[0] = 0.0f;
- } else {
- Wi[0] += dWi[0];
+ dWj[1] += 0.5 * mindt * a_grav_j[0];
+ dWj[2] += 0.5 * mindt * a_grav_j[1];
+ dWj[3] += 0.5 * mindt * a_grav_j[2];
}
+
+ Wi[0] += dWi[0];
Wi[1] += dWi[1];
Wi[2] += dWi[2];
Wi[3] += dWi[3];
- if (-dWi[4] > Wi[4]) {
- Wi[4] = 0.0f;
- } else {
- Wi[4] += dWi[4];
- }
+ Wi[4] += dWi[4];
- if (-dWj[0] > Wj[0]) {
- Wj[0] = 0.0f;
- } else {
- Wj[0] += dWj[0];
- }
+ Wj[0] += dWj[0];
Wj[1] += dWj[1];
Wj[2] += dWj[2];
Wj[3] += dWj[3];
- if (-dWj[4] > Wj[4]) {
- Wj[4] = 0.0f;
- } else {
- Wj[4] += dWj[4];
- }
+ Wj[4] += dWj[4];
+
+ gizmo_check_physical_quantity("density", Wi[0]);
+ gizmo_check_physical_quantity("pressure", Wi[4]);
+ gizmo_check_physical_quantity("density", Wj[0]);
+ gizmo_check_physical_quantity("pressure", Wj[4]);
}
#endif // SWIFT_HYDRO_GRADIENTS_H
diff --git a/src/hydro/Gizmo/hydro_gradients_gizmo.h b/src/hydro/Gizmo/hydro_gradients_gizmo.h
index aa6e4406b94e7a5cafcd0ca556162476003477de..ee3ad6919f81f042ceacc5db8b4e818d63c90266 100644
--- a/src/hydro/Gizmo/hydro_gradients_gizmo.h
+++ b/src/hydro/Gizmo/hydro_gradients_gizmo.h
@@ -45,6 +45,18 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_init(
p->primitives.gradients.P[1] = 0.0f;
p->primitives.gradients.P[2] = 0.0f;
+ p->gravity.grad_a[0][0] = 0.0f;
+ p->gravity.grad_a[0][1] = 0.0f;
+ p->gravity.grad_a[0][2] = 0.0f;
+
+ p->gravity.grad_a[1][0] = 0.0f;
+ p->gravity.grad_a[1][1] = 0.0f;
+ p->gravity.grad_a[1][2] = 0.0f;
+
+ p->gravity.grad_a[2][0] = 0.0f;
+ p->gravity.grad_a[2][1] = 0.0f;
+ p->gravity.grad_a[2][2] = 0.0f;
+
hydro_slope_limit_cell_init(p);
}
@@ -93,56 +105,146 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
- /* 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]);
+ 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]);
+
+ pi->gravity.grad_a[0][0] +=
+ (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->gravity.grad_a[0][1] +=
+ (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->gravity.grad_a[0][2] +=
+ (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->gravity.grad_a[1][0] +=
+ (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->gravity.grad_a[1][1] +=
+ (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->gravity.grad_a[1][2] +=
+ (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->gravity.grad_a[2][0] +=
+ (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->gravity.grad_a[2][1] +=
+ (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->gravity.grad_a[2][2] +=
+ (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ } else {
+ /* The gradient matrix was not well-behaved, switch to SPH gradients */
+
+ pi->primitives.gradients.rho[0] -=
+ wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pi->primitives.gradients.rho[1] -=
+ wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pi->primitives.gradients.rho[2] -=
+ wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
+
+ pi->primitives.gradients.v[0][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pi->primitives.gradients.v[0][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pi->primitives.gradients.v[0][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+
+ pi->primitives.gradients.v[1][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pi->primitives.gradients.v[1][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pi->primitives.gradients.v[1][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+
+ pi->primitives.gradients.v[2][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pi->primitives.gradients.v[2][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pi->primitives.gradients.v[2][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+
+ pi->primitives.gradients.P[0] -=
+ wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
+ pi->primitives.gradients.P[1] -=
+ wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
+ pi->primitives.gradients.P[2] -=
+ wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
+
+ pi->gravity.grad_a[0][0] -=
+ wi_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+ pi->gravity.grad_a[0][1] -=
+ wi_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+ pi->gravity.grad_a[0][2] -=
+ wi_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+
+ pi->gravity.grad_a[1][0] -=
+ wi_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+ pi->gravity.grad_a[1][1] -=
+ wi_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+ pi->gravity.grad_a[1][2] -=
+ wi_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+
+ pi->gravity.grad_a[2][0] -=
+ wi_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+ pi->gravity.grad_a[2][1] -=
+ wi_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+ pi->gravity.grad_a[2][2] -=
+ wi_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+ }
hydro_slope_limit_cell_collect(pi, pj, r);
@@ -151,57 +253,146 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
- /* 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]);
+ 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]);
+
+ pj->gravity.grad_a[0][0] +=
+ (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->gravity.grad_a[0][1] +=
+ (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->gravity.grad_a[0][2] +=
+ (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+ pj->gravity.grad_a[1][0] +=
+ (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->gravity.grad_a[1][1] +=
+ (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->gravity.grad_a[1][2] +=
+ (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+ pj->gravity.grad_a[2][0] +=
+ (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->gravity.grad_a[2][1] +=
+ (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->gravity.grad_a[2][2] +=
+ (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+ } else {
+ /* SPH gradients */
+
+ pj->primitives.gradients.rho[0] -=
+ wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pj->primitives.gradients.rho[1] -=
+ wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pj->primitives.gradients.rho[2] -=
+ wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
+
+ pj->primitives.gradients.v[0][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pj->primitives.gradients.v[0][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pj->primitives.gradients.v[0][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+
+ pj->primitives.gradients.v[1][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pj->primitives.gradients.v[1][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pj->primitives.gradients.v[1][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pj->primitives.gradients.v[2][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pj->primitives.gradients.v[2][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pj->primitives.gradients.v[2][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+
+ pj->primitives.gradients.P[0] -=
+ wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
+ pj->primitives.gradients.P[1] -=
+ wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
+ pj->primitives.gradients.P[2] -=
+ wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
+
+ pj->gravity.grad_a[0][0] -=
+ wj_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+ pj->gravity.grad_a[0][1] -=
+ wj_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+ pj->gravity.grad_a[0][2] -=
+ wj_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+
+ pj->gravity.grad_a[1][0] -=
+ wj_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+ pj->gravity.grad_a[1][1] -=
+ wj_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+ pj->gravity.grad_a[1][2] -=
+ wj_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+
+ pj->gravity.grad_a[2][0] -=
+ wj_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+ pj->gravity.grad_a[2][1] -=
+ wj_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+ pj->gravity.grad_a[2][2] -=
+ wj_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+ }
hydro_slope_limit_cell_collect(pj, pi, r);
}
@@ -250,56 +441,145 @@ hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
- /* 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]);
+ 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]);
+
+ pi->gravity.grad_a[0][0] +=
+ (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->gravity.grad_a[0][1] +=
+ (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->gravity.grad_a[0][2] +=
+ (pi->gravity.old_a[0] - pj->gravity.old_a[0]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->gravity.grad_a[1][0] +=
+ (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->gravity.grad_a[1][1] +=
+ (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->gravity.grad_a[1][2] +=
+ (pi->gravity.old_a[1] - pj->gravity.old_a[1]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->gravity.grad_a[2][0] +=
+ (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->gravity.grad_a[2][1] +=
+ (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->gravity.grad_a[2][2] +=
+ (pi->gravity.old_a[2] - pj->gravity.old_a[2]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ } else {
+ /* Gradient matrix is not well-behaved, switch to SPH gradients */
+
+ pi->primitives.gradients.rho[0] -=
+ wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pi->primitives.gradients.rho[1] -=
+ wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) / r;
+ pi->primitives.gradients.rho[2] -=
+ wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) / r;
+
+ pi->primitives.gradients.v[0][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pi->primitives.gradients.v[0][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pi->primitives.gradients.v[0][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) / r;
+ pi->primitives.gradients.v[1][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pi->primitives.gradients.v[1][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+ pi->primitives.gradients.v[1][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) / r;
+
+ pi->primitives.gradients.v[2][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pi->primitives.gradients.v[2][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+ pi->primitives.gradients.v[2][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) / r;
+
+ pi->primitives.gradients.P[0] -=
+ wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) / r;
+ pi->primitives.gradients.P[1] -=
+ wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) / r;
+ pi->primitives.gradients.P[2] -=
+ wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) / r;
+
+ pi->gravity.grad_a[0][0] -=
+ wi_dx * dx[0] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+ pi->gravity.grad_a[0][1] -=
+ wi_dx * dx[1] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+ pi->gravity.grad_a[0][2] -=
+ wi_dx * dx[2] * (pi->gravity.old_a[0] - pj->gravity.old_a[0]) / r;
+
+ pi->gravity.grad_a[1][0] -=
+ wi_dx * dx[0] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+ pi->gravity.grad_a[1][1] -=
+ wi_dx * dx[1] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+ pi->gravity.grad_a[1][2] -=
+ wi_dx * dx[2] * (pi->gravity.old_a[1] - pj->gravity.old_a[1]) / r;
+
+ pi->gravity.grad_a[2][0] -=
+ wi_dx * dx[0] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+ pi->gravity.grad_a[2][1] -=
+ wi_dx * dx[1] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+ pi->gravity.grad_a[2][2] -=
+ wi_dx * dx[2] * (pi->gravity.old_a[2] - pj->gravity.old_a[2]) / r;
+ }
hydro_slope_limit_cell_collect(pi, pj, r);
}
@@ -319,23 +599,73 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
ih = 1.0f / h;
const float ihdim = pow_dimension(ih);
- 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;
+ 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;
+
+ p->gravity.grad_a[0][0] *= ihdim;
+ p->gravity.grad_a[0][1] *= ihdim;
+ p->gravity.grad_a[0][2] *= ihdim;
+
+ p->gravity.grad_a[1][0] *= ihdim;
+ p->gravity.grad_a[1][1] *= ihdim;
+ p->gravity.grad_a[1][2] *= ihdim;
+
+ p->gravity.grad_a[2][0] *= ihdim;
+ p->gravity.grad_a[2][1] *= ihdim;
+ p->gravity.grad_a[2][2] *= ihdim;
+ } else {
+ const float ihdimp1 = pow_dimension_plus_one(ih);
+
+ float volume = p->geometry.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->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;
+
+ p->gravity.grad_a[0][0] *= ihdimp1 * volume;
+ p->gravity.grad_a[0][1] *= ihdimp1 * volume;
+ p->gravity.grad_a[0][2] *= ihdimp1 * volume;
+
+ p->gravity.grad_a[1][0] *= ihdimp1 * volume;
+ p->gravity.grad_a[1][1] *= ihdimp1 * volume;
+ p->gravity.grad_a[1][2] *= ihdimp1 * volume;
+
+ p->gravity.grad_a[2][0] *= ihdimp1 * volume;
+ p->gravity.grad_a[2][1] *= ihdimp1 * volume;
+ p->gravity.grad_a[2][2] *= ihdimp1 * volume;
+ }
hydro_slope_limit_cell(p);
}
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index d707e0ee1b5707086393ea206ea9f0f60f9c1853..8798dc859a790a83ab7a3b6f1709b1302f574581 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -23,6 +23,8 @@
#include "hydro_gradients.h"
#include "riemann.h"
+#define GIZMO_VOLUME_CORRECTION
+
/**
* @brief Calculate the volume interaction between particle i and particle j
*
@@ -62,6 +64,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
for (k = 0; k < 3; k++)
for (l = 0; l < 3; l++) pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
+ pi->geometry.centroid[0] -= dx[0] * wi;
+ pi->geometry.centroid[1] -= dx[1] * wi;
+ pi->geometry.centroid[2] -= dx[2] * wi;
+
/* Compute density of pj. */
h_inv = 1.0 / hj;
xj = r * h_inv;
@@ -74,6 +80,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pj->geometry.volume += wj;
for (k = 0; k < 3; k++)
for (l = 0; l < 3; l++) pj->geometry.matrix_E[k][l] += dx[k] * dx[l] * wj;
+
+ pj->geometry.centroid[0] += dx[0] * wj;
+ pj->geometry.centroid[1] += dx[1] * wj;
+ pj->geometry.centroid[2] += dx[2] * wj;
}
/**
@@ -117,6 +127,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
pi->geometry.volume += wi;
for (k = 0; k < 3; k++)
for (l = 0; l < 3; l++) pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
+
+ pi->geometry.centroid[0] -= dx[0] * wi;
+ pi->geometry.centroid[1] -= dx[1] * wi;
+ pi->geometry.centroid[2] -= dx[2] * wi;
}
/**
@@ -325,14 +339,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* calculate the maximal signal velocity */
if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
-#ifdef EOS_ISOTHERMAL_GAS
- /* we use a value that is slightly higher than necessary, since the correct
- value does not always work */
- vmax = 2.5 * const_isothermal_soundspeed;
-#else
vmax =
sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
-#endif
} else {
vmax = 0.0f;
}
@@ -381,23 +389,63 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* Compute area */
/* eqn. (7) */
Anorm = 0.0f;
- for (k = 0; k < 3; k++) {
- /* we add a minus sign since dx is pi->x - pj->x */
- A[k] = -Vi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) * wi *
- hi_inv_dim -
- Vj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) * wj *
- hj_inv_dim;
- Anorm += A[k] * A[k];
+ if (pi->density.wcorr > const_gizmo_min_wcorr &&
+ pj->density.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
+ different
+ smoothing lengths), then the expressions below are more stable. */
+ float Xi = Vi;
+ float Xj = Vj;
+#ifdef GIZMO_VOLUME_CORRECTION
+ if (fabsf(Vi - Vj) / fminf(Vi, Vj) > 1.5 * hydro_dimension) {
+ Xi = (Vi * hj + Vj * hi) / (hi + hj);
+ Xj = Xi;
+ }
+#endif
+ for (k = 0; k < 3; k++) {
+ /* we add a minus sign since dx is pi->x - pj->x */
+ A[k] = -Xi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) *
+ wj * hj_inv_dim -
+ Xj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) *
+ wi * hi_inv_dim;
+ Anorm += A[k] * A[k];
+ }
+ } else {
+ /* ill condition gradient matrix: revert to SPH face area */
+ Anorm = -(hidp1 * Vi * Vi * wi_dx + hjdp1 * Vj * Vj * wj_dx) * ri;
+ A[0] = -Anorm * dx[0];
+ A[1] = -Anorm * dx[1];
+ A[2] = -Anorm * dx[2];
+ Anorm *= Anorm * r2;
}
- if (!Anorm) {
+ if (Anorm == 0.) {
/* if the interface has no area, nothing happens and we return */
/* continuing results in dividing by zero and NaN's... */
return;
}
- /* compute the normal vector of the interface */
Anorm = sqrtf(Anorm);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* For stability reasons, we do require A and dx to have opposite
+ directions (basically meaning that the surface normal for the surface
+ always points from particle i to particle j, as it would in a real
+ moving-mesh code). If not, our scheme is no longer upwind and hence can
+ become unstable. */
+ float dA_dot_dx = A[0] * dx[0] + A[1] * dx[1] + A[2] * dx[2];
+ /* In GIZMO, Phil Hopkins reverts to an SPH integration scheme if this
+ happens. We curently just ignore this case and display a message. */
+ const float rdim = pow_dimension(r);
+ if (dA_dot_dx > 1.e-6 * rdim) {
+ message("Ill conditioned gradient matrix (%g %g %g %g %g)!", dA_dot_dx,
+ Anorm, Vi, Vj, r);
+ }
+#endif
+
+ /* compute the normal vector of the interface */
for (k = 0; k < 3; k++) n_unit[k] = A[k] / Anorm;
/* Compute interface position (relative to pi, since we don't need the actual
@@ -436,43 +484,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* we don't need to rotate, we can use the unit vector in the Riemann problem
* itself (see GIZMO) */
- if (Wi[0] < 0.0f || Wj[0] < 0.0f || Wi[4] < 0.0f || Wj[4] < 0.0f) {
- printf("mindt: %g\n", mindt);
- printf("WL: %g %g %g %g %g\n", pi->primitives.rho, pi->primitives.v[0],
- pi->primitives.v[1], pi->primitives.v[2], pi->primitives.P);
-#ifdef USE_GRADIENTS
- printf("dWL: %g %g %g %g %g\n", dWi[0], dWi[1], dWi[2], dWi[3], dWi[4]);
-#endif
- printf("gradWL[0]: %g %g %g\n", pi->primitives.gradients.rho[0],
- pi->primitives.gradients.rho[1], pi->primitives.gradients.rho[2]);
- printf("gradWL[1]: %g %g %g\n", pi->primitives.gradients.v[0][0],
- pi->primitives.gradients.v[0][1], pi->primitives.gradients.v[0][2]);
- printf("gradWL[2]: %g %g %g\n", pi->primitives.gradients.v[1][0],
- pi->primitives.gradients.v[1][1], pi->primitives.gradients.v[1][2]);
- printf("gradWL[3]: %g %g %g\n", pi->primitives.gradients.v[2][0],
- pi->primitives.gradients.v[2][1], pi->primitives.gradients.v[2][2]);
- printf("gradWL[4]: %g %g %g\n", pi->primitives.gradients.P[0],
- pi->primitives.gradients.P[1], pi->primitives.gradients.P[2]);
- printf("WL': %g %g %g %g %g\n", Wi[0], Wi[1], Wi[2], Wi[3], Wi[4]);
- printf("WR: %g %g %g %g %g\n", pj->primitives.rho, pj->primitives.v[0],
- pj->primitives.v[1], pj->primitives.v[2], pj->primitives.P);
-#ifdef USE_GRADIENTS
- printf("dWR: %g %g %g %g %g\n", dWj[0], dWj[1], dWj[2], dWj[3], dWj[4]);
-#endif
- printf("gradWR[0]: %g %g %g\n", pj->primitives.gradients.rho[0],
- pj->primitives.gradients.rho[1], pj->primitives.gradients.rho[2]);
- printf("gradWR[1]: %g %g %g\n", pj->primitives.gradients.v[0][0],
- pj->primitives.gradients.v[0][1], pj->primitives.gradients.v[0][2]);
- printf("gradWR[2]: %g %g %g\n", pj->primitives.gradients.v[1][0],
- pj->primitives.gradients.v[1][1], pj->primitives.gradients.v[1][2]);
- printf("gradWR[3]: %g %g %g\n", pj->primitives.gradients.v[2][0],
- pj->primitives.gradients.v[2][1], pj->primitives.gradients.v[2][2]);
- printf("gradWR[4]: %g %g %g\n", pj->primitives.gradients.P[0],
- pj->primitives.gradients.P[1], pj->primitives.gradients.P[2]);
- printf("WR': %g %g %g %g %g\n", Wj[0], Wj[1], Wj[2], Wj[3], Wj[4]);
- error("Negative density or pressure!\n");
- }
-
float totflux[5];
riemann_solve_for_flux(Wi, Wj, n_unit, vij, totflux);
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index 236106a1fb04cc2e5b84f997a2389d583ce17cff..3d58be2f47c4e1904aaac5f69d1862f1d453e488 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -127,7 +127,7 @@ float convert_Etot(struct engine* e, struct part* p) {
void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) {
- *num_fields = 14;
+ *num_fields = 11;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -143,22 +143,16 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
parts, primitives.P, convert_u);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
- list[6] = io_make_output_field("Acceleration", FLOAT, 3,
- UNIT_CONV_ACCELERATION, parts, a_hydro);
- list[7] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
+ list[6] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
primitives.rho);
- list[8] = io_make_output_field("Volume", FLOAT, 1, UNIT_CONV_VOLUME, parts,
- geometry.volume);
- list[9] = io_make_output_field("GradDensity", FLOAT, 3, UNIT_CONV_DENSITY,
- parts, primitives.gradients.rho);
- list[10] = io_make_output_field_convert_part(
+ list[7] = io_make_output_field_convert_part(
"Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, primitives.P, convert_A);
- list[11] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
- parts, primitives.P);
- list[12] =
+ list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
+ parts, primitives.P);
+ list[9] =
io_make_output_field_convert_part("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
parts, conserved.energy, convert_Etot);
- list[13] = io_make_output_field("GravAcceleration", FLOAT, 3,
+ list[10] = io_make_output_field("GravAcceleration", FLOAT, 3,
UNIT_CONV_ACCELERATION, parts, gravity.old_a);
}
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index d552a3f7e86031311098293845f1aa11270c417f..6c96004847ae23b46ec3f5182f742e0e84f1118d 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/Gizmo/hydro_part.h
@@ -148,6 +148,9 @@ struct part {
/* Total surface area of the particle. */
float Atot;
+ /* Centroid of the "cell". */
+ float centroid[3];
+
} geometry;
/* Variables used for timestep calculation (currently not used). */
@@ -201,6 +204,8 @@ struct part {
/* Previous value of the gravitational acceleration. */
float old_a[3];
+ float grad_a[3][3];
+
/* Previous value of the mass flux vector. */
float old_mflux[3];
diff --git a/src/hydro/Gizmo/hydro_slope_limiters_face.h b/src/hydro/Gizmo/hydro_slope_limiters_face.h
index 7ae5dd2eb073d9aae8ab6f2efffdf8df15b4bb4a..ba96063d661a93a4efc4069ff7e7269a4ac58c3b 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters_face.h
+++ b/src/hydro/Gizmo/hydro_slope_limiters_face.h
@@ -53,14 +53,22 @@ hydro_slope_limit_face_quantity(float phi_i, float phi_j, float phi_mid0,
if ((phimax + delta1) * phimax > 0.0f) {
phiplus = phimax + delta1;
} else {
- phiplus = phimax / (1.0f + delta1 / fabs(phimax));
+ if (phimax != 0.) {
+ phiplus = phimax / (1.0f + delta1 / fabs(phimax));
+ } else {
+ phiplus = 0.;
+ }
}
/* if sign(phimin-delta1) == sign(phimin) */
if ((phimin - delta1) * phimin > 0.0f) {
phiminus = phimin - delta1;
} else {
- phiminus = phimin / (1.0f + delta1 / fabs(phimin));
+ if (phimin != 0.) {
+ phiminus = phimin / (1.0f + delta1 / fabs(phimin));
+ } else {
+ phiminus = 0.;
+ }
}
if (phi_i < phi_j) {
diff --git a/src/hydro/Gizmo/hydro_unphysical.h b/src/hydro/Gizmo/hydro_unphysical.h
new file mode 100644
index 0000000000000000000000000000000000000000..517e3e0918ad340580e270477c0a166590546850
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_unphysical.h
@@ -0,0 +1,55 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_HYDRO_UNPHYSICAL_H
+#define SWIFT_HYDRO_UNPHYSICAL_H
+
+#if defined(GIZMO_UNPHYSICAL_ERROR) || defined(GIZMO_UNPHYSICAL_RESCUE)
+
+#if defined(GIZMO_UNPHYSICAL_ERROR)
+
+/*! @brief Crash whenever an unphysical value is detected. */
+#define gizmo_unphysical_message(name, quantity) \
+ error("Unphysical " name " detected (%g)!", quantity);
+
+#elif defined(GIZMO_UNPHYSICAL_WARNING)
+
+/*! @brief Show a warning whenever an unphysical value is detected. */
+#define gizmo_unphysical_message(name, quantity) \
+ message("Unphysical " name " detected (%g), reset to 0!", quantity);
+
+#else
+
+/*! @brief Don't tell anyone an unphysical value was detected. */
+#define gizmo_unphysical_message(name, quantity)
+
+#endif
+
+#define gizmo_check_physical_quantity(name, quantity) \
+ if (quantity < 0.) { \
+ gizmo_unphysical_message(name, quantity); \
+ quantity = 0.; \
+ }
+
+#else // defined(GIZMO_UNPHYSICAL_ERROR) || defined(GIZMO_UNPHYSICAL_RESCUE)
+
+#define gizmo_check_physical_quantity(name, quantity)
+
+#endif // defined(GIZMO_UNPHYSICAL_ERROR) || defined(GIZMO_UNPHYSICAL_RESCUE)
+
+#endif // SWIFT_HYDRO_UNPHYSICAL_H
diff --git a/src/hydro/Gizmo/hydro_velocities.h b/src/hydro/Gizmo/hydro_velocities.h
new file mode 100644
index 0000000000000000000000000000000000000000..a0c896d600aaf971bc12732c294d8e656698ed4c
--- /dev/null
+++ b/src/hydro/Gizmo/hydro_velocities.h
@@ -0,0 +1,159 @@
+/*******************************************************************************
+ * 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.;
+ p->v[1] = 0.;
+ p->v[2] = 0.;
+#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 particel data to act upon.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_velocities_prepare_force(struct part* restrict p,
+ const struct xpart* restrict xp) {
+
+#ifndef GIZMO_FIX_PARTICLES
+ p->force.v_full[0] = xp->v_full[0];
+ p->force.v_full[1] = xp->v_full[1];
+ p->force.v_full[2] = xp->v_full[2];
+#endif
+}
+
+/**
+ * @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.;
+ p->v[1] = 0.;
+ p->v[2] = 0.;
+
+#else // GIZMO_FIX_PARTICLES
+
+ if (p->conserved.mass > 0. && p->primitives.rho > 0.) {
+ /* Normal case: set particle velocity to fluid velocity. */
+ p->v[0] = p->conserved.momentum[0] / p->conserved.mass;
+ p->v[1] = p->conserved.momentum[1] / p->conserved.mass;
+ p->v[2] = p->conserved.momentum[2] / p->conserved.mass;
+ } else {
+ /* Vacuum particles have no fluid velocity. */
+ p->v[0] = 0.;
+ p->v[1] = 0.;
+ p->v[2] = 0.;
+ }
+
+#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.25;
+ const float etaR = eta * R;
+ const float xi = 1.;
+ 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.9 * etaR) {
+ float fac = xi * soundspeed / d;
+ if (d < 1.1 * etaR) {
+ fac *= 5. * (d - 0.9 * etaR) / etaR;
+ }
+ p->v[0] -= ds[0] * fac;
+ p->v[1] -= ds[1] * fac;
+ p->v[2] -= ds[2] * fac;
+ }
+
+#endif // GIZMO_STEER_MOTION
+
+#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.h b/src/riemann.h
index 685d40708e598249151f6cbe13be016edea79553..ab6d162514326778e8d6478e07c9bae2947a7c2a 100644
--- a/src/riemann.h
+++ b/src/riemann.h
@@ -25,10 +25,8 @@
#if defined(RIEMANN_SOLVER_EXACT)
#define RIEMANN_SOLVER_IMPLEMENTATION "Exact Riemann solver (Toro 2009)"
-#if defined(EOS_IDEAL_GAS)
+#if defined(EOS_IDEAL_GAS) || defined(EOS_ISOTHERMAL_GAS)
#include "riemann/riemann_exact.h"
-#elif defined(EOS_ISOTHERMAL_GAS)
-#include "riemann/riemann_exact_isothermal.h"
#else
#error "The Exact Riemann solver is incompatible with this equation of state!"
#endif