diff --git a/src/hydro/Shadowswift/hydro_io.h b/src/hydro/Shadowswift/hydro_io.h
index d192ea9b8d4b22b29ab2ae5685fa8cecd5554775..a25f6847505f5a7ee23c4e4e76bf6ad29b18a07f 100644
--- a/src/hydro/Shadowswift/hydro_io.h
+++ b/src/hydro/Shadowswift/hydro_io.h
@@ -170,7 +170,7 @@ INLINE static void hydro_write_particles(const struct part* parts,
struct io_props* list,
int* num_fields) {
- *num_fields = 13;
+ *num_fields = 14;
/* List what we want to write */
list[0] = io_make_output_field_convert_part(
@@ -226,6 +226,10 @@ INLINE static void hydro_write_particles(const struct part* parts,
list[12] = io_make_output_field_convert_part(
"TotalEnergies", FLOAT, 1, UNIT_CONV_ENERGY, -3.f * hydro_gamma_minus_one,
parts, xparts, convert_Etot, "Total (co-moving) energy of the particles");
+
+ list[13] = io_make_output_field("NumFaceInteractions", INT, 1,
+ UNIT_CONV_NO_UNITS, 0.f, parts, voronoi.nface,
+ "Number of Voronoi face interactions");
}
/**
diff --git a/src/hydro/Shadowswift/hydro_part.h b/src/hydro/Shadowswift/hydro_part.h
index 027b0e2ff72d60c3e207a51f7ee180c86b6774dd..ac079a0753892a7fae010ebab435edebfd12446d 100644
--- a/src/hydro/Shadowswift/hydro_part.h
+++ b/src/hydro/Shadowswift/hydro_part.h
@@ -223,6 +223,9 @@ struct part {
tessellation. */
int flag;
+ /* Number of times this cell interacted with other cells. */
+ int nface;
+
} voronoi;
} SWIFT_STRUCT_ALIGN;
diff --git a/src/runner_doiact_functions_hydro.h b/src/runner_doiact_functions_hydro.h
index 2d8bf4fbb958ae407b089d9cce6d69e8031ed43d..97103ef0506c79a8fa0e0dbe20251897d1393e96 100644
--- a/src/runner_doiact_functions_hydro.h
+++ b/src/runner_doiact_functions_hydro.h
@@ -1065,8 +1065,13 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
const float a = cosmo->a;
const float H = cosmo->H;
-#if defined(SHADOWFAX_SPH) && (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
- cell_shadowfax_do_pair1(e, ci, cj, sid, shift);
+#if defined(SHADOWFAX_SPH)
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ cell_shadowfax_do_pair1_density(e, ci, cj, sid, shift);
+#elif (FUNCTION_TASK_LOOP == TASK_LOOP_GRADIENT)
+#else
+ error("Using wrong pair function!");
+#endif
#endif
if (cell_is_active_hydro(ci, e)) {
@@ -1408,6 +1413,15 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
const float a = cosmo->a;
const float H = cosmo->H;
+#if defined(SHADOWFAX_SPH)
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) || \
+ (FUNCTION_TASK_LOOP == TASK_LOOP_GRADIENT)
+ error("Using wrong pair function!");
+#else
+ cell_shadowfax_do_pair2_force(e, ci, cj, sid, shift);
+#endif
+#endif
+
/* Maximal displacement since last rebuild */
const double dx_max = (ci->hydro.dx_max_sort + cj->hydro.dx_max_sort);
@@ -1986,8 +2000,13 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
TIMER_TIC;
-#if defined(SHADOWFAX_SPH) && (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
- cell_shadowfax_do_self1(e, c);
+#if defined(SHADOWFAX_SPH)
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+ cell_shadowfax_do_self1_density(e, c);
+#elif (FUNCTION_TASK_LOOP == TASK_LOOP_GRADIENT)
+#else
+ error("Using wrong self function!");
+#endif
#endif
struct part *restrict parts = c->hydro.parts;
@@ -2217,6 +2236,15 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
TIMER_TIC;
+#if defined(SHADOWFAX_SPH)
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) || \
+ (FUNCTION_TASK_LOOP == TASK_LOOP_GRADIENT)
+ error("Using wrong self function!");
+#else
+ cell_shadowfax_do_self2_force(e, c);
+#endif
+#endif
+
struct part *restrict parts = c->hydro.parts;
const int count = c->hydro.count;
diff --git a/src/shadowfax/cell_shadowfax.h b/src/shadowfax/cell_shadowfax.h
index 2f369dc1de823def81a20fa1dbd737e29cb72ebf..79762bade884929eb1eac70940249e67e36b5fda 100644
--- a/src/shadowfax/cell_shadowfax.h
+++ b/src/shadowfax/cell_shadowfax.h
@@ -4,6 +4,7 @@
#include "active.h"
#include "cell.h"
#include "delaunay.h"
+#include "hydro_shadowfax.h"
#include "voronoi.h"
__attribute__((always_inline)) INLINE static void shadowfax_flag_particle_added(
@@ -34,13 +35,15 @@ cell_malloc_delaunay_tessellation(struct cell *c,
/* Get a pointer to the ith particle. */
struct part *restrict p = &parts[pd];
p->voronoi.flag = 0;
+ p->voronoi.nface = 0;
}
c->hydro.shadowfax_enabled = 1;
}
-__attribute__((always_inline)) INLINE static void cell_shadowfax_do_self1(
- const struct engine *e, struct cell *restrict c) {
+__attribute__((always_inline)) INLINE static void
+cell_shadowfax_do_self1_density(const struct engine *e,
+ struct cell *restrict c) {
const int count = c->hydro.count;
struct part *restrict parts = c->hydro.parts;
@@ -56,8 +59,19 @@ __attribute__((always_inline)) INLINE static void cell_shadowfax_do_self1(
shadowfax_flag_particle_added(p, 0);
}
}
+}
+
+__attribute__((always_inline)) INLINE static void cell_shadowfax_do_self2_force(
+ const struct engine *e, struct cell *restrict c) {
- /* delaunay_consolidate(&c->hydro.deltess);*/
+ struct part *restrict parts = c->hydro.parts;
+
+ struct voronoi *vortess = &c->hydro.vortess;
+ for (int i = 0; i < vortess->pair_index[0]; ++i) {
+ struct voronoi_pair *pair = &vortess->pairs[0][i];
+ hydro_shadowfax_flux_exchange(&parts[pair->left], &parts[pair->right],
+ pair->midpoint, pair->surface_area);
+ }
}
__attribute__((always_inline)) INLINE static void cell_shadowfax_do_pair_naive(
@@ -240,9 +254,11 @@ __attribute__((always_inline)) INLINE static void cell_shadowfax_do_pair_subset(
}
}
-__attribute__((always_inline)) INLINE static void cell_shadowfax_do_pair1(
- const struct engine *e, struct cell *restrict ci, struct cell *restrict cj,
- int sid, const double *shift) {
+__attribute__((always_inline)) INLINE static void
+cell_shadowfax_do_pair1_density(const struct engine *e,
+ struct cell *restrict ci,
+ struct cell *restrict cj, int sid,
+ const double *shift) {
if (ci == cj) error("Interacting cell with itself!");
@@ -377,6 +393,21 @@ __attribute__((always_inline)) INLINE static void cell_shadowfax_do_pair1(
} /* Cell cj is active */
}
+__attribute__((always_inline)) INLINE static void cell_shadowfax_do_pair2_force(
+ const struct engine *e, struct cell *restrict ci, struct cell *restrict cj,
+ int sid, const double *shift) {
+
+ struct part *restrict parts_i = ci->hydro.parts;
+ struct part *restrict parts_j = cj->hydro.parts;
+
+ struct voronoi *vortess = &ci->hydro.vortess;
+ for (int i = 0; i < vortess->pair_index[1 + sid]; ++i) {
+ struct voronoi_pair *pair = &vortess->pairs[1 + sid][i];
+ hydro_shadowfax_flux_exchange(&parts_i[pair->left], &parts_j[pair->right],
+ pair->midpoint, pair->surface_area);
+ }
+}
+
__attribute__((always_inline)) INLINE static void cell_shadowfax_do_pair2(
const struct engine *e, struct cell *restrict ci, struct cell *restrict cj,
int sid, const double *shift) {
diff --git a/src/shadowfax/hydro_shadowfax.h b/src/shadowfax/hydro_shadowfax.h
new file mode 100644
index 0000000000000000000000000000000000000000..235cb122220354e2869e32e1b14e16c59fc9fdd1
--- /dev/null
+++ b/src/shadowfax/hydro_shadowfax.h
@@ -0,0 +1,11 @@
+#ifndef SWIFT_HYDRO_SHADOWFAX_H
+#define SWIFT_HYDRO_SHADOWFAX_H
+
+__attribute__((always_inline)) INLINE static void hydro_shadowfax_flux_exchange(
+ struct part *restrict pi, struct part *restrict pj, double *midpoint, double surface_area) {
+
+ ++pi->voronoi.nface;
+ ++pj->voronoi.nface;
+}
+
+#endif /* SWIFT_HYDRO_SHADOWFAX_H */
diff --git a/src/shadowfax/voronoi.h b/src/shadowfax/voronoi.h
index e9077c0812776d79fa9c52eb71ff66bbeb345b47..70d4ab5191524a916398322e8e9607cdd9bf0449 100644
--- a/src/shadowfax/voronoi.h
+++ b/src/shadowfax/voronoi.h
@@ -1,6 +1,6 @@
/*******************************************************************************
* This file is part of cVoronoi.
- * Copyright (c) 2020 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ * Copyright (c) 2020, 2021 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
@@ -68,6 +68,7 @@ struct voronoi_pair {
#ifdef VORONOI_STORE_CONNECTIONS
/*! First vertex of the interface. */
double a[2];
+
/*! Second vertex of the interface. */
double b[2];
#endif
@@ -87,6 +88,7 @@ struct voronoi_cell_new {
double centroid[2];
#ifdef VORONOI_STORE_GENERATORS
+ /*! Position of the cell generator. */
double generator[2];
#endif
@@ -109,11 +111,20 @@ struct voronoi {
/*! @brief Voronoi cells. */
struct voronoi_cell_new *cells;
- /*! @brief Number of cells (and number of generators). */
+ /*! @brief Number of cells. */
int number_of_cells;
+ /*! @brief Voronoi cell pairs. We store these per (SWIFT) cell, i.e. pairs[0]
+ * contains all pairs that are completely contained within this cell, while
+ * pairs[1] corresponds to pairs crossing the boundary between this cell and
+ * the cell with coordinates that are lower in all coordinate directions (the
+ * cell to the left, front, bottom, sid=0), and so on. */
struct voronoi_pair *pairs[27];
+
+ /*! @brief Current number of pairs per cell index. */
int pair_index[27];
+
+ /*! @brief Allocated number of pairs per cell index. */
int pair_size[27];
};
@@ -219,6 +230,7 @@ static inline double voronoi_compute_centroid_volume_triangle(
*
* @param v Voronoi grid.
* @param d Delaunay tessellation (read-only).
+ * @param parts Local cell generators (read-only).
*/
static inline void voronoi_init(struct voronoi *restrict v,
const struct delaunay *restrict d,
@@ -456,25 +468,49 @@ static inline void voronoi_check_grid(const struct voronoi *restrict v) {
printf("Total volume: %g\n", V);
}
+/**
+ * @brief Write the Voronoi grid information to the given file.
+ *
+ * The output depends on the configuration. The maximal output contains 3
+ * different types of output lines:
+ * - "G\tgx\tgx: x and y position of a single grid generator (optional).
+ * - "C\tcx\tcy\tV\tnface": centroid position, volume and (optionally) number
+ * of faces for a single Voronoi cell.
+ * - "F\tax\tay\tbx\tby\tleft\tngb\tright\tA\tmx\tmy": edge positions
+ * (optional), left and right generator index (and ngb cell index), surface
+ * area and midpoint position for a single two-pair interface.
+ *
+ * @param v Voronoi grid.
+ * @param file File to write to.
+ */
static inline void voronoi_write_grid(const struct voronoi *restrict v,
FILE *file) {
-#ifdef VORONOI_STORE_GENERATORS
+ /* first write the cells (and generators, if those are stored) */
for (int i = 0; i < v->number_of_cells; ++i) {
struct voronoi_cell_new *this_cell = &v->cells[i];
+#ifdef VORONOI_STORE_GENERATORS
fprintf(file, "G\t%g\t%g\n", this_cell->generator[0],
this_cell->generator[1]);
- }
#endif
+ fprintf(file, "C\t%g\t%g\t%g", this_cell->centroid[0],
+ this_cell->centroid[1], this_cell->volume);
+#ifdef VORONOI_STORE_CELL_STATS
+ fprintf(file, "\t%i", this_cell->nface);
+#endif
+ fprintf(file, "\n");
+ }
+ /* now write the pairs */
for (int ngb = 0; ngb < 27; ++ngb) {
for (int i = 0; i < v->pair_index[ngb]; ++i) {
struct voronoi_pair *pair = &v->pairs[ngb][i];
+ fprintf(file, "F\t");
#ifdef VORONOI_STORE_CONNECTIONS
- fprintf(file, "F\t%g\t%g\t%g\t%g\t%i\t%i\t%i\n", pair->a[0], pair->a[1],
- pair->b[0], pair->b[1], pair->left, ngb, pair->right);
-#else
- fprintf(file, "F\t%i\t%i\t%i\n", pair->left, ngb, pair->right);
+ fprintf(file, "%g\t%g\t%g\t%g\t", pair->a[0], pair->a[1], pair->b[0],
+ pair->b[1]);
#endif
+ fprintf(file, "%i\t%i\t%i\t%g\t%g\t%g\n", pair->left, ngb, pair->right,
+ pair->surface_area, pair->midpoint[0], pair->midpoint[1]);
}
}
}
@@ -482,16 +518,6 @@ static inline void voronoi_write_grid(const struct voronoi *restrict v,
/**
* @brief Print the Voronoi grid to a file with the given name.
*
- * The grid is output as follows:
- * 1. First, each generator is output, together with all its connections.
- * The generator is output as "G\tx\ty\n", where x and y are the coordinates
- * of the generator. The centroid of the corresponding cell is output as
- * "M\tx\ty\n".
- * The connections are output as "C\tindex\tindex\n", where the two indices
- * are the indices of two vertices of the grid, in the order output by 2.
- * The midpoint of each edge is output as "F\tx\ty\n".
- * 2. Next, all vertices of the grid are output, in the format "V\tx\ty\n".
- *
* @param v Voronoi grid (read-only).
* @param file_name Name of the output file.
*/