From 2f5068ee203f4a1c82a99393c8cf8932b3afbe18 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Mon, 25 Jul 2016 22:15:12 +0100
Subject: [PATCH] test_125 now checks all particles against the similar
brute-force call
---
src/hydro/Gadget2/hydro_part.h | 52 +--
src/tools.c | 100 +++++
src/tools.h | 2 +
tests/difffloat.py | 1 +
tests/test125cells.c | 421 ++++++++++++++++------
tests/test125cells.sh.in | 17 +-
tests/test27cells.c | 18 +-
tests/test27cells.sh.in | 11 +-
tests/test27cellsPerturbed.sh.in | 4 +-
tests/testPair.sh.in | 2 +-
tests/testPairPerturbed.sh.in | 2 +-
tests/tolerance_125.dat | 3 +
tests/{tolerance.dat => tolerance_27.dat} | 0
tests/tolerance_pair.dat | 3 +
14 files changed, 471 insertions(+), 165 deletions(-)
create mode 100644 tests/tolerance_125.dat
rename tests/{tolerance.dat => tolerance_27.dat} (100%)
create mode 100644 tests/tolerance_pair.dat
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index b2e5e53094..afca176400 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -61,46 +61,46 @@ struct part {
/* Derivative of the density with respect to smoothing length. */
float rho_dh;
- // union {
+ union {
- struct {
+ struct {
- /* Number of neighbours. */
- float wcount;
+ /* Number of neighbours. */
+ float wcount;
- /* Number of neighbours spatial derivative. */
- float wcount_dh;
+ /* Number of neighbours spatial derivative. */
+ float wcount_dh;
- /* Particle velocity curl. */
- float rot_v[3];
+ /* Particle velocity curl. */
+ float rot_v[3];
- /* Particle velocity divergence. */
- float div_v;
+ /* Particle velocity divergence. */
+ float div_v;
- } density;
+ } density;
- struct {
+ struct {
- /* Balsara switch */
- float balsara;
+ /* Balsara switch */
+ float balsara;
- /* Pressure over density squared (including drho/dh term) */
- float P_over_rho2;
+ /* Pressure over density squared (including drho/dh term) */
+ float P_over_rho2;
- /* Particle sound speed. */
- float soundspeed;
+ /* Particle sound speed. */
+ float soundspeed;
- /* Signal velocity. */
- float v_sig;
+ /* Signal velocity. */
+ float v_sig;
- /* Entropy time derivative */
- float entropy_dt;
+ /* Entropy time derivative */
+ float entropy_dt;
- /* Time derivative of the smoothing length */
- float h_dt;
+ /* Time derivative of the smoothing length */
+ float h_dt;
- } force;
- // };
+ } force;
+ };
/* Particle ID. */
long long id;
diff --git a/src/tools.c b/src/tools.c
index 710efd7e4d..50ad48af38 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -240,6 +240,70 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
}
}
+void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
+
+ float r2, hi, hj, hig2, hjg2, dx[3];
+ struct part *pi, *pj;
+
+ /* Implements a double-for loop and checks every interaction */
+ for (int i = 0; i < ci->count; ++i) {
+
+ pi = &ci->parts[i];
+ hi = pi->h;
+ hig2 = hi * hi * kernel_gamma2;
+
+ for (int j = 0; j < cj->count; ++j) {
+
+ pj = &cj->parts[j];
+ hj = pj->h;
+ hjg2 = hj * hj * kernel_gamma2;
+
+ /* Pairwise distance */
+ r2 = 0.0f;
+ for (int k = 0; k < 3; k++) {
+ dx[k] = ci->parts[i].x[k] - cj->parts[j].x[k];
+ r2 += dx[k] * dx[k];
+ }
+
+ /* Hit or miss? */
+ if (r2 < hig2 || r2 < hjg2) {
+
+ /* Interact */
+ runner_iact_nonsym_force(r2, dx, hi, hj, pi, pj);
+ }
+ }
+ }
+
+ /* Reverse double-for loop and checks every interaction */
+ for (int j = 0; j < cj->count; ++j) {
+
+ pj = &cj->parts[j];
+ hj = pj->h;
+ hjg2 = hj * hj * kernel_gamma2;
+
+ for (int i = 0; i < ci->count; ++i) {
+
+ pi = &ci->parts[i];
+ hi = pi->h;
+ hig2 = hi * hi * kernel_gamma2;
+
+ /* Pairwise distance */
+ r2 = 0.0f;
+ for (int k = 0; k < 3; k++) {
+ dx[k] = cj->parts[j].x[k] - ci->parts[i].x[k];
+ r2 += dx[k] * dx[k];
+ }
+
+ /* Hit or miss? */
+ if (r2 < hjg2 || r2 < hig2) {
+
+ /* Interact */
+ runner_iact_nonsym_force(r2, dx, hj, pi->h, pj, pi);
+ }
+ }
+ }
+}
+
void self_all_density(struct runner *r, struct cell *ci) {
float r2, hi, hj, hig2, hjg2, dxi[3]; //, dxj[3];
struct part *pi, *pj;
@@ -287,6 +351,42 @@ void self_all_density(struct runner *r, struct cell *ci) {
}
}
+void self_all_force(struct runner *r, struct cell *ci) {
+ float r2, hi, hj, hig2, hjg2, dxi[3]; //, dxj[3];
+ struct part *pi, *pj;
+
+ /* Implements a double-for loop and checks every interaction */
+ for (int i = 0; i < ci->count; ++i) {
+
+ pi = &ci->parts[i];
+ hi = pi->h;
+ hig2 = hi * hi * kernel_gamma2;
+
+ for (int j = i + 1; j < ci->count; ++j) {
+
+ pj = &ci->parts[j];
+ hj = pj->h;
+ hjg2 = hj * hj * kernel_gamma2;
+
+ if (pi == pj) continue;
+
+ /* Pairwise distance */
+ r2 = 0.0f;
+ for (int k = 0; k < 3; k++) {
+ dxi[k] = ci->parts[i].x[k] - ci->parts[j].x[k];
+ r2 += dxi[k] * dxi[k];
+ }
+
+ /* Hit or miss? */
+ if (r2 < hig2 || r2 < hjg2) {
+
+ /* Interact */
+ runner_iact_force(r2, dxi, hi, hj, pi, pj);
+ }
+ }
+ }
+}
+
void pairs_single_grav(double *dim, long long int pid,
struct gpart *restrict gparts, const struct part *parts,
int N, int periodic) {
diff --git a/src/tools.h b/src/tools.h
index 97f036994b..ba5b2f8c7f 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -40,6 +40,8 @@ void pairs_single_density(double *dim, long long int pid,
void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj);
void self_all_density(struct runner *r, struct cell *ci);
+void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj);
+void self_all_force(struct runner *r, struct cell *ci);
void pairs_n2(double *dim, struct part *restrict parts, int N, int periodic);
diff --git a/tests/difffloat.py b/tests/difffloat.py
index d4b48d54cb..7b0a27f668 100644
--- a/tests/difffloat.py
+++ b/tests/difffloat.py
@@ -89,6 +89,7 @@ for i in range(n_lines_to_check):
abs_diff = abs(data1[i,j] - data2[i,j])
sum = abs(data1[i,j] + data2[i,j])
+ if sum < 1e8: continue
if sum > 0:
rel_diff = abs(data1[i,j] - data2[i,j]) / sum
else:
diff --git a/tests/test125cells.c b/tests/test125cells.c
index bb59df5a8e..3d96c6d370 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -30,6 +30,153 @@
/* Local headers. */
#include "swift.h"
+enum velocity_field {
+ velocity_zero,
+ velocity_const,
+ velocity_divergent,
+ velocity_rotating
+};
+
+enum pressure_field { pressure_const, pressure_gradient, pressure_divergent };
+
+void set_velocity(struct part *part, enum velocity_field vel, float size) {
+
+ switch (vel) {
+ case velocity_zero:
+ part->v[0] = 0.f;
+ part->v[1] = 0.f;
+ part->v[2] = 0.f;
+ break;
+ case velocity_const:
+ part->v[0] = 1.f;
+ part->v[1] = 0.f;
+ part->v[2] = 0.f;
+ break;
+ case velocity_divergent:
+ part->v[0] = part->x[0] - 2.5 * size;
+ part->v[1] = part->x[1] - 2.5 * size;
+ part->v[2] = part->x[2] - 2.5 * size;
+ break;
+ case velocity_rotating:
+ part->v[0] = part->x[1];
+ part->v[1] = -part->x[0];
+ part->v[2] = 0.f;
+ break;
+ }
+}
+
+float get_pressure(double x[3], enum pressure_field press, float size) {
+
+ float r2 = 0.;
+ float dx[3] = {0.f};
+
+ switch (press) {
+ case pressure_const:
+ return 1.5f;
+ break;
+ case pressure_gradient:
+ return 1.5f * x[0]; /* gradient along x */
+ break;
+ case pressure_divergent:
+ dx[0] = x[0] - 2.5 * size;
+ dx[1] = x[1] - 2.5 * size;
+ dx[2] = x[2] - 2.5 * size;
+ r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+ return sqrt(r2) + 1.5f;
+ break;
+ }
+ return 0.f;
+}
+
+void set_energy_state(struct part *part, enum pressure_field press, float size,
+ float density) {
+
+ const float pressure = get_pressure(part->x, press, size);
+
+#if defined(GADGET2_SPH)
+ part->entropy = pressure / pow_gamma(density);
+#else
+ error("Need to define pressure here !");
+#endif
+}
+
+struct solution_part {
+
+ long long id;
+ double x[3];
+ float v[3];
+ float a_hydro[3];
+ float h;
+ float rho;
+ float div_v;
+ float S;
+ float u;
+ float P;
+ float c;
+ float h_dt;
+ float v_sig;
+ float S_dt;
+ float u_dt;
+};
+
+void get_solution(const struct cell *main_cell, struct solution_part *solution,
+ float density, enum velocity_field vel,
+ enum pressure_field press, float size) {
+
+ for (size_t i = 0; i < main_cell->count; ++i) {
+
+ solution[i].id = main_cell->parts[i].id;
+
+ solution[i].x[0] = main_cell->parts[i].x[0];
+ solution[i].x[1] = main_cell->parts[i].x[1];
+ solution[i].x[2] = main_cell->parts[i].x[2];
+
+ solution[i].v[0] = main_cell->parts[i].v[0];
+ solution[i].v[1] = main_cell->parts[i].v[1];
+ solution[i].v[2] = main_cell->parts[i].v[2];
+
+ solution[i].h = main_cell->parts[i].h;
+
+ solution[i].rho = density;
+
+ solution[i].P = get_pressure(solution[i].x, press, size);
+ solution[i].u = solution[i].P / (solution[i].rho * hydro_gamma_minus_one);
+ solution[i].S = solution[i].P / pow_gamma(solution[i].rho);
+ solution[i].c = sqrt(hydro_gamma * solution[i].P / solution[i].rho);
+
+ if (vel == velocity_divergent)
+ solution[i].div_v = 3.f;
+ else
+ solution[i].div_v = 0.f;
+
+ solution[i].h_dt = solution[i].h * solution[i].div_v / 3.;
+
+ float gradP[3] = {0.f};
+ if (press == pressure_gradient) {
+ gradP[0] = 1.5f;
+ gradP[1] = 0.f;
+ gradP[2] = 0.f;
+ } else if (press == pressure_divergent) {
+ float dx[3];
+ dx[0] = solution[i].x[0] - 2.5 * size;
+ dx[1] = solution[i].x[1] - 2.5 * size;
+ dx[2] = solution[i].x[2] - 2.5 * size;
+ float r = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
+ if (r > 0.) {
+ gradP[0] = dx[0] / r;
+ gradP[1] = dx[1] / r;
+ gradP[2] = dx[2] / r;
+ }
+ }
+
+ solution[i].a_hydro[0] = -gradP[0] / solution[i].rho;
+ solution[i].a_hydro[1] = -gradP[1] / solution[i].rho;
+ solution[i].a_hydro[2] = -gradP[2] / solution[i].rho;
+
+ solution[i].v_sig = 2.f * solution[i].c;
+ }
+}
+
/**
* @brief Constructs a cell and all of its particle in a valid state prior to
* a SPH time-step.
@@ -41,11 +188,10 @@
* separation.
* @param density The density of the fluid.
* @param partId The running counter of IDs.
- * @param pert The perturbation to apply to the particles in the cell in units
- * of the inter-particle separation.
*/
struct cell *make_cell(size_t n, const double offset[3], double size, double h,
- double density, long long *partId, double pert) {
+ double density, long long *partId,
+ enum velocity_field vel, enum pressure_field press) {
const size_t count = n * n * n;
const double volume = size * size * size;
@@ -67,28 +213,24 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
for (size_t x = 0; x < n; ++x) {
for (size_t y = 0; y < n; ++y) {
for (size_t z = 0; z < n; ++z) {
- part->x[0] =
- offset[0] +
- size * (x + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
- part->x[1] =
- offset[1] +
- size * (y + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
- part->x[2] =
- offset[2] +
- size * (z + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
- part->v[0] = 0.f;
- part->v[1] = 0.f;
- part->v[2] = 0.f;
+ part->x[0] = offset[0] + size * (x + 0.5) / (float)n;
+ part->x[1] = offset[1] + size * (y + 0.5) / (float)n;
+ part->x[2] = offset[2] + size * (z + 0.5) / (float)n;
part->h = size * h / (float)n;
- part->entropy = 1.f;
- part->id = ++(*partId);
part->mass = density * volume / count;
+
+ set_velocity(part, vel, size);
+ set_energy_state(part, press, size, density);
+
+ part->id = ++(*partId);
part->ti_begin = 0;
part->ti_end = 1;
+
xpart->v_full[0] = part->v[0];
xpart->v_full[1] = part->v[1];
xpart->v_full[2] = part->v[2];
++part;
+ ++xpart;
}
}
}
@@ -129,78 +271,59 @@ void clean_up(struct cell *ci) {
* @brief Dump all the particles to a file
*/
void dump_particle_fields(char *fileName, struct cell *main_cell,
- struct cell **cells) {
+ struct solution_part *solution, int with_solution) {
FILE *file = fopen(fileName, "w");
/* Write header */
fprintf(file,
- "# %4s %10s %10s %10s %10s %10s %10s %13s %13s %13s %13s %13s "
- "%13s %13s %13s\n",
- "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");
+ "# %4s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s "
+ "%8s %8s %8s %8s %8s\n",
+ "ID", "pos_x", "pos_y", "pos_z", "v_x", "v_y", "v_z", "h", "rho",
+ "div_v", "S", "u", "P", "c", "a_x", "a_y", "a_z", "h_dt", "v_sig",
+ "dS/dt", "du/dt");
fprintf(file, "# Main cell --------------------------------------------\n");
/* Write main cell */
for (size_t pid = 0; pid < main_cell->count; pid++) {
fprintf(file,
- "%6llu %10f %10f %10f %10f %10f %10f %13e %13e %13e %13e %13e "
- "%13e %13e %13e\n",
+ "%6llu %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f "
+ "%8.5f "
+ "%8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f\n",
main_cell->parts[pid].id, main_cell->parts[pid].x[0],
main_cell->parts[pid].x[1], main_cell->parts[pid].x[2],
main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
- main_cell->parts[pid].v[2], main_cell->parts[pid].rho,
- main_cell->parts[pid].rho_dh, main_cell->parts[pid].density.wcount,
- main_cell->parts[pid].density.wcount_dh,
-#if defined(GADGET2_SPH)
- main_cell->parts[pid].div_v, main_cell->parts[pid].density.rot_v[0],
- main_cell->parts[pid].density.rot_v[1],
- main_cell->parts[pid].density.rot_v[2]
-#elif defined(DEFAULT_SPH)
- main_cell->parts[pid].density.div_v,
- main_cell->parts[pid].density.rot_v[0],
- main_cell->parts[pid].density.rot_v[1],
- main_cell->parts[pid].density.rot_v[2]
-#else
- 0., 0., 0., 0.
-#endif
- );
+ main_cell->parts[pid].v[2], main_cell->parts[pid].h,
+ main_cell->parts[pid].rho, main_cell->parts[pid].density.div_v,
+ hydro_get_entropy(&main_cell->parts[pid], 0.f),
+ hydro_get_internal_energy(&main_cell->parts[pid], 0.f),
+ hydro_get_pressure(&main_cell->parts[pid], 0.f),
+ hydro_get_soundspeed(&main_cell->parts[pid], 0.f),
+ 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].force.entropy_dt, 0.f);
}
- /* Write all other cells */
- for (int i = 0; i < 5; ++i) {
- for (int j = 0; j < 5; ++j) {
- for (int k = 0; k < 5; ++k) {
- struct cell *cj = cells[i * 25 + j * 5 + k];
- if (cj == main_cell) continue;
-
- fprintf(file,
- "# Offset: [%2d %2d %2d] -----------------------------------\n",
- i - 2, j - 2, k - 2);
-
- for (size_t pjd = 0; pjd < cj->count; pjd++) {
- fprintf(
- file,
- "%6llu %10f %10f %10f %10f %10f %10f %13e %13e %13e %13e %13e "
- "%13e %13e %13e\n",
- cj->parts[pjd].id, cj->parts[pjd].x[0], cj->parts[pjd].x[1],
- cj->parts[pjd].x[2], cj->parts[pjd].v[0], cj->parts[pjd].v[1],
- cj->parts[pjd].v[2], cj->parts[pjd].rho, cj->parts[pjd].rho_dh,
- cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
-#if defined(GADGET2_SPH)
- cj->parts[pjd].div_v, cj->parts[pjd].density.rot_v[0],
- cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
-#elif defined(DEFAULT_SPH)
- cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0],
- cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
-#else
- 0., 0., 0., 0.
-#endif
- );
- }
- }
+ if (with_solution) {
+
+ fprintf(file, "# Solution ---------------------------------------------\n");
+
+ for (size_t pid = 0; pid < main_cell->count; pid++) {
+ fprintf(file,
+ "%6llu %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f "
+ "%8.5f %8.5f "
+ "%8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f\n",
+ solution[pid].id, solution[pid].x[0], solution[pid].x[1],
+ solution[pid].x[2], solution[pid].v[0], solution[pid].v[1],
+ solution[pid].v[2], solution[pid].h, solution[pid].rho,
+ solution[pid].div_v, solution[pid].S, solution[pid].u,
+ solution[pid].P, solution[pid].c, solution[pid].a_hydro[0],
+ solution[pid].a_hydro[1], solution[pid].a_hydro[2],
+ solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt, 0.f);
}
}
+
fclose(file);
}
@@ -214,10 +337,11 @@ void runner_doself2_force(struct runner *r, struct cell *ci);
int main(int argc, char *argv[]) {
size_t runs = 0, particles = 0;
- double h = 1.2348, size = 1., rho = 1.;
- double perturbation = 0.;
+ double h = 1.23485, size = 1., rho = 2.5;
char outputFileNameExtension[200] = "";
char outputFileName[200] = "";
+ enum velocity_field vel = velocity_zero;
+ enum pressure_field press = pressure_const;
/* Initialize CPU frequency, this also starts time. */
unsigned long long cpufreq = 0;
@@ -230,7 +354,7 @@ int main(int argc, char *argv[]) {
srand(0);
char c;
- while ((c = getopt(argc, argv, "m:s:h:p:r:t:d:f:v:")) != -1) {
+ while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:p:")) != -1) {
switch (c) {
case 'h':
sscanf(optarg, "%lf", &h);
@@ -238,21 +362,24 @@ int main(int argc, char *argv[]) {
case 's':
sscanf(optarg, "%lf", &size);
break;
- case 'p':
+ case 'n':
sscanf(optarg, "%zu", &particles);
break;
case 'r':
sscanf(optarg, "%zu", &runs);
break;
- case 'd':
- sscanf(optarg, "%lf", &perturbation);
- break;
case 'm':
sscanf(optarg, "%lf", &rho);
break;
case 'f':
strcpy(outputFileNameExtension, optarg);
break;
+ case 'v':
+ sscanf(optarg, "%d", (int *)&vel);
+ break;
+ case 'p':
+ sscanf(optarg, "%d", (int *)&press);
+ break;
case '?':
error("Unknown option.");
break;
@@ -261,28 +388,38 @@ int main(int argc, char *argv[]) {
if (h < 0 || particles == 0 || runs == 0) {
printf(
- "\nUsage: %s -p PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
- "\nGenerates a cell pair, filled with particles on a Cartesian grid."
- "\nThese are then interacted using runner_dopair1_density."
+ "\nUsage: %s -n PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
+ "\nGenerates 125 cells, filled with particles on a Cartesian grid."
+ "\nThese are then interacted using runner_dopair1_density() and "
+ "runner_doself1_density() followed by runner_dopair2_force() and "
+ "runner_doself2_force()"
"\n\nOptions:"
"\n-h DISTANCE=1.2348 - Smoothing length in units of <x>"
"\n-m rho - Physical density in the cell"
"\n-s size - Physical size of the cell"
- "\n-d pert - Perturbation to apply to the particles [0,1["
+ "\n-v type (0,1,2,3) - Velocity field: (zero, constant, divergent, "
"rotating)"
+ "\n-p type (0,1,2) - Pressure field: (constant, gradient divergent)"
"\n-f fileName - Part of the file name used to save the dumps\n",
argv[0]);
exit(1);
}
/* Help users... */
+ message("Adiabatic index: ga = %f", hydro_gamma);
message("Smoothing length: h = %f", h * size);
message("Kernel: %s", kernel_name);
- message("Neighbour target: N = %f",
- h * h * h * 4.0 * M_PI * kernel_gamma3 / 3.0);
+ message("Neighbour target: N = %f", h * h * h * kernel_norm);
message("Density target: rho = %f", rho);
- // message("div_v target: div = %f", vel == 2 ? 3.f : 0.f);
- // message("curl_v target: curl = [0., 0., %f]", vel == 3 ? -2.f : 0.f);
+ message("div_v target: div = %f", vel == 2 ? 3.f : 0.f);
+ message("curl_v target: curl = [0., 0., %f]", vel == 3 ? -2.f : 0.f);
+ if (press == pressure_const)
+ message("P field constant");
+ else if (press == pressure_gradient)
+ message("P field gradient");
+ else
+ message("P field divergent");
+
printf("\n");
/* Build the infrastructure */
@@ -291,10 +428,9 @@ int main(int argc, char *argv[]) {
space.h_max = h;
struct hydro_props hp;
- hp.target_neighbours = h * h * h * 4.0 * M_PI * kernel_gamma3 / 3.0;
- hp.delta_neighbours = 1;
- hp.max_smoothing_iterations =
- 1; /* We construct correct h values, 1 is enough */
+ hp.target_neighbours = h * h * h * kernel_norm;
+ hp.delta_neighbours = 1.;
+ hp.max_smoothing_iterations = 1;
struct engine engine;
engine.hydro_properties = &hp;
@@ -320,7 +456,7 @@ int main(int argc, char *argv[]) {
/* Construct it */
cells[i * 25 + j * 5 + k] =
- make_cell(particles, offset, size, h, rho, &partId, perturbation);
+ make_cell(particles, offset, size, h, rho, &partId, vel, press);
/* Store the inner cells */
if (i > 0 && i < 4 && j > 0 && j < 4 && k > 0 && k < 4) {
@@ -334,6 +470,12 @@ int main(int argc, char *argv[]) {
/* Store the main cell for future use */
main_cell = cells[62];
+ /* Construct the real solution */
+ struct solution_part *solution =
+ malloc(main_cell->count * sizeof(struct solution_part));
+ get_solution(main_cell, solution, rho, vel, press, size);
+
+ /* Start the test */
ticks time = 0;
for (size_t i = 0; i < runs; ++i) {
@@ -387,8 +529,6 @@ int main(int argc, char *argv[]) {
/* Ghost to finish everything on the central cells */
for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j]);
- message("N_ngb = %f", main_cell->parts[0].density.wcount);
-
/* Do the force calculation */
#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
@@ -415,50 +555,97 @@ int main(int argc, char *argv[]) {
time += toc - tic;
/* Dump if necessary */
- if (i % 50 == 0) {
+ if (i == 0) {
sprintf(outputFileName, "swift_dopair_125_%s.dat",
outputFileNameExtension);
- dump_particle_fields(outputFileName, main_cell, cells);
+ dump_particle_fields(outputFileName, main_cell, solution, 0);
}
}
/* Output timing */
message("SWIFT calculation took : %15lli ticks.", time / runs);
- /* Now perform a brute-force version for accuracy tests */
+ /* NOW BRUTE-FORCE CALCULATION */
+
+ const ticks tic = getticks();
+
+ /* Initialise the particles */
+ for (int j = 0; j < 125; ++j) runner_do_init(&runner, cells[j], 0);
+
+/* Do the density calculation */
+#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
+
+ /* Run all the pairs (only once !)*/
+ for (int i = 0; i < 5; i++) {
+ for (int j = 0; j < 5; j++) {
+ for (int k = 0; k < 5; k++) {
+
+ struct cell *ci = cells[i * 25 + j * 5 + k];
+
+ for (int ii = -1; ii < 2; ii++) {
+ int iii = i + ii;
+ if (iii < 0 || iii >= 5) continue;
+ iii = (iii + 5) % 5;
+ for (int jj = -1; jj < 2; jj++) {
+ int jjj = j + jj;
+ if (jjj < 0 || jjj >= 5) continue;
+ jjj = (jjj + 5) % 5;
+ for (int kk = -1; kk < 2; kk++) {
+ int kkk = k + kk;
+ if (kkk < 0 || kkk >= 5) continue;
+ kkk = (kkk + 5) % 5;
+
+ struct cell *cj = cells[iii * 25 + jjj * 5 + kkk];
+
+ if (cj > ci) pairs_all_density(&runner, ci, cj);
+ }
+ }
+ }
+ }
+ }
+ }
- /* /\* Zero the fields *\/ */
- /* for (int i = 0; i < 125; ++i) zero_particle_fields(cells[i]); */
+ /* And now the self-interaction for the central cells*/
+ for (int j = 0; j < 27; ++j) self_all_density(&runner, inner_cells[j]);
- /* const ticks tic = getticks(); */
+#endif
- /* #if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION) */
+ /* Ghost to finish everything on the central cells */
+ for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j]);
- /* /\* Run all the brute-force pairs *\/ */
- /* for (int j = 0; j < 125; ++j) */
- /* if (cells[j] != main_cell) pairs_all_density(&runner, main_cell,
- * cells[j]); */
+/* Do the force calculation */
+#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
- /* /\* And now the self-interaction *\/ */
- /* self_all_density(&runner, main_cell); */
+ /* Do the pairs (for the central 27 cells) */
+ for (int i = 1; i < 4; i++) {
+ for (int j = 1; j < 4; j++) {
+ for (int k = 1; k < 4; k++) {
- /* #endif */
+ struct cell *cj = cells[i * 25 + j * 5 + k];
+
+ if (main_cell != cj) pairs_all_force(&runner, main_cell, cj);
+ }
+ }
+ }
+
+ /* And now the self-interaction for the main cell */
+ self_all_force(&runner, main_cell);
+#endif
- /* const ticks toc = getticks(); */
+ /* Finally, give a gentle kick */
+ runner_do_kick(&runner, main_cell, 0);
- /* /\* Let's get physical ! *\/ */
- /* end_calculation(main_cell); */
+ const ticks toc = getticks();
- /* /\* Dump *\/ */
- /* sprintf(outputFileName, "brute_force_125_%s.dat",
- * outputFileNameExtension); */
- /* dump_particle_fields(outputFileName, main_cell, cells); */
+ /* Output timing */
+ message("Brute force calculation took : %15lli ticks.", toc - tic);
- /* /\* Output timing *\/ */
- /* message("Brute force calculation took : %15lli ticks.", toc - tic); */
+ sprintf(outputFileName, "brute_force_125_%s.dat", outputFileNameExtension);
+ dump_particle_fields(outputFileName, main_cell, solution, 0);
/* Clean things to make the sanitizer happy ... */
for (int i = 0; i < 125; ++i) clean_up(cells[i]);
+ free(solution);
return 0;
}
diff --git a/tests/test125cells.sh.in b/tests/test125cells.sh.in
index 838bbef973..ae8eb4ec39 100755
--- a/tests/test125cells.sh.in
+++ b/tests/test125cells.sh.in
@@ -1,8 +1,13 @@
#!/bin/bash
-rm brute_force_125_standard.dat swift_dopair_125_standard.dat
-
-./test125cells -p 4 -r 1 -d 0 -f standard
-
-#python @srcdir@/difffloat.py brute_force_125_standard.dat swift_dopair_125_standard.dat @srcdir@/tolerance.dat 4
-
+for v in {0..3}
+do
+ for p in {0..2}
+ do
+ rm brute_force_125_standard.dat swift_dopair_125_standard.dat
+ ./test125cells -n 6 -r 1 -v $v -p $p -f standard
+ python @srcdir@/difffloat.py brute_force_125_standard.dat swift_dopair_125_standard.dat @srcdir@/tolerance_125.dat 6
+ echo ""
+ done
+done
+
exit $?
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 5d9cd0e2c8..c6c0c0f10d 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -249,7 +249,7 @@ void runner_doself1_density(struct runner *r, struct cell *ci);
/* And go... */
int main(int argc, char *argv[]) {
size_t runs = 0, particles = 0;
- double h = 1.2348, size = 1., rho = 1.;
+ double h = 1.23485, size = 1., rho = 1.;
double perturbation = 0.;
char outputFileNameExtension[200] = "";
char outputFileName[200] = "";
@@ -266,7 +266,7 @@ int main(int argc, char *argv[]) {
srand(0);
char c;
- while ((c = getopt(argc, argv, "m:s:h:p:r:t:d:f:v:")) != -1) {
+ while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:")) != -1) {
switch (c) {
case 'h':
sscanf(optarg, "%lf", &h);
@@ -274,7 +274,7 @@ int main(int argc, char *argv[]) {
case 's':
sscanf(optarg, "%lf", &size);
break;
- case 'p':
+ case 'n':
sscanf(optarg, "%zu", &particles);
break;
case 'r':
@@ -300,9 +300,10 @@ int main(int argc, char *argv[]) {
if (h < 0 || particles == 0 || runs == 0) {
printf(
- "\nUsage: %s -p PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
- "\nGenerates a cell pair, filled with particles on a Cartesian grid."
- "\nThese are then interacted using runner_dopair1_density."
+ "\nUsage: %s -n PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
+ "\nGenerates 27 cells, filled with particles on a Cartesian grid."
+ "\nThese are then interacted using runner_dopair1_density() and "
+ "runner_doself1_density()."
"\n\nOptions:"
"\n-h DISTANCE=1.2348 - Smoothing length in units of <x>"
"\n-m rho - Physical density in the cell"
@@ -316,13 +317,14 @@ int main(int argc, char *argv[]) {
}
/* Help users... */
+ message("Adiabatic index: ga = %f", hydro_gamma);
message("Smoothing length: h = %f", h * size);
message("Kernel: %s", kernel_name);
- message("Neighbour target: N = %f",
- h * h * h * 4.0 * M_PI * kernel_gamma3 / 3.0);
+ message("Neighbour target: N = %f", h * h * h * kernel_norm);
message("Density target: rho = %f", rho);
message("div_v target: div = %f", vel == 2 ? 3.f : 0.f);
message("curl_v target: curl = [0., 0., %f]", vel == 3 ? -2.f : 0.f);
+
printf("\n");
/* Build the infrastructure */
diff --git a/tests/test27cells.sh.in b/tests/test27cells.sh.in
index b556576647..546b9cd27f 100755
--- a/tests/test27cells.sh.in
+++ b/tests/test27cells.sh.in
@@ -1,8 +1,11 @@
#!/bin/bash
-rm brute_force_27_standard.dat swift_dopair_27_standard.dat
-./test27cells -p 6 -r 1 -d 0 -f standard
-
-python @srcdir@/difffloat.py brute_force_27_standard.dat swift_dopair_27_standard.dat @srcdir@/tolerance.dat 6
+for v in {0..3}
+do
+ rm brute_force_27_standard.dat swift_dopair_27_standard.dat
+ ./test27cells -n 6 -r 1 -d 0 -f standard -v $v
+ python @srcdir@/difffloat.py brute_force_27_standard.dat swift_dopair_27_standard.dat @srcdir@/tolerance_27.dat 6
+ echo ""
+done
exit $?
diff --git a/tests/test27cellsPerturbed.sh.in b/tests/test27cellsPerturbed.sh.in
index 3d7574f9e9..c003b046a4 100755
--- a/tests/test27cellsPerturbed.sh.in
+++ b/tests/test27cellsPerturbed.sh.in
@@ -1,8 +1,8 @@
#!/bin/bash
rm brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat
-./test27cells -p 6 -r 1 -d 0.1 -f perturbed
+./test27cells -n 6 -r 1 -d 0.1 -f perturbed
-python @srcdir@/difffloat.py brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat @srcdir@/tolerance.dat 6
+python @srcdir@/difffloat.py brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat @srcdir@/tolerance_27.dat 6
exit $?
diff --git a/tests/testPair.sh.in b/tests/testPair.sh.in
index e719aeec4a..086db14fdd 100755
--- a/tests/testPair.sh.in
+++ b/tests/testPair.sh.in
@@ -3,6 +3,6 @@ rm brute_force_standard.dat swift_dopair_standard.dat
./testPair -p 6 -r 1 -d 0 -f standard
-python @srcdir@/difffloat.py brute_force_standard.dat swift_dopair_standard.dat @srcdir@/tolerance.dat
+python @srcdir@/difffloat.py brute_force_standard.dat swift_dopair_standard.dat @srcdir@/tolerance_pair.dat
exit $?
diff --git a/tests/testPairPerturbed.sh.in b/tests/testPairPerturbed.sh.in
index d94a040d42..f17d7559e0 100755
--- a/tests/testPairPerturbed.sh.in
+++ b/tests/testPairPerturbed.sh.in
@@ -3,6 +3,6 @@ rm brute_force_perturbed.dat swift_dopair_perturbed.dat
./testPair -p 6 -r 1 -d 0.1 -f perturbed
-python @srcdir@/difffloat.py brute_force_perturbed.dat swift_dopair_perturbed.dat @srcdir@/tolerance.dat
+python @srcdir@/difffloat.py brute_force_perturbed.dat swift_dopair_perturbed.dat @srcdir@/tolerance_pair.dat
exit $?
diff --git a/tests/tolerance_125.dat b/tests/tolerance_125.dat
new file mode 100644
index 0000000000..4870b20d96
--- /dev/null
+++ b/tests/tolerance_125.dat
@@ -0,0 +1,3 @@
+# ID pos_x pos_y pos_z v_x v_y v_z h rho div_v S u P c a_x a_y a_z h_dt v_sig dS/dt du/dt
+ 0 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5
+ 0 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5
diff --git a/tests/tolerance.dat b/tests/tolerance_27.dat
similarity index 100%
rename from tests/tolerance.dat
rename to tests/tolerance_27.dat
diff --git a/tests/tolerance_pair.dat b/tests/tolerance_pair.dat
new file mode 100644
index 0000000000..f5031c5f47
--- /dev/null
+++ b/tests/tolerance_pair.dat
@@ -0,0 +1,3 @@
+# 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 1e-5 1e-5 2e-5 3e-2 1e-5 1e-5 1e-5 1e-5
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-5 1.2e-5 1e-5 1e-2 1e-4 1e-4 1e-4 1e-4
--
GitLab