diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index e84a434230d83adfd1c9369a947feea48eaef44b..a15a09b4f06e6e1b6e082a48dfaa3371049775de 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -43,11 +43,11 @@ Statistics:
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
- delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
- max_volume_change: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step
+ h_tolerance: 1e-4 # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths.
h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
+ max_volume_change: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step.
+ max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
# Parameters for the self-gravity scheme
Gravity:
diff --git a/src/cell.c b/src/cell.c
index a9cc9cc7e1ca05404b24fc196f725ece26e5a3fc..be071f570c20e6a13363e75d2356e574e51a58a6 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1718,7 +1718,7 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
else
error("Drift task missing !");
if (t->type == task_type_pair) scheduler_activate(s, cj->drift_part);
-
+
if (cell_is_active(cj, e)) {
for (l = cj->send_rho; l != NULL && l->t->cj->nodeID != ci->nodeID;
@@ -1870,11 +1870,10 @@ void cell_drift_part(struct cell *c, const struct engine *e) {
float dx_max = 0.f, dx2_max = 0.f;
float dx_max_sort = 0.0f, dx2_max_sort = 0.f;
float cell_h_max = 0.f;
-
+
#ifdef SWIFT_DEBUG_CHECKS
/* Check that we only drift local cells. */
- if (c->nodeID != engine_rank)
- error("Drifting a foreign cell is nope.");
+ if (c->nodeID != engine_rank) error("Drifting a foreign cell is nope.");
/* Check that we are actually going to move forward. */
if (ti_current < ti_old_part) error("Attempt to drift to the past");
diff --git a/src/dimension.h b/src/dimension.h
index 7c58b5b846490e8f5227a232eaaeb3df5995f795..4ee9377b4f9608b2549f3dfee2a683b61a71c76a 100644
--- a/src/dimension.h
+++ b/src/dimension.h
@@ -118,6 +118,34 @@ __attribute__((always_inline)) INLINE static float pow_dimension_plus_one(
#endif
}
+/**
+ * @brief Returns the argument to the power given by the dimension minus one
+ *
+ * Computes \f$x^{d-1}\f$.
+ */
+__attribute__((always_inline)) INLINE static float pow_dimension_minus_one(
+ float x) {
+
+#if defined(HYDRO_DIMENSION_3D)
+
+ return x * x;
+
+#elif defined(HYDRO_DIMENSION_2D)
+
+ return x;
+
+#elif defined(HYDRO_DIMENSION_1D)
+
+ return 1.f;
+
+#else
+
+ error("The dimension is not defined !");
+ return 0.f;
+
+#endif
+}
+
/**
* @brief Inverts the given dimension by dimension matrix (in place)
*
diff --git a/src/engine.c b/src/engine.c
index 68a34cc25150bb2d5a71fa32750abe222fbea676..8e11fe14beb0c8f52db2ceae6c0f6ec3f8ab1720 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -4266,7 +4266,7 @@ void engine_init(struct engine *e, struct space *s,
"Version: %s \n# "
"Number of threads: %d\n# Number of MPI ranks: %d\n# Hydrodynamic "
"scheme: %s\n# Hydrodynamic kernel: %s\n# No. of neighbours: %.2f "
- "+/- %.2f\n# Eta: %f\n",
+ "+/- %.4f\n# Eta: %f\n",
hostname(), git_branch(), git_revision(), compiler_name(),
compiler_version(), e->nr_threads, e->nr_nodes, SPH_IMPLEMENTATION,
kernel_name, e->hydro_properties->target_neighbours,
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 91626749a89ede387547b6351dce59fa3569307a..66a475f32ec06eb40ff2bc890bc156f76e3b7b9f 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -206,12 +206,13 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->rho += p->mass * kernel_root;
p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
p->density.wcount += kernel_root;
+ p->density.wcount_dh -= hydro_dimension * kernel_root;
/* Finish the calculation by inserting the missing h-factors */
p->rho *= h_inv_dim;
p->density.rho_dh *= h_inv_dim_plus_one;
- p->density.wcount *= kernel_norm;
- p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
+ p->density.wcount *= h_inv_dim;
+ p->density.wcount_dh *= h_inv_dim_plus_one;
const float rho_inv = 1.f / p->rho;
@@ -224,6 +225,31 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.div_v *= h_inv_dim_plus_one * rho_inv;
}
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
+ struct part *restrict p, struct xpart *restrict xp) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+
+ /* Re-set problematic values */
+ p->rho = p->mass * kernel_root * h_inv_dim;
+ p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+ p->density.rho_dh = 0.f;
+ p->density.wcount_dh = 0.f;
+ p->density.div_v = 0.f;
+ p->density.rot_v[0] = 0.f;
+ p->density.rot_v[1] = 0.f;
+ p->density.rot_v[2] = 0.f;
+}
+
/**
* @brief Prepare a particle for the force calculation.
*
@@ -239,6 +265,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float fac_mu = 1.f; /* Will change with cosmological integration */
+ /* Inverse of the physical density */
+ const float rho_inv = 1.f / p->rho;
+
/* Compute the norm of the curl */
const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
p->density.rot_v[1] * p->density.rot_v[1] +
@@ -254,7 +283,6 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
/* Divide the pressure by the density squared to get the SPH term */
- const float rho_inv = 1.f / p->rho;
const float P_over_rho2 = pressure * rho_inv * rho_inv;
/* Compute the Balsara switch */
@@ -262,11 +290,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);
/* Compute the "grad h" term */
- const float grad_h_term =
+ const float omega_inv =
1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
/* Update variables. */
- p->force.f = grad_h_term;
+ p->force.f = omega_inv;
p->force.P_over_rho2 = P_over_rho2;
p->force.soundspeed = soundspeed;
p->force.balsara = balsara;
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index ae4ee27f5ed83350aaf1a0bfd0069d7b85701854..25eb1e4b76d088d7dcd3d4ba57fe5901188be000 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -64,7 +64,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the number of neighbours */
pi->density.wcount += wi;
- pi->density.wcount_dh -= ui * wi_dx;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
/* Compute the kernel function for pj */
const float hj_inv = 1.f / hj;
@@ -77,7 +77,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the number of neighbours */
pj->density.wcount += wj;
- pj->density.wcount_dh -= uj * wj_dx;
+ pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
const float faci = mj * wi_dx * r_inv;
const float facj = mi * wj_dx * r_inv;
@@ -112,9 +112,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
struct part **pj) {
-#ifdef WITH_VECTORIZATION
+#ifdef WITH_OLD_VECTORIZATION
- vector r, ri, r2, xi, xj, hi, hj, hi_inv, hj_inv, wi, wj, wi_dx, wj_dx;
+ vector r, ri, r2, ui, uj, hi, hj, hi_inv, hj_inv, wi, wj, wi_dx, wj_dx;
vector rhoi, rhoj, rhoi_dh, rhoj_dh, wcounti, wcountj, wcounti_dh, wcountj_dh;
vector mi, mj;
vector dx[3], dv[3];
@@ -161,15 +161,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
hi.v = vec_load(Hi);
hi_inv = vec_reciprocal(hi);
- xi.v = r.v * hi_inv.v;
+ ui.v = r.v * hi_inv.v;
hj.v = vec_load(Hj);
hj_inv = vec_reciprocal(hj);
- xj.v = r.v * hj_inv.v;
+ uj.v = r.v * hj_inv.v;
/* Compute the kernel function. */
- kernel_deval_vec(&xi, &wi, &wi_dx);
- kernel_deval_vec(&xj, &wj, &wj_dx);
+ kernel_deval_vec(&ui, &wi, &wi_dx);
+ kernel_deval_vec(&uj, &wj, &wj_dx);
/* Compute dv. */
dv[0].v = vi[0].v - vj[0].v;
@@ -188,17 +188,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
/* Compute density of pi. */
rhoi.v = mj.v * wi.v;
- rhoi_dh.v = mj.v * (vec_set1(hydro_dimension) * wi.v + xi.v * wi_dx.v);
+ rhoi_dh.v = mj.v * (vec_set1(hydro_dimension) * wi.v + ui.v * wi_dx.v);
wcounti.v = wi.v;
- wcounti_dh.v = xi.v * wi_dx.v;
+ wcounti_dh.v = (vec_set1(hydro_dimension) * wi.v + ui.v * wi_dx.v);
div_vi.v = mj.v * dvdr.v * wi_dx.v;
for (k = 0; k < 3; k++) curl_vi[k].v = mj.v * curlvr[k].v * wi_dx.v;
/* Compute density of pj. */
rhoj.v = mi.v * wj.v;
- rhoj_dh.v = mi.v * (vec_set1(hydro_dimension) * wj.v + xj.v * wj_dx.v);
+ rhoj_dh.v = mi.v * (vec_set1(hydro_dimension) * wj.v + uj.v * wj_dx.v);
wcountj.v = wj.v;
- wcountj_dh.v = xj.v * wj_dx.v;
+ wcountj_dh.v = (vec_set1(hydro_dimension) * wj.v + uj.v * wj_dx.v);
div_vj.v = mi.v * dvdr.v * wj_dx.v;
for (k = 0; k < 3; k++) curl_vj[k].v = mi.v * curlvr[k].v * wj_dx.v;
@@ -241,7 +241,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* Get r and r inverse. */
const float r = sqrtf(r2);
- const float ri = 1.0f / r;
+ const float r_inv = 1.0f / r;
/* Compute the kernel function */
const float hi_inv = 1.0f / hi;
@@ -254,9 +254,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* Compute contribution to the number of neighbours */
pi->density.wcount += wi;
- pi->density.wcount_dh -= ui * wi_dx;
+ pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
- const float fac = mj * wi_dx * ri;
+ const float fac = mj * wi_dx * r_inv;
/* Compute dv dot r */
dv[0] = pi->v[0] - pj->v[0];
@@ -282,9 +282,9 @@ __attribute__((always_inline)) INLINE static void
runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
struct part **pi, struct part **pj) {
-#ifdef WITH_VECTORIZATION
+#ifdef WITH_OLD_VECTORIZATION
- vector r, ri, r2, xi, hi, hi_inv, wi, wi_dx;
+ vector r, ri, r2, ui, hi, hi_inv, wi, wi_dx;
vector rhoi, rhoi_dh, wcounti, wcounti_dh, div_vi;
vector mj;
vector dx[3], dv[3];
@@ -328,9 +328,9 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
hi.v = vec_load(Hi);
hi_inv = vec_reciprocal(hi);
- xi.v = r.v * hi_inv.v;
+ ui.v = r.v * hi_inv.v;
- kernel_deval_vec(&xi, &wi, &wi_dx);
+ kernel_deval_vec(&ui, &wi, &wi_dx);
/* Compute dv. */
dv[0].v = vi[0].v - vj[0].v;
@@ -349,9 +349,9 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
/* Compute density of pi. */
rhoi.v = mj.v * wi.v;
- rhoi_dh.v = mj.v * (vec_set1(hydro_dimension) * wi.v + xi.v * wi_dx.v);
+ rhoi_dh.v = mj.v * (vec_set1(hydro_dimension) * wi.v + ui.v * wi_dx.v);
wcounti.v = wi.v;
- wcounti_dh.v = xi.v * wi_dx.v;
+ wcounti_dh.v = (vec_set1(hydro_dimension) * wi.v + ui.v * wi_dx.v);
div_vi.v = mj.v * dvdr.v * wi_dx.v;
for (k = 0; k < 3; k++) curl_vi[k].v = mj.v * curlvr[k].v * wi_dx.v;
@@ -390,7 +390,7 @@ runner_iact_nonsym_1_vec_density(vector *r2, vector *dx, vector *dy, vector *dz,
vector *curlvySum, vector *curlvzSum,
vector mask, int knlMask) {
- vector r, ri, xi, wi, wi_dx;
+ vector r, ri, ui, wi, wi_dx;
vector mj;
vector dvx, dvy, dvz;
vector vjx, vjy, vjz;
@@ -407,10 +407,10 @@ runner_iact_nonsym_1_vec_density(vector *r2, vector *dx, vector *dy, vector *dz,
ri = vec_reciprocal_sqrt(*r2);
r.v = vec_mul(r2->v, ri.v);
- xi.v = vec_mul(r.v, hi_inv.v);
+ ui.v = vec_mul(r.v, hi_inv.v);
/* Calculate the kernel for two particles. */
- kernel_deval_1_vec(&xi, &wi, &wi_dx);
+ kernel_deval_1_vec(&ui, &wi, &wi_dx);
/* Compute dv. */
dvx.v = vec_sub(vix.v, vjx.v);
@@ -440,12 +440,13 @@ runner_iact_nonsym_1_vec_density(vector *r2, vector *dx, vector *dy, vector *dz,
rho_dhSum->v =
_mm512_mask_sub_ps(rho_dhSum->v, knlMask, rho_dhSum->v,
vec_mul(mj.v, vec_fma(vec_set1(hydro_dimension), wi.v,
- vec_mul(xi.v, wi_dx.v))));
+ vec_mul(ui.v, wi_dx.v))));
wcountSum->v = _mm512_mask_add_ps(wcountSum->v, knlMask, wi.v, wcountSum->v);
- wcount_dhSum->v = _mm512_mask_sub_ps(wcount_dhSum->v, knlMask,
- wcount_dhSum->v, vec_mul(xi.v, wi_dx.v));
+ wcount_dhSum->v = _mm512_mask_sub_ps(
+ rho_dhSum->v, knlMask, rho_dhSum->v,
+ vec_fma(vec_set1(hydro_dimension), wi.v, vec_mul(ui.v, wi_dx.v)));
div_vSum->v = _mm512_mask_sub_ps(div_vSum->v, knlMask, div_vSum->v,
vec_mul(mj.v, vec_mul(dvdr.v, wi_dx.v)));
@@ -464,10 +465,11 @@ runner_iact_nonsym_1_vec_density(vector *r2, vector *dx, vector *dy, vector *dz,
#else
rhoSum->v += vec_and(vec_mul(mj.v, wi.v), mask.v);
rho_dhSum->v -= vec_and(vec_mul(mj.v, vec_fma(vec_set1(hydro_dimension), wi.v,
- vec_mul(xi.v, wi_dx.v))),
+ vec_mul(ui.v, wi_dx.v))),
mask.v);
wcountSum->v += vec_and(wi.v, mask.v);
- wcount_dhSum->v -= vec_and(vec_mul(xi.v, wi_dx.v), mask.v);
+ wcount_dhSum->v -= vec_and(
+ vec_fma(vec_set1(hydro_dimension), wi.v, vec_mul(ui.v, wi_dx.v)), mask.v);
div_vSum->v -= vec_and(vec_mul(mj.v, vec_mul(dvdr.v, wi_dx.v)), mask.v);
curlvxSum->v += vec_and(vec_mul(mj.v, vec_mul(curlvrx.v, wi_dx.v)), mask.v);
curlvySum->v += vec_and(vec_mul(mj.v, vec_mul(curlvry.v, wi_dx.v)), mask.v);
@@ -487,13 +489,13 @@ runner_iact_nonsym_2_vec_density(
vector *div_vSum, vector *curlvxSum, vector *curlvySum, vector *curlvzSum,
vector mask, vector mask2, int knlMask, int knlMask2) {
- vector r, ri, r2, xi, wi, wi_dx;
+ vector r, ri, r2, ui, wi, wi_dx;
vector mj;
vector dx, dy, dz, dvx, dvy, dvz;
vector vjx, vjy, vjz;
vector dvdr;
vector curlvrx, curlvry, curlvrz;
- vector r_2, ri2, r2_2, xi2, wi2, wi_dx2;
+ vector r_2, ri2, r2_2, ui2, wi2, wi_dx2;
vector mj2;
vector dx2, dy2, dz2, dvx2, dvy2, dvz2;
vector vjx2, vjy2, vjz2;
@@ -524,11 +526,11 @@ runner_iact_nonsym_2_vec_density(
r.v = vec_mul(r2.v, ri.v);
r_2.v = vec_mul(r2_2.v, ri2.v);
- xi.v = vec_mul(r.v, hi_inv.v);
- xi2.v = vec_mul(r_2.v, hi_inv.v);
+ ui.v = vec_mul(r.v, hi_inv.v);
+ ui2.v = vec_mul(r_2.v, hi_inv.v);
/* Calculate the kernel for two particles. */
- kernel_deval_2_vec(&xi, &wi, &wi_dx, &xi2, &wi2, &wi_dx2);
+ kernel_deval_2_vec(&ui, &wi, &wi_dx, &ui2, &wi2, &wi_dx2);
/* Compute dv. */
dvx.v = vec_sub(vix.v, vjx.v);
@@ -575,20 +577,23 @@ runner_iact_nonsym_2_vec_density(
rho_dhSum->v =
_mm512_mask_sub_ps(rho_dhSum->v, knlMask, rho_dhSum->v,
vec_mul(mj.v, vec_fma(vec_set1(hydro_dimension), wi.v,
- vec_mul(xi.v, wi_dx.v))));
- rho_dhSum->v = _mm512_mask_sub_ps(
- rho_dhSum->v, knlMask2, rho_dhSum->v,
- vec_mul(mj2.v, vec_fma(vec_set1(hydro_dimension), wi2.v,
- vec_mul(xi2.v, wi_dx2.v))));
+ vec_mul(ui.v, wi_dx.v))));
+ rho_dhSum->v =
+ _mm512_mask_sub_ps(rho_dhSum->v, knlMask2, rho_dhSum->v,
+ vec_mul(mj.v, vec_fma(vec_set1(hydro_dimension), wi2.v,
+ vec_mul(ui2.v, wi_dx2.v))));
wcountSum->v = _mm512_mask_add_ps(wcountSum->v, knlMask, wi.v, wcountSum->v);
wcountSum->v =
_mm512_mask_add_ps(wcountSum->v, knlMask2, wi2.v, wcountSum->v);
- wcount_dhSum->v = _mm512_mask_sub_ps(wcount_dhSum->v, knlMask,
- wcount_dhSum->v, vec_mul(xi.v, wi_dx.v));
wcount_dhSum->v = _mm512_mask_sub_ps(
- wcount_dhSum->v, knlMask2, wcount_dhSum->v, vec_mul(xi2.v, wi_dx2.v));
+ rho_dhSum->v, knlMask, rho_dhSum->v,
+ vec_fma(vec_set1(hydro_dimension), wi.v, vec_mul(ui.v, wi_dx.v)));
+
+ wcount_dhSum->v = _mm512_mask_sub_ps(
+ rho_dhSum->v, knlMask2, rho_dhSum->v,
+ vec_fma(vec_set1(hydro_dimension), wi2.v, vec_mul(ui2.v, wi_dx2.v)));
div_vSum->v = _mm512_mask_sub_ps(div_vSum->v, knlMask, div_vSum->v,
vec_mul(mj.v, vec_mul(dvdr.v, wi_dx.v)));
@@ -619,16 +624,19 @@ runner_iact_nonsym_2_vec_density(
rhoSum->v += vec_and(vec_mul(mj.v, wi.v), mask.v);
rhoSum->v += vec_and(vec_mul(mj2.v, wi2.v), mask2.v);
rho_dhSum->v -= vec_and(vec_mul(mj.v, vec_fma(vec_set1(hydro_dimension), wi.v,
- vec_mul(xi.v, wi_dx.v))),
+ vec_mul(ui.v, wi_dx.v))),
mask.v);
rho_dhSum->v -=
vec_and(vec_mul(mj2.v, vec_fma(vec_set1(hydro_dimension), wi2.v,
- vec_mul(xi2.v, wi_dx2.v))),
+ vec_mul(ui2.v, wi_dx2.v))),
mask2.v);
wcountSum->v += vec_and(wi.v, mask.v);
wcountSum->v += vec_and(wi2.v, mask2.v);
- wcount_dhSum->v -= vec_and(vec_mul(xi.v, wi_dx.v), mask.v);
- wcount_dhSum->v -= vec_and(vec_mul(xi2.v, wi_dx2.v), mask2.v);
+ wcount_dhSum->v -= vec_and(
+ vec_fma(vec_set1(hydro_dimension), wi.v, vec_mul(ui.v, wi_dx.v)), mask.v);
+ wcount_dhSum->v -= vec_and(
+ vec_fma(vec_set1(hydro_dimension), wi2.v, vec_mul(ui2.v, wi_dx2.v)),
+ mask2.v);
div_vSum->v -= vec_and(vec_mul(mj.v, vec_mul(dvdr.v, wi_dx.v)), mask.v);
div_vSum->v -= vec_and(vec_mul(mj2.v, vec_mul(dvdr2.v, wi_dx2.v)), mask2.v);
curlvxSum->v += vec_and(vec_mul(mj.v, vec_mul(curlvrx.v, wi_dx.v)), mask.v);
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index c9b0e5ff7f26aa433a9e3492a4282f94cf2f6624..276d6d837bcc164bbc006906b2e632041e033cae 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -35,14 +35,24 @@
#define hydro_props_default_max_iterations 30
#define hydro_props_default_volume_change 1.4f
#define hydro_props_default_h_max FLT_MAX
+#define hydro_props_default_h_tolerance 1e-4
void hydro_props_init(struct hydro_props *p,
const struct swift_params *params) {
/* Kernel properties */
p->eta_neighbours = parser_get_param_float(params, "SPH:resolution_eta");
+
+ /* Tolerance for the smoothing length Newton-Raphson scheme */
+ p->h_tolerance = parser_get_opt_param_float(params, "SPH:h_tolerance",
+ hydro_props_default_h_tolerance);
+
+ /* Get derived properties */
p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
- p->delta_neighbours = parser_get_param_float(params, "SPH:delta_neighbours");
+ const float delta_eta = p->eta_neighbours * (1.f + p->h_tolerance);
+ p->delta_neighbours =
+ (pow_dimension(delta_eta) - pow_dimension(p->eta_neighbours)) *
+ kernel_norm;
#ifdef SHADOWFAX_SPH
/* change the meaning of target_neighbours and delta_neighbours */
@@ -81,9 +91,11 @@ void hydro_props_print(const struct hydro_props *p) {
message("Hydrodynamic scheme: %s in %dD.", SPH_IMPLEMENTATION,
(int)hydro_dimension);
- message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours (eta=%f).",
- kernel_name, p->target_neighbours, p->delta_neighbours,
- p->eta_neighbours);
+ message("Hydrodynamic kernel: %s with eta=%f (%.2f neighbours).", kernel_name,
+ p->eta_neighbours, p->target_neighbours);
+
+ message("Hydrodynamic tolerance in h: %.5f (+/- %.4f neighbours).",
+ p->h_tolerance, p->delta_neighbours);
message("Hydrodynamic integration: CFL parameter: %.4f.", p->CFL_condition);
@@ -110,6 +122,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
io_write_attribute_f(h_grpsph, "Kernel target N_ngb", p->target_neighbours);
io_write_attribute_f(h_grpsph, "Kernel delta N_ngb", p->delta_neighbours);
io_write_attribute_f(h_grpsph, "Kernel eta", p->eta_neighbours);
+ io_write_attribute_f(h_grpsph, "Smoothing length tolerance", p->h_tolerance);
io_write_attribute_f(h_grpsph, "Maximal smoothing length", p->h_max);
io_write_attribute_f(h_grpsph, "CFL parameter", p->CFL_condition);
io_write_attribute_f(h_grpsph, "Volume log(max(delta h))",
diff --git a/src/hydro_properties.h b/src/hydro_properties.h
index 716c4c060c21eb95d05f9d50e13d4681a958a6fd..a887ccb6df13b649cd1ef1009059c6f08908669c 100644
--- a/src/hydro_properties.h
+++ b/src/hydro_properties.h
@@ -16,10 +16,14 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-
#ifndef SWIFT_HYDRO_PROPERTIES
#define SWIFT_HYDRO_PROPERTIES
+/**
+ * @file hydro_properties.h
+ * @brief Contains all the constants and parameters of the hydro scheme
+ */
+
/* Config parameters. */
#include "../config.h"
@@ -35,19 +39,28 @@
*/
struct hydro_props {
- /* Kernel properties */
+ /*! Resolution parameter */
float eta_neighbours;
+
+ /*! Target weightd number of neighbours (for info only)*/
float target_neighbours;
+
+ /*! Smoothing length tolerance */
+ float h_tolerance;
+
+ /*! Tolerance on neighbour number (for info only)*/
float delta_neighbours;
- /* Maximal smoothing length */
+ /*! Maximal smoothing length */
float h_max;
- /* Number of iterations to converge h */
+ /*! Maximal number of iterations to converge h */
int max_smoothing_iterations;
- /* Time integration properties */
+ /*! Time integration properties */
float CFL_condition;
+
+ /*! Maximal change of h over one time-step */
float log_max_h_change;
};
diff --git a/src/runner.c b/src/runner.c
index dc0879b5ef6195a04ce2ee7b7672a5789da190c1..4588383bf3e3218a36d567402877ad1930b03649 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -622,11 +622,9 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const struct space *s = e->s;
const float hydro_h_max = e->hydro_properties->h_max;
- const float target_wcount = e->hydro_properties->target_neighbours;
- const float max_wcount =
- target_wcount + e->hydro_properties->delta_neighbours;
- const float min_wcount =
- target_wcount - e->hydro_properties->delta_neighbours;
+ const float eps = e->hydro_properties->h_tolerance;
+ const float hydro_eta_dim =
+ pow_dimension(e->hydro_properties->eta_neighbours);
const int max_smoothing_iter = e->hydro_properties->max_smoothing_iterations;
int redo = 0, count = 0;
@@ -670,28 +668,47 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
if (!part_is_active(p, e)) error("Ghost applied to inactive particle");
#endif
- /* Finish the density calculation */
- hydro_end_density(p);
+ /* Get some useful values */
+ const float h_old = p->h;
+ const float h_old_dim = pow_dimension(h_old);
+ const float h_old_dim_minus_one = pow_dimension_minus_one(h_old);
+ float h_new;
- /* Did we get the right number of neighbours? */
- if (p->density.wcount > max_wcount || p->density.wcount < min_wcount) {
+ if (p->density.wcount == 0.f) { /* No neighbours case */
- float h_corr = 0.f;
+ /* Double h and try again */
+ h_new = 2.f * h_old;
+ } else {
- /* If no derivative, double the smoothing length. */
- if (p->density.wcount_dh == 0.0f) h_corr = p->h;
+ /* Finish the density calculation */
+ hydro_end_density(p);
- /* Otherwise, compute the smoothing length update (Newton step). */
- else {
- h_corr = (target_wcount - p->density.wcount) / p->density.wcount_dh;
+ /* Compute one step of the Newton-Raphson scheme */
+ const float n_sum = p->density.wcount * h_old_dim;
+ const float n_target = hydro_eta_dim;
+ const float f = n_sum - n_target;
+ const float f_prime =
+ p->density.wcount_dh * h_old_dim +
+ hydro_dimension * p->density.wcount * h_old_dim_minus_one;
- /* Truncate to the range [ -p->h/2 , p->h ]. */
- h_corr = (h_corr < p->h) ? h_corr : p->h;
- h_corr = (h_corr > -0.5f * p->h) ? h_corr : -0.5f * p->h;
- }
+ h_new = h_old - f / f_prime;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if ((f > 0.f && h_new > h_old) || (f < 0.f && h_new < h_old))
+ error(
+ "Smoothing length correction not going in the right direction");
+#endif
+
+ /* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */
+ h_new = min(h_new, 2.f * h_old);
+ h_new = max(h_new, 0.5f * h_old);
+ }
+
+ /* Check whether the particle has an inappropriate smoothing length */
+ if (fabsf(h_new - h_old) > eps * h_old) {
/* Ok, correct then */
- p->h += h_corr;
+ p->h = h_new;
/* If below the absolute maximum, try again */
if (p->h < hydro_h_max) {
@@ -709,6 +726,10 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* Ok, this particle is a lost cause... */
p->h = hydro_h_max;
+
+ /* Do some damage control if no neighbours at all were found */
+ if (p->density.wcount == kernel_root * kernel_norm)
+ hydro_part_has_no_neighbours(p, xp);
}
}
diff --git a/tests/test125cells.c b/tests/test125cells.c
index ed1bc6617db685de0091e4baeb9143cdf3dccaf1..f76fb3f348f9b9a59994020305c215fd0a25faca 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -549,8 +549,7 @@ int main(int argc, char *argv[]) {
prog_const.const_newton_G = 1.f;
struct hydro_props hp;
- hp.target_neighbours = pow_dimension(h) * kernel_norm;
- hp.delta_neighbours = 4.;
+ hp.h_tolerance = 1e0;
hp.h_max = FLT_MAX;
hp.max_smoothing_iterations = 1;
hp.CFL_condition = 0.1;
diff --git a/tests/tolerance_27_normal.dat b/tests/tolerance_27_normal.dat
index 3ec58173c9c5a0465382de9094b5d15b81b88a2a..2664972b165444249a05e5a3381995d37780f713 100644
--- a/tests/tolerance_27_normal.dat
+++ b/tests/tolerance_27_normal.dat
@@ -1,4 +1,4 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 4e-5 2e-4 2e-3 1e-5 6e-6 6e-6 6e-6
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 4e-5 3e-4 1e-2 1e-5 6e-6 6e-6 6e-6
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.2e-4 1e-4 1e-4 2e-4 1e-4 1e-4 1e-4
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6
diff --git a/tests/tolerance_27_perturbed.dat b/tests/tolerance_27_perturbed.dat
index ed7df2ca9f55dd6211ed6321ffc40faa37d98092..2bc26dc8ecf590ffa674f3df1535c548e4ecbfa9 100644
--- a/tests/tolerance_27_perturbed.dat
+++ b/tests/tolerance_27_perturbed.dat
@@ -1,4 +1,4 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 1e-4 2e-4 3e-3 1e-5 3e-6 3e-6 7e-6
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-3 1e-5 1e-4 6e-5 2e-3 2e-3 2e-3
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 4e-4 1e-6 1e-6 1e-6 2e-6 2e-6 2e-6
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 1e-4 2e-4 1e-2 1e-5 3e-6 3e-6 7e-6
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-3 1e-5 1e-3 6e-5 2e-3 2e-3 2e-3
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 4e-4 1e-6 1e0 1e-6 2e-6 2e-6 2e-6
diff --git a/tests/tolerance_27_perturbed_h.dat b/tests/tolerance_27_perturbed_h.dat
index 8a83dcb8c460a1bf673e1918e082de446248ec69..4b9ce9562d77ffceee38e217d965c252338e767e 100644
--- a/tests/tolerance_27_perturbed_h.dat
+++ b/tests/tolerance_27_perturbed_h.dat
@@ -1,4 +1,4 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 1e-4 2e-4 4e-3 1e-5 3e-6 3e-6 7e-6
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.4e-2 1e-5 1e-4 2.5e-4 3e-3 3e-3 3e-3
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 4e-6 4e-6 4e-6
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 1e-4 4e-4 1.2e-2 1e-5 3e-6 3e-6 7e-6
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.4e-2 1e-5 1e-3 2.5e-4 3e-3 3e-3 3e-3
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e0 1e-6 4e-6 4e-6 4e-6
diff --git a/tests/tolerance_periodic_BC_normal.dat b/tests/tolerance_periodic_BC_normal.dat
index ce24743d0668a049ac13c53ba5c921e7a46ed099..823e4af488b343f57e3c90e89ee2d4f13d3ca94b 100644
--- a/tests/tolerance_periodic_BC_normal.dat
+++ b/tests/tolerance_periodic_BC_normal.dat
@@ -1,4 +1,4 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 4e-6 4e-5 3e-4 3e-3 2e-4 2e-4 2e-4 2e-4
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 2e-4 1e-4 1e-4 6e-4 2e-3 2e-3 2e-3
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 4e-6 4e-5 1e-3 1e-2 2e-4 2e-4 2e-4 2e-4
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 2e-4 1e-4 2e-4 6e-4 2e-3 2e-3 2e-3
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-4 1e-6 1e-4 5e-4 2e-4 2e-4 2e-4
diff --git a/tests/tolerance_periodic_BC_perturbed.dat b/tests/tolerance_periodic_BC_perturbed.dat
index 78e573ee01811bfdc5c72f59b20f17fe00a54348..df5ee6458ba05eed08006586514467fcdb715990 100644
--- a/tests/tolerance_periodic_BC_perturbed.dat
+++ b/tests/tolerance_periodic_BC_perturbed.dat
@@ -1,4 +1,4 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 3e-6 4e-5 2e-4 3e-3 2e-4 1e-4 1e-4 1e-4
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 6e-3 1e-4 1e-4 1e-2 6e-3 6e-3 6e-3
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-3 1e-6 1e-4 5e-4 3e-3 3e-3 3e-3
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 3e-6 4e-5 1e-3 1e-2 2e-4 1e-4 1e-4 1e-4
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 6e-3 1e-4 3e-3 1e-2 6e-3 6e-3 6e-3
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-3 1e-6 1e0 5e-4 3e-3 3e-3 3e-3