diff --git a/.gitignore b/.gitignore
index e883fac300bb92a5fade0c778b513e709cbd2854..83abbfcc834bb5cac391a960fb25c982a04c40b4 100644
--- a/.gitignore
+++ b/.gitignore
@@ -78,6 +78,7 @@ tests/testRiemannTRRS
tests/testRiemannHLLC
tests/testMatrixInversion
tests/testVoronoi1D
+tests/threadpool_test
theory/latex/swift.pdf
theory/kernel/kernels.pdf
diff --git a/examples/SedovBlast_1D/sedov.yml b/examples/SedovBlast_1D/sedov.yml
index 6f519835d26ff5aa851ffb8999e650815c522cd3..4748d8fd665a65f87934f666cd880277daf92da6 100644
--- a/examples/SedovBlast_1D/sedov.yml
+++ b/examples/SedovBlast_1D/sedov.yml
@@ -11,7 +11,7 @@ TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 5e-2 # The end time of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
+ dt_max: 1e-5 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
diff --git a/src/engine.c b/src/engine.c
index bd838cf3c337e05cd6f1b3eec4d7ba2f11592490..84e5655dd53e42b75104c9b68744ede74ac378c3 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2724,7 +2724,7 @@ void engine_step(struct engine *e) {
mask |= 1 << task_type_sub_pair;
mask |= 1 << task_type_ghost;
mask |= 1 << task_type_extra_ghost; /* Adding unnecessary things to the mask
- does not harm */
+ does not harm */
submask |= 1 << task_subtype_density;
submask |= 1 << task_subtype_gradient;
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 3700d4b3dd426d06859adf6feb9a6f0af85e3af7..258115f48a8c5939599ba3271e19c0d88c8f8b89 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -36,7 +36,9 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const float CFL_condition = hydro_properties->CFL_condition;
- return CFL_condition * p->h / fabsf(p->timestepvars.vmax);
+ float R = get_radius_dimension_sphere(p->cell.volume);
+
+ return CFL_condition * R / fabsf(p->timestepvars.vmax);
}
/**
@@ -79,7 +81,6 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->density.wcount = 0.0f;
p->density.wcount_dh = 0.0f;
p->rho = 0.0f;
- p->geometry.volume = 0.0f;
voronoi_cell_init(&p->cell, p->x);
}
@@ -133,7 +134,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.wcount_dh = 1.0f / (1.1f * hnew - p->h);
}
volume = p->cell.volume;
- p->geometry.volume = volume;
/* compute primitive variables */
/* eqns (3)-(5) */
@@ -239,7 +239,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
float volume;
float m;
float momentum[3];
- volume = p->geometry.volume;
+ volume = p->cell.volume;
p->conserved.mass = m = p->mass;
p->primitives.rho = m / volume;
@@ -343,9 +343,9 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
/* We know the desired velocity; now make sure the particle is accelerated
accordingly during the next time step */
- p->a_hydro[0] = (vcell[0] - p->primitives.v[0]) / p->force.dt;
- p->a_hydro[1] = (vcell[1] - p->primitives.v[1]) / p->force.dt;
- p->a_hydro[2] = (vcell[2] - p->primitives.v[2]) / p->force.dt;
+ p->a_hydro[0] = (vcell[0] - p->force.v_full[0]) / p->force.dt;
+ p->a_hydro[1] = (vcell[1] - p->force.v_full[1]) / p->force.dt;
+ p->a_hydro[2] = (vcell[2] - p->force.v_full[2]) / p->force.dt;
} else {
p->a_hydro[0] = 0.0f;
p->a_hydro[1] = 0.0f;
diff --git a/src/hydro/Shadowswift/hydro_debug.h b/src/hydro/Shadowswift/hydro_debug.h
index db672bab16662961a99558b563b9320d4e75ecc4..fa6c25dd301fd1ed400964a07a9b61a2f2b89131 100644
--- a/src/hydro/Shadowswift/hydro_debug.h
+++ b/src/hydro/Shadowswift/hydro_debug.h
@@ -19,65 +19,65 @@
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
const struct part* p, const struct xpart* xp) {
- printf(
- "x=[%.16e,%.16e,%.16e], "
- "v=[%.3e,%.3e,%.3e], "
- "a=[%.3e,%.3e,%.3e], "
- "h=%.3e, "
- "ti_begin=%d, "
- "ti_end=%d, "
- "primitives={"
- "v=[%.3e,%.3e,%.3e], "
- "rho=%.3e, "
- "P=%.3e, "
- "gradients={"
- "rho=[%.3e,%.3e,%.3e], "
- "v=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]], "
- "P=[%.3e,%.3e,%.3e]}, "
- "limiter={"
- "rho=[%.3e,%.3e], "
- "v=[[%.3e,%.3e],[%.3e,%.3e],[%.3e,%.3e]], "
- "P=[%.3e,%.3e], "
- "maxr=%.3e}}, "
- "conserved={"
- "momentum=[%.3e,%.3e,%.3e], "
- "mass=%.3e, "
- "energy=%.3e}, "
- "geometry={"
- "volume=%.3e, "
- "matrix_E=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]]}, "
- "timestepvars={"
- "vmax=%.3e}, "
- "density={"
- "div_v=%.3e, "
- "wcount_dh=%.3e, "
- "curl_v=[%.3e,%.3e,%.3e], "
- "wcount=%.3e}, "
- "mass=%.3e\n",
- p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
- p->a_hydro[1], p->a_hydro[2], p->h, p->ti_begin, p->ti_end,
- p->primitives.v[0], p->primitives.v[1], p->primitives.v[2],
- p->primitives.rho, p->primitives.P, p->primitives.gradients.rho[0],
- p->primitives.gradients.rho[1], p->primitives.gradients.rho[2],
- p->primitives.gradients.v[0][0], p->primitives.gradients.v[0][1],
- p->primitives.gradients.v[0][2], p->primitives.gradients.v[1][0],
- p->primitives.gradients.v[1][1], p->primitives.gradients.v[1][2],
- p->primitives.gradients.v[2][0], p->primitives.gradients.v[2][1],
- p->primitives.gradients.v[2][2], p->primitives.gradients.P[0],
- p->primitives.gradients.P[1], p->primitives.gradients.P[2],
- p->primitives.limiter.rho[0], p->primitives.limiter.rho[1],
- p->primitives.limiter.v[0][0], p->primitives.limiter.v[0][1],
- p->primitives.limiter.v[1][0], p->primitives.limiter.v[1][1],
- p->primitives.limiter.v[2][0], p->primitives.limiter.v[2][1],
- p->primitives.limiter.P[0], p->primitives.limiter.P[1],
- p->primitives.limiter.maxr, p->conserved.momentum[0],
- p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.mass,
- p->conserved.energy, p->geometry.volume, p->geometry.matrix_E[0][0],
- p->geometry.matrix_E[0][1], p->geometry.matrix_E[0][2],
- p->geometry.matrix_E[1][0], p->geometry.matrix_E[1][1],
- p->geometry.matrix_E[1][2], p->geometry.matrix_E[2][0],
- p->geometry.matrix_E[2][1], p->geometry.matrix_E[2][2],
- p->timestepvars.vmax, p->density.div_v, p->density.wcount_dh,
- p->density.curl_v[0], p->density.curl_v[1], p->density.curl_v[2],
- p->density.wcount, p->mass);
+ /* printf(
+ "x=[%.16e,%.16e,%.16e], "
+ "v=[%.3e,%.3e,%.3e], "
+ "a=[%.3e,%.3e,%.3e], "
+ "h=%.3e, "
+ "ti_begin=%d, "
+ "ti_end=%d, "
+ "primitives={"
+ "v=[%.3e,%.3e,%.3e], "
+ "rho=%.3e, "
+ "P=%.3e, "
+ "gradients={"
+ "rho=[%.3e,%.3e,%.3e], "
+ "v=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]], "
+ "P=[%.3e,%.3e,%.3e]}, "
+ "limiter={"
+ "rho=[%.3e,%.3e], "
+ "v=[[%.3e,%.3e],[%.3e,%.3e],[%.3e,%.3e]], "
+ "P=[%.3e,%.3e], "
+ "maxr=%.3e}}, "
+ "conserved={"
+ "momentum=[%.3e,%.3e,%.3e], "
+ "mass=%.3e, "
+ "energy=%.3e}, "
+ "geometry={"
+ "volume=%.3e, "
+ "matrix_E=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]]}, "
+ "timestepvars={"
+ "vmax=%.3e}, "
+ "density={"
+ "div_v=%.3e, "
+ "wcount_dh=%.3e, "
+ "curl_v=[%.3e,%.3e,%.3e], "
+ "wcount=%.3e}, "
+ "mass=%.3e\n",
+ p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
+ p->a_hydro[1], p->a_hydro[2], p->h, p->ti_begin, p->ti_end,
+ p->primitives.v[0], p->primitives.v[1], p->primitives.v[2],
+ p->primitives.rho, p->primitives.P, p->primitives.gradients.rho[0],
+ p->primitives.gradients.rho[1], p->primitives.gradients.rho[2],
+ p->primitives.gradients.v[0][0], p->primitives.gradients.v[0][1],
+ p->primitives.gradients.v[0][2], p->primitives.gradients.v[1][0],
+ p->primitives.gradients.v[1][1], p->primitives.gradients.v[1][2],
+ p->primitives.gradients.v[2][0], p->primitives.gradients.v[2][1],
+ p->primitives.gradients.v[2][2], p->primitives.gradients.P[0],
+ p->primitives.gradients.P[1], p->primitives.gradients.P[2],
+ p->primitives.limiter.rho[0], p->primitives.limiter.rho[1],
+ p->primitives.limiter.v[0][0], p->primitives.limiter.v[0][1],
+ p->primitives.limiter.v[1][0], p->primitives.limiter.v[1][1],
+ p->primitives.limiter.v[2][0], p->primitives.limiter.v[2][1],
+ p->primitives.limiter.P[0], p->primitives.limiter.P[1],
+ p->primitives.limiter.maxr, p->conserved.momentum[0],
+ p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.mass,
+ p->conserved.energy, p->geometry.volume, p->geometry.matrix_E[0][0],
+ p->geometry.matrix_E[0][1], p->geometry.matrix_E[0][2],
+ p->geometry.matrix_E[1][0], p->geometry.matrix_E[1][1],
+ p->geometry.matrix_E[1][2], p->geometry.matrix_E[2][0],
+ p->geometry.matrix_E[2][1], p->geometry.matrix_E[2][2],
+ p->timestepvars.vmax, p->density.div_v, p->density.wcount_dh,
+ p->density.curl_v[0], p->density.curl_v[1], p->density.curl_v[2],
+ p->density.wcount, p->mass);*/
}
diff --git a/src/hydro/Shadowswift/hydro_gradients.h b/src/hydro/Shadowswift/hydro_gradients.h
index 152993966684d11711a0f9527bc1a0e53881a6dc..a8b62ebd42e13bd6f49d01a9fd317a02158f2d50 100644
--- a/src/hydro/Shadowswift/hydro_gradients.h
+++ b/src/hydro/Shadowswift/hydro_gradients.h
@@ -91,14 +91,13 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
float dWi[5], dWj[5];
float xij_j[3];
- float A[3];
- if (!voronoi_get_face(&pj->cell, pi->id, A, xij_j)) {
- /* Should never happen */
- message("pi->ti_end: %u, pj->ti_end: %u", pi->ti_end, pj->ti_end);
- error("pj is neighbour of pi, but pi not of pj. This should never happen!");
- return;
- }
+ /* xij_j = real_midpoint - pj->x
+ = xij_i + pi->x - pj->x
+ = xij_i + dx */
+ xij_j[0] = xij_i[0] + dx[0];
+ xij_j[1] = xij_i[1] + dx[1];
+ xij_j[2] = xij_i[2] + dx[2];
dWi[0] = pi->primitives.gradients.rho[0] * xij_i[0] +
pi->primitives.gradients.rho[1] * xij_i[1] +
@@ -198,6 +197,21 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
Wj[2] += dWj[2];
Wj[3] += dWj[3];
Wj[4] += dWj[4];
+
+ /* Sanity check: if density or pressure becomes negative after the
+ interpolation, just reset them */
+ if (Wi[0] < 0.0f) {
+ Wi[0] -= dWi[0];
+ }
+ if (Wi[4] < 0.0f) {
+ Wi[4] -= dWi[4];
+ }
+ if (Wj[0] < 0.0f) {
+ Wj[0] -= dWj[0];
+ }
+ if (Wj[4] < 0.0f) {
+ Wj[4] -= dWj[4];
+ }
}
#endif // SWIFT_HYDRO_GRADIENTS_H
diff --git a/src/hydro/Shadowswift/hydro_gradients_shadowfax.h b/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
index a8a62e076393359bf35b42315dff82b4317a2a49..4a7df6560f8409a279e349c3226dcbc60b554b05 100644
--- a/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
+++ b/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
@@ -185,7 +185,7 @@ hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
__attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
struct part *p) {
- float volume = p->geometry.volume;
+ float volume = p->cell.volume;
p->primitives.gradients.rho[0] /= volume;
p->primitives.gradients.rho[1] /= volume;
diff --git a/src/hydro/Shadowswift/hydro_iact.h b/src/hydro/Shadowswift/hydro_iact.h
index ee3045f80db7e4efe53c735f299c305001c14b7c..40d5da95958fd2109571e083b5c5ecb53de110b9 100644
--- a/src/hydro/Shadowswift/hydro_iact.h
+++ b/src/hydro/Shadowswift/hydro_iact.h
@@ -2,7 +2,7 @@
* This file is part of SWIFT.
* Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
* Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- * Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ * 2016 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
diff --git a/src/hydro/Shadowswift/hydro_io.h b/src/hydro/Shadowswift/hydro_io.h
index aaae075fe8b8614333a95032424cda746bc52313..b3e01bf1206cbf8ccdbed9183a01b754e560dd19 100644
--- a/src/hydro/Shadowswift/hydro_io.h
+++ b/src/hydro/Shadowswift/hydro_io.h
@@ -63,6 +63,16 @@ float convert_A(struct engine* e, struct part* p) {
return p->primitives.P / pow_gamma(p->primitives.rho);
}
+float convert_Etot(struct engine* e, struct part* p) {
+ float momentum2;
+
+ momentum2 = p->conserved.momentum[0] * p->conserved.momentum[0] +
+ p->conserved.momentum[1] * p->conserved.momentum[1] +
+ p->conserved.momentum[2] * p->conserved.momentum[2];
+
+ return p->conserved.energy + 0.5f * momentum2 / p->conserved.mass;
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
@@ -94,15 +104,16 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
list[7] = 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);
+ cell.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(
"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] = io_make_output_field("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
- parts, conserved.energy);
+ list[12] =
+ io_make_output_field_convert_part("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
+ parts, conserved.energy, convert_Etot);
}
/**
diff --git a/src/hydro/Shadowswift/hydro_part.h b/src/hydro/Shadowswift/hydro_part.h
index 73c89811f740839e432562db2a611dc440b155ea..f0a360c3e46c55f88f49285fa810d8fc562fa538 100644
--- a/src/hydro/Shadowswift/hydro_part.h
+++ b/src/hydro/Shadowswift/hydro_part.h
@@ -129,18 +129,6 @@ struct part {
} conserved;
- /* Geometrical quantities used for hydro. */
- struct {
-
- /* Volume of the particle. */
- float volume;
-
- /* Geometrical shear matrix used to calculate second order accurate
- gradients */
- float matrix_E[3][3];
-
- } geometry;
-
/* Variables used for timestep calculation (currently not used). */
struct {
@@ -152,15 +140,9 @@ struct part {
/* Quantities used during the volume (=density) loop. */
struct {
- /* Particle velocity divergence. */
- float div_v;
-
/* Derivative of particle number density. */
float wcount_dh;
- /* Particle velocity curl. */
- float curl_v[3];
-
/* Particle number density. */
float wcount;
@@ -192,9 +174,6 @@ struct part {
/* Variables needed for the code to compile (should be removed/replaced). */
float rho;
- /* Old internal energy flux */
- float du_dt;
-
/* Voronoi cell. */
struct voronoi_cell cell;
diff --git a/src/hydro/Shadowswift/hydro_slope_limiters_cell.h b/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
index aa99b43721f669f47a7888a5da0b1933ca1ebd62..2e62bce4fc49c97593d73b71ab1c9c5df2353783 100644
--- a/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
+++ b/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
@@ -80,94 +80,51 @@ hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
pi->primitives.limiter.maxr = fmax(r, pi->primitives.limiter.maxr);
}
-/**
- * @brief Slope limit cell gradients
- *
- * @param p Particle.
- */
-__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
- struct part* p) {
+__attribute__((always_inline)) INLINE static void
+hydro_slope_limit_cell_quantity(float* grad, float qval, float qmin, float qmax,
+ float maxr) {
- float gradrho[3], gradv[3][3], gradP[3];
float gradtrue, gradmax, gradmin, alpha;
- gradrho[0] = p->primitives.gradients.rho[0];
- gradrho[1] = p->primitives.gradients.rho[1];
- gradrho[2] = p->primitives.gradients.rho[2];
-
- gradv[0][0] = p->primitives.gradients.v[0][0];
- gradv[0][1] = p->primitives.gradients.v[0][1];
- gradv[0][2] = p->primitives.gradients.v[0][2];
-
- gradv[1][0] = p->primitives.gradients.v[1][0];
- gradv[1][1] = p->primitives.gradients.v[1][1];
- gradv[1][2] = p->primitives.gradients.v[1][2];
-
- gradv[2][0] = p->primitives.gradients.v[2][0];
- gradv[2][1] = p->primitives.gradients.v[2][1];
- gradv[2][2] = p->primitives.gradients.v[2][2];
-
- gradP[0] = p->primitives.gradients.P[0];
- gradP[1] = p->primitives.gradients.P[1];
- gradP[2] = p->primitives.gradients.P[2];
-
- gradtrue = sqrtf(gradrho[0] * gradrho[0] + gradrho[1] * gradrho[1] +
- gradrho[2] * gradrho[2]);
+ gradtrue = sqrtf(grad[0] * grad[0] + grad[1] * grad[1] + grad[2] * grad[2]);
if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.rho[1] - p->primitives.rho;
- gradmin = p->primitives.rho - p->primitives.limiter.rho[0];
+ gradtrue *= maxr;
+ gradmax = qmax - qval;
+ gradmin = qval - qmin;
alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
- p->primitives.gradients.rho[0] *= alpha;
- p->primitives.gradients.rho[1] *= alpha;
- p->primitives.gradients.rho[2] *= alpha;
- }
-
- gradtrue = sqrtf(gradv[0][0] * gradv[0][0] + gradv[0][1] * gradv[0][1] +
- gradv[0][2] * gradv[0][2]);
- if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.v[0][1] - p->primitives.v[0];
- gradmin = p->primitives.v[0] - p->primitives.limiter.v[0][0];
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
- p->primitives.gradients.v[0][0] *= alpha;
- p->primitives.gradients.v[0][1] *= alpha;
- p->primitives.gradients.v[0][2] *= alpha;
- }
-
- gradtrue = sqrtf(gradv[1][0] * gradv[1][0] + gradv[1][1] * gradv[1][1] +
- gradv[1][2] * gradv[1][2]);
- if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.v[1][1] - p->primitives.v[1];
- gradmin = p->primitives.v[1] - p->primitives.limiter.v[1][0];
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
- p->primitives.gradients.v[1][0] *= alpha;
- p->primitives.gradients.v[1][1] *= alpha;
- p->primitives.gradients.v[1][2] *= alpha;
+ grad[0] *= alpha;
+ grad[1] *= alpha;
+ grad[2] *= alpha;
}
+}
- gradtrue = sqrtf(gradv[2][0] * gradv[2][0] + gradv[2][1] * gradv[2][1] +
- gradv[2][2] * gradv[2][2]);
- if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.v[2][1] - p->primitives.v[2];
- gradmin = p->primitives.v[2] - p->primitives.limiter.v[2][0];
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
- p->primitives.gradients.v[2][0] *= alpha;
- p->primitives.gradients.v[2][1] *= alpha;
- p->primitives.gradients.v[2][2] *= alpha;
- }
+/**
+ * @brief Slope limit cell gradients
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
+ struct part* p) {
- gradtrue =
- sqrtf(gradP[0] * gradP[0] + gradP[1] * gradP[1] + gradP[2] * gradP[2]);
- if (gradtrue) {
- gradtrue *= p->primitives.limiter.maxr;
- gradmax = p->primitives.limiter.P[1] - p->primitives.P;
- gradmin = p->primitives.P - p->primitives.limiter.P[0];
- alpha = fmin(1.0f, fmin(gradmax / gradtrue, gradmin / gradtrue));
- p->primitives.gradients.P[0] *= alpha;
- p->primitives.gradients.P[1] *= alpha;
- p->primitives.gradients.P[2] *= alpha;
- }
+ hydro_slope_limit_cell_quantity(
+ p->primitives.gradients.rho, p->primitives.rho,
+ p->primitives.limiter.rho[0], p->primitives.limiter.rho[1],
+ p->primitives.limiter.maxr);
+
+ hydro_slope_limit_cell_quantity(
+ p->primitives.gradients.v[0], p->primitives.v[0],
+ p->primitives.limiter.v[0][0], p->primitives.limiter.v[0][1],
+ p->primitives.limiter.maxr);
+ hydro_slope_limit_cell_quantity(
+ p->primitives.gradients.v[1], p->primitives.v[1],
+ p->primitives.limiter.v[1][0], p->primitives.limiter.v[1][1],
+ p->primitives.limiter.maxr);
+ hydro_slope_limit_cell_quantity(
+ p->primitives.gradients.v[2], p->primitives.v[2],
+ p->primitives.limiter.v[2][0], p->primitives.limiter.v[2][1],
+ p->primitives.limiter.maxr);
+
+ hydro_slope_limit_cell_quantity(
+ p->primitives.gradients.P, p->primitives.P, p->primitives.limiter.P[0],
+ p->primitives.limiter.P[1], p->primitives.limiter.maxr);
}
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 7e8069f6ad3251420f92daedf0e209561f3507b2..fed0fb133ff56497d1a142e4f76cd1adbb97f3d1 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -197,7 +197,12 @@ void reset_particles(struct cell *c, enum velocity_field vel,
set_energy_state(p, press, size, density);
#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
- p->geometry.volume = p->mass / density;
+ float volume = p->mass / density;
+#if defined(GIZMO_SPH)
+ p->geometry.volume = volume;
+#else
+ p->cell.volume = volume;
+#endif
p->primitives.rho = density;
p->primitives.v[0] = p->v[0];
p->primitives.v[1] = p->v[1];
@@ -207,7 +212,7 @@ void reset_particles(struct cell *c, enum velocity_field vel,
p->conserved.momentum[1] = p->conserved.mass * p->v[1];
p->conserved.momentum[2] = p->conserved.mass * p->v[2];
p->conserved.energy =
- p->primitives.P / hydro_gamma_minus_one * p->geometry.volume +
+ p->primitives.P / hydro_gamma_minus_one * volume +
0.5f * (p->conserved.momentum[0] * p->conserved.momentum[0] +
p->conserved.momentum[1] * p->conserved.momentum[1] +
p->conserved.momentum[2] * p->conserved.momentum[2]) /
@@ -270,7 +275,12 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
hydro_first_init_part(part, xpart);
#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
- part->geometry.volume = part->mass / density;
+ float volume = part->mass / density;
+#ifdef GIZMO_SPH
+ part->geometry.volume = volume;
+#else
+ part->cell.volume = volume;
+#endif
part->primitives.rho = density;
part->primitives.v[0] = part->v[0];
part->primitives.v[1] = part->v[1];
@@ -280,7 +290,7 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
part->conserved.momentum[1] = part->conserved.mass * part->v[1];
part->conserved.momentum[2] = part->conserved.mass * part->v[2];
part->conserved.energy =
- part->primitives.P / hydro_gamma_minus_one * part->geometry.volume +
+ part->primitives.P / hydro_gamma_minus_one * volume +
0.5f * (part->conserved.momentum[0] * part->conserved.momentum[0] +
part->conserved.momentum[1] * part->conserved.momentum[1] +
part->conserved.momentum[2] * part->conserved.momentum[2]) /
@@ -353,7 +363,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
main_cell->parts[pid].v[2], main_cell->parts[pid].h,
main_cell->parts[pid].rho,
-#ifdef MINIMAL_SPH
+#if defined(MINIMAL_SPH) || defined(SHADOWSWIFT)
0.f,
#else
main_cell->parts[pid].density.div_v,