From 1aeb605a585176f4f84fe417ff5d79b6e25f993c Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Tue, 26 Jul 2016 17:25:37 +0100
Subject: [PATCH] Tests now also succesful for DEFAULT_SPH
---
src/const.h | 4 ++--
src/hydro/Default/hydro.h | 42 ++++++++++++++++++++++++++++++++--
src/hydro/Default/hydro_iact.h | 28 +++++++++++------------
src/hydro/Default/hydro_part.h | 2 +-
tests/test125cells.c | 13 ++++++++++-
tests/testNonSymInt.c | 10 +-------
tests/testSPHStep.c | 14 ++++++++----
tests/testSingle.c | 4 ++--
tests/testSymInt.c | 10 +-------
9 files changed, 82 insertions(+), 45 deletions(-)
diff --git a/src/const.h b/src/const.h
index 50d72e3339..270474d05e 100644
--- a/src/const.h
+++ b/src/const.h
@@ -56,8 +56,8 @@
/* SPH variant to use */
//#define MINIMAL_SPH
-#define GADGET2_SPH
-//#define DEFAULT_SPH
+//#define GADGET2_SPH
+#define DEFAULT_SPH
/* External potential properties */
#define EXTERNAL_POTENTIAL_POINTMASS
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 6e6e4f18d2..f1b414be2b 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -20,6 +20,8 @@
#include "adiabatic_index.h"
#include "approx_math.h"
+#include <float.h>
+
/**
* @brief Computes the hydro time-step of a given particle
*
@@ -145,7 +147,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute this particle's sound speed. */
const float u = p->u;
- const float fc = p->force.c = sqrtf(hydro_gamma * hydro_gamma_minus_one * u);
+ const float fc = p->force.soundspeed = sqrtf(hydro_gamma * hydro_gamma_minus_one * u);
/* Compute the P/Omega/rho2. */
xp->omega = 1.0f + 0.3333333333f * h * p->rho_dh / p->rho;
@@ -155,7 +157,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->force.balsara = normDiv_v / (normDiv_v + normRot_v + 0.0001f * fc * ih);
/* Viscosity parameter decay time */
- const float tau = h / (2.f * const_viscosity_length * p->force.c);
+ const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
/* Viscosity source term */
const float S = fmaxf(-normDiv_v, 0.f);
@@ -262,3 +264,39 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
return p->u;
}
+
+/**
+ * @brief Returns the pressure of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+ const struct part *restrict p, float dt) {
+
+ return p->force.POrho2 * p->rho * p->rho / p->rho_dh;
+}
+
+/**
+ * @brief Returns the entropy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+ const struct part *restrict p, float dt) {
+
+ return hydro_gamma_minus_one * p->u * pow_minus_gamma_minus_one(p->rho);
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
+ const struct part *restrict p, float dt) {
+
+ return p->force.soundspeed;
+}
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index 41bf3d6fbc..0e1616f889 100644
--- a/src/hydro/Default/hydro_iact.h
+++ b/src/hydro/Default/hydro_iact.h
@@ -398,7 +398,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
omega_ij = fminf(dvdr, 0.f);
/* Compute signal velocity */
- v_sig = pi->force.c + pj->force.c - 2.0f * omega_ij;
+ v_sig = pi->force.soundspeed + pj->force.soundspeed - 2.0f * omega_ij;
/* Compute viscosity parameter */
alpha_ij = -0.5f * (pi->alpha + pj->alpha);
@@ -496,11 +496,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
pju.v = vec_set(pj[0]->u, pj[1]->u, pj[2]->u, pj[3]->u, pj[4]->u, pj[5]->u,
pj[6]->u, pj[7]->u);
ci.v =
- vec_set(pi[0]->force.c, pi[1]->force.c, pi[2]->force.c, pi[3]->force.c,
- pi[4]->force.c, pi[5]->force.c, pi[6]->force.c, pi[7]->force.c);
+ vec_set(pi[0]->force.soundspeed, pi[1]->force.soundspeed, pi[2]->force.soundspeed, pi[3]->force.soundspeed,
+ pi[4]->force.soundspeed, pi[5]->force.soundspeed, pi[6]->force.soundspeed, pi[7]->force.soundspeed);
cj.v =
- vec_set(pj[0]->force.c, pj[1]->force.c, pj[2]->force.c, pj[3]->force.c,
- pj[4]->force.c, pj[5]->force.c, pj[6]->force.c, pj[7]->force.c);
+ vec_set(pj[0]->force.soundspeed, pj[1]->force.soundspeed, pj[2]->force.soundspeed, pj[3]->force.soundspeed,
+ pj[4]->force.soundspeed, pj[5]->force.soundspeed, pj[6]->force.soundspeed, pj[7]->force.soundspeed);
vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
pi[3]->force.v_sig, pi[4]->force.v_sig, pi[5]->force.v_sig,
pi[6]->force.v_sig, pi[7]->force.v_sig);
@@ -539,9 +539,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
piu.v = vec_set(pi[0]->u, pi[1]->u, pi[2]->u, pi[3]->u);
pju.v = vec_set(pj[0]->u, pj[1]->u, pj[2]->u, pj[3]->u);
ci.v =
- vec_set(pi[0]->force.c, pi[1]->force.c, pi[2]->force.c, pi[3]->force.c);
+ vec_set(pi[0]->force.soundspeed, pi[1]->force.soundspeed, pi[2]->force.soundspeed, pi[3]->force.soundspeed);
cj.v =
- vec_set(pj[0]->force.c, pj[1]->force.c, pj[2]->force.c, pj[3]->force.c);
+ vec_set(pj[0]->force.soundspeed, pj[1]->force.soundspeed, pj[2]->force.soundspeed, pj[3]->force.soundspeed);
vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
pi[3]->force.v_sig);
vj_sig.v = vec_set(pj[0]->force.v_sig, pj[1]->force.v_sig, pj[2]->force.v_sig,
@@ -711,7 +711,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
omega_ij = fminf(dvdr, 0.f);
/* Compute signal velocity */
- v_sig = pi->force.c + pj->force.c - 2.0f * omega_ij;
+ v_sig = pi->force.soundspeed + pj->force.soundspeed - 2.0f * omega_ij;
/* Compute viscosity parameter */
alpha_ij = -0.5f * (pi->alpha + pj->alpha);
@@ -801,11 +801,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
pju.v = vec_set(pj[0]->u, pj[1]->u, pj[2]->u, pj[3]->u, pj[4]->u, pj[5]->u,
pj[6]->u, pj[7]->u);
ci.v =
- vec_set(pi[0]->force.c, pi[1]->force.c, pi[2]->force.c, pi[3]->force.c,
- pi[4]->force.c, pi[5]->force.c, pi[6]->force.c, pi[7]->force.c);
+ vec_set(pi[0]->force.soundspeed, pi[1]->force.soundspeed, pi[2]->force.soundspeed, pi[3]->force.soundspeed,
+ pi[4]->force.soundspeed, pi[5]->force.soundspeed, pi[6]->force.soundspeed, pi[7]->force.soundspeed);
cj.v =
- vec_set(pj[0]->force.c, pj[1]->force.c, pj[2]->force.c, pj[3]->force.c,
- pj[4]->force.c, pj[5]->force.c, pj[6]->force.c, pj[7]->force.c);
+ vec_set(pj[0]->force.soundspeed, pj[1]->force.soundspeed, pj[2]->force.soundspeed, pj[3]->force.soundspeed,
+ pj[4]->force.soundspeed, pj[5]->force.soundspeed, pj[6]->force.soundspeed, pj[7]->force.soundspeed);
vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
pi[3]->force.v_sig, pi[4]->force.v_sig, pi[5]->force.v_sig,
pi[6]->force.v_sig, pi[7]->force.v_sig);
@@ -843,9 +843,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
piu.v = vec_set(pi[0]->u, pi[1]->u, pi[2]->u, pi[3]->u);
pju.v = vec_set(pj[0]->u, pj[1]->u, pj[2]->u, pj[3]->u);
ci.v =
- vec_set(pi[0]->force.c, pi[1]->force.c, pi[2]->force.c, pi[3]->force.c);
+ vec_set(pi[0]->force.soundspeed, pi[1]->force.soundspeed, pi[2]->force.soundspeed, pi[3]->force.soundspeed);
cj.v =
- vec_set(pj[0]->force.c, pj[1]->force.c, pj[2]->force.c, pj[3]->force.c);
+ vec_set(pj[0]->force.soundspeed, pj[1]->force.soundspeed, pj[2]->force.soundspeed, pj[3]->force.soundspeed);
vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
pi[3]->force.v_sig);
vj_sig.v = vec_set(pj[0]->force.v_sig, pj[1]->force.v_sig, pj[2]->force.v_sig,
diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h
index 7730133fa4..25ffb18053 100644
--- a/src/hydro/Default/hydro_part.h
+++ b/src/hydro/Default/hydro_part.h
@@ -99,7 +99,7 @@ struct part {
float v_sig;
/* Sound speed */
- float c;
+ float soundspeed;
/* Change in smoothing length over time. */
float h_dt;
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 5bc3d85ca0..a5fca13bec 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -95,6 +95,8 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
#if defined(GADGET2_SPH)
part->entropy = pressure / pow_gamma(density);
+#elif defined(DEFAULT_SPH)
+ part->u = pressure / (hydro_gamma_minus_one * density);
#else
error("Need to define pressure here !");
#endif
@@ -188,6 +190,8 @@ void get_solution(const struct cell *main_cell, struct solution_part *solution,
* separation.
* @param density The density of the fluid.
* @param partId The running counter of IDs.
+ * @param vel The type of velocity field.
+ * @param press The type of pressure field.
*/
struct cell *make_cell(size_t n, const double offset[3], double size, double h,
double density, long long *partId,
@@ -302,7 +306,14 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->parts[pid].a_hydro[0], main_cell->parts[pid].a_hydro[1],
main_cell->parts[pid].a_hydro[2], main_cell->parts[pid].force.h_dt,
main_cell->parts[pid].force.v_sig,
- main_cell->parts[pid].entropy_dt, 0.f);
+#if defined(GADGET2_SPH)
+ main_cell->parts[pid].entropy_dt, 0.f
+#elif defined(DEFAULT_SPH)
+ 0.f, main_cell->parts[pid].force.u_dt
+#else
+ 0.f, 0.f
+#endif
+ );
}
if (with_solution) {
diff --git a/tests/testNonSymInt.c b/tests/testNonSymInt.c
index da98c9efce..53f8510644 100644
--- a/tests/testNonSymInt.c
+++ b/tests/testNonSymInt.c
@@ -79,16 +79,8 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
"%13e %13e %13e\n",
p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->rho,
p->rho_dh, p->density.wcount, p->density.wcount_dh,
-#if defined(GADGET2_SPH)
p->density.div_v, p->density.rot_v[0], p->density.rot_v[1],
- p->density.rot_v[2]
-#elif defined(DEFAULT_SPH)
- p->density.div_v, p->density.curl_v[0], p->density.curl_v[1],
- p->density.curl_v[2]
-#else
- 0., 0., 0., 0.
-#endif
- );
+ p->density.rot_v[2]);
fclose(file);
}
diff --git a/tests/testSPHStep.c b/tests/testSPHStep.c
index dade82f1b4..ebc517f478 100644
--- a/tests/testSPHStep.c
+++ b/tests/testSPHStep.c
@@ -61,9 +61,9 @@ struct cell *make_cell(size_t N, float cellSize, int offset[3], int id_offset) {
offset[2] * cellSize + z * cellSize / N + cellSize / (2 * N);
part->h = h;
part->id = x * N * N + y * N + z + id_offset;
- ++part;
part->ti_begin = 0;
part->ti_end = 1;
+ ++part;
}
}
}
@@ -193,16 +193,20 @@ int main() {
runner_do_ghost(&r, ci);
message("h=%f rho=%f N_ngb=%f", p->h, p->rho, p->density.wcount);
- message("c=%f", p->force.c);
+ message("c=%f", p->force.soundspeed);
runner_doself2_force(&r, ci);
runner_do_kick(&r, ci, 1);
message("ti_end=%d", p->ti_end);
- free(ci->parts);
- free(ci->xparts);
-
+ for (int j = 0; j < 27; ++j) {
+ free(cells[j]->parts);
+ free(cells[j]->xparts);
+ free(cells[j]->sort);
+ free(cells[j]);
+ }
+
return 0;
#endif
}
diff --git a/tests/testSingle.c b/tests/testSingle.c
index 6f6a8c8ad6..79e1630b35 100644
--- a/tests/testSingle.c
+++ b/tests/testSingle.c
@@ -81,9 +81,9 @@ int main(int argc, char *argv[]) {
p2.rho = 1.0f;
p2.mass = 9.7059e-4;
p2.h = 0.222871287 / 2;
- p1.force.c = 0.0040824829f;
+ p1.force.soundspeed = 0.0040824829f;
p1.force.balsara = 0.0f;
- p2.force.c = 58.8972740361f;
+ p2.force.soundspeed = 58.8972740361f;
p2.force.balsara = 0.0f;
p1.u = 1.e-5 / (hydro_gamma_minus_one * p1.rho);
p2.u = 1.e-5 / (hydro_gamma_minus_one * p2.rho) + 100.0f / (33 * p2.mass);
diff --git a/tests/testSymInt.c b/tests/testSymInt.c
index 73fa8ed9b9..15d497515b 100644
--- a/tests/testSymInt.c
+++ b/tests/testSymInt.c
@@ -76,16 +76,8 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
"%13e %13e %13e\n",
p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->rho,
p->rho_dh, p->density.wcount, p->density.wcount_dh,
-#if defined(GADGET2_SPH)
p->density.div_v, p->density.rot_v[0], p->density.rot_v[1],
- p->density.rot_v[2]
-#elif defined(DEFAULT_SPH)
- p->density.div_v, p->density.curl_v[0], p->density.curl_v[1],
- p->density.curl_v[2]
-#else
- 0., 0., 0., 0.
-#endif
- );
+ p->density.rot_v[2]);
fclose(file);
}
--
GitLab