Commit a66c7479 authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Replaced global variables by a new hydro_space struct that is part of space....

Replaced global variables by a new hydro_space struct that is part of space. This struct is then passed on to hydro_init_part, and used to initialize the Voronoi cells. For non SHADOWFAX_SPH schemes, this struct is empty. Replaced random_uniform declarations in some unit tests with the tools.h version. Added some more comments to 3D Voronoi code; it is now almost completely documented (meaning I almost completely understand it now).
parent 81997283
......@@ -101,10 +101,6 @@ void print_help_message() {
"parameter file.\n");
}
#if defined(SHADOWFAX_SPH)
VORONOI_DECLARE_GLOBAL_VARIABLES()
#endif
/**
* @brief Main routine that loads a few particles and generates some output.
*
......@@ -390,13 +386,6 @@ int main(int argc, char *argv[]) {
struct gravity_props gravity_properties;
if (with_self_gravity) gravity_props_init(&gravity_properties, params);
#if defined(SHADOWFAX_SPH)
/* Override the variables governing the density iteration
(see src/hydro/Shadowswift/hydro.h for full explanation) */
hydro_properties.target_neighbours = 1.0f;
hydro_properties.delta_neighbours = 0.0f;
#endif
/* Read particles and space information from (GADGET) ICs */
char ICfileName[200] = "";
parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
......@@ -439,27 +428,6 @@ int main(int argc, char *argv[]) {
fflush(stdout);
}
#if defined(SHADOWFAX_SPH)
/* set the *global* box dimensions */
float box_anchor[3], box_side[3];
if (periodic) {
box_anchor[0] = -0.5f * dim[0];
box_anchor[1] = -0.5f * dim[1];
box_anchor[2] = -0.5f * dim[2];
box_side[0] = 2.0f * dim[0];
box_side[1] = 2.0f * dim[1];
box_side[2] = 2.0f * dim[2];
} else {
box_anchor[0] = 0.0f;
box_anchor[1] = 0.0f;
box_anchor[2] = 0.0f;
box_side[0] = dim[0];
box_side[1] = dim[1];
box_side[2] = dim[2];
}
voronoi_set_box(box_anchor, box_side);
#endif
#ifdef SWIFT_DEBUG_CHECKS
/* Check once and for all that we don't have unwanted links */
if (!with_stars) {
......
......@@ -45,7 +45,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
physical_constants.h physical_constants_cgs.h potential.h version.h \
hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
hydro_space.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
......@@ -55,7 +56,8 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
physical_constants.c potential.c hydro_properties.c \
runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \
statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
part_type.c xmf.c gravity_properties.c gravity.c
part_type.c xmf.c gravity_properties.c gravity.c \
hydro_space.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
......
......@@ -22,6 +22,7 @@
#include "adiabatic_index.h"
#include "approx_math.h"
#include "equation_of_state.h"
#include "hydro_space.h"
#include "minmax.h"
#include <float.h>
......@@ -165,9 +166,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
* the variaous density tasks
*
* @param p The particle to act upon
* @param hs #hydro_space containing hydro specific space information.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part *restrict p) {
struct part *restrict p, const struct hydro_space *hs) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
......@@ -400,7 +402,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->u_full = p->u;
hydro_reset_acceleration(p);
hydro_init_part(p);
hydro_init_part(p, NULL);
}
#endif /* SWIFT_DEFAULT_HYDRO_H */
......@@ -36,6 +36,7 @@
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "hydro_space.h"
#include "kernel_hydro.h"
#include "minmax.h"
......@@ -169,9 +170,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
* the variaous density tasks
*
* @param p The particle to act upon
* @param hs #hydro_space containing hydro specific space information.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part *restrict p) {
struct part *restrict p, const struct hydro_space *hs) {
p->rho = 0.f;
p->density.wcount = 0.f;
......@@ -456,7 +458,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->entropy_full = p->entropy;
hydro_reset_acceleration(p);
hydro_init_part(p);
hydro_init_part(p, NULL);
}
#endif /* SWIFT_GADGET2_HYDRO_H */
......@@ -23,6 +23,7 @@
#include "approx_math.h"
#include "equation_of_state.h"
#include "hydro_gradients.h"
#include "hydro_space.h"
#include "minmax.h"
#include "riemann.h"
......@@ -145,9 +146,10 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
* Simply makes sure all necessary variables are initialized to zero.
*
* @param p The particle to act upon
* @param hs #hydro_space containing hydro specific space information.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part* p) {
struct part* p, const struct hydro_space* hs) {
p->density.wcount = 0.0f;
p->density.wcount_dh = 0.0f;
......
......@@ -38,6 +38,7 @@
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "hydro_space.h"
#include "kernel_hydro.h"
#include "minmax.h"
......@@ -183,9 +184,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
* density sub-structure of a particle get zeroed in here.
*
* @param p The particle to act upon
* @param hs #hydro_space containing hydro specific space information.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part *restrict p) {
struct part *restrict p, const struct hydro_space *hs) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
......@@ -429,7 +431,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->u_full = p->u;
hydro_reset_acceleration(p);
hydro_init_part(p);
hydro_init_part(p, NULL);
}
#endif /* SWIFT_MINIMAL_HYDRO_H */
......@@ -36,6 +36,7 @@
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "hydro_space.h"
#include "kernel_hydro.h"
#include "minmax.h"
......@@ -169,9 +170,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
* the variaous density tasks
*
* @param p The particle to act upon
* @param hs #hydro_space containing hydro specific space information.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part *restrict p) {
struct part *restrict p, const struct hydro_space *hs) {
p->rho = 0.f;
p->rho_bar = 0.f;
......@@ -474,7 +476,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->v_full[2] = p->v[2];
hydro_reset_acceleration(p);
hydro_init_part(p);
hydro_init_part(p, NULL);
}
#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_H */
......@@ -22,6 +22,7 @@
#include "approx_math.h"
#include "equation_of_state.h"
#include "hydro_gradients.h"
#include "hydro_space.h"
#include "voronoi_algorithm.h"
/**
......@@ -128,14 +129,15 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
* Initializes the Voronoi cell.
*
* @param p The particle to act upon
* @param hs #hydro_space containing extra information about the space.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part* p) {
struct part* p, const struct hydro_space* hs) {
p->density.wcount = 0.0f;
p->density.wcount_dh = 0.0f;
voronoi_cell_init(&p->cell, p->x);
voronoi_cell_init(&p->cell, p->x, hs->anchor, hs->side);
/* Set the active flag to active. */
p->force.active = 1;
......
......@@ -20,15 +20,12 @@
#ifndef SWIFT_VORONOIXD_ALGORITHM_H
#define SWIFT_VORONOIXD_ALGORITHM_H
#include <float.h>
#include <math.h>
#include <stdlib.h>
#include "error.h"
#include "inline.h"
#include "voronoi1d_cell.h"
#define VORONOI_DECLARE_GLOBAL_VARIABLES()
/**
* @brief Store the extents of the simulation box in the global variables.
*
......@@ -47,12 +44,15 @@ __attribute__((always_inline)) INLINE static void voronoi_set_box(
*
* @param cell 1D Voronoi cell to initialize.
* @param x Position of the generator of the cell.
* @param anchor Anchor of the simulation box.
* @param side Side lengths of the simulation box.
*/
__attribute__((always_inline)) INLINE void voronoi_cell_init(
struct voronoi_cell *cell, const double *x) {
struct voronoi_cell *cell, const double *x, const double *anchor,
const double *side) {
cell->x = x[0];
cell->xL = -DBL_MAX;
cell->xR = DBL_MAX;
cell->xL = anchor[0] - cell->x;
cell->xR = anchor[0] + side[0] - cell->x;
cell->idL = 0;
cell->idR = 0;
cell->volume = 0.0f;
......
......@@ -43,48 +43,17 @@
#define VORONOI2D_BOX_TOP 18446744073709551604llu
#define VORONOI2D_BOX_BOTTOM 18446744073709551605llu
/* Global variables used to store the extent of the simulation box */
extern float global_voronoi_box_anchor[2];
extern float global_voronoi_box_side[2];
/* Macro used to declare the global variables in any piece of code that uses
this header file */
#define VORONOI_DECLARE_GLOBAL_VARIABLES() \
float global_voronoi_box_anchor[2]; \
float global_voronoi_box_side[2];
/* Macros that make it easier to change the box extents that are used throughout
the code */
#define VORONOI2D_BOX_ANCHOR_X global_voronoi_box_anchor[0]
#define VORONOI2D_BOX_ANCHOR_Y global_voronoi_box_anchor[1]
#define VORONOI2D_BOX_SIDE_X global_voronoi_box_side[0]
#define VORONOI2D_BOX_SIDE_Y global_voronoi_box_side[1]
/**
* @brief Set the values of the global variables that store the simulation box
* extents.
*
* @param anchor Corner of the box with the lowest coordinate values.
* @param side Side lengths of the box.
*/
__attribute__((always_inline)) INLINE static void voronoi_set_box(
const float *anchor, const float *side) {
global_voronoi_box_anchor[0] = anchor[0];
global_voronoi_box_anchor[1] = anchor[1];
global_voronoi_box_side[0] = side[0];
global_voronoi_box_side[1] = side[1];
}
/**
* @brief Initialize a 2D Voronoi cell.
*
* @param cell 2D Voronoi cell to initialize.
* @param x Position of the generator of the cell.
* @param anchor Anchor of the simulation box containing all particles.
* @param side Side lengths of the simulation box containing all particles.
*/
__attribute__((always_inline)) INLINE void voronoi_cell_init(
struct voronoi_cell *cell, const double *x) {
struct voronoi_cell *cell, const double *x, const double *anchor,
const double *side) {
/* Set the position of the generator of the cell (for reference) */
cell->x[0] = x[0];
......@@ -94,21 +63,17 @@ __attribute__((always_inline)) INLINE void voronoi_cell_init(
(note: all vertex coordinates are relative w.r.t. the cell generator) */
cell->nvert = 4;
cell->vertices[0][0] = VORONOI2D_BOX_ANCHOR_X - cell->x[0];
cell->vertices[0][1] = VORONOI2D_BOX_ANCHOR_Y - cell->x[1];
cell->vertices[0][0] = anchor[0] - cell->x[0];
cell->vertices[0][1] = anchor[1] - cell->x[1];
cell->vertices[1][0] = VORONOI2D_BOX_ANCHOR_X - cell->x[0];
cell->vertices[1][1] =
VORONOI2D_BOX_ANCHOR_Y + VORONOI2D_BOX_SIDE_Y - cell->x[1];
cell->vertices[1][0] = anchor[0] - cell->x[0];
cell->vertices[1][1] = anchor[1] + side[1] - cell->x[1];
cell->vertices[2][0] =
VORONOI2D_BOX_ANCHOR_X + VORONOI2D_BOX_SIDE_X - cell->x[0];
cell->vertices[2][1] =
VORONOI2D_BOX_ANCHOR_Y + VORONOI2D_BOX_SIDE_Y - cell->x[1];
cell->vertices[2][0] = anchor[0] + side[0] - cell->x[0];
cell->vertices[2][1] = anchor[1] + side[1] - cell->x[1];
cell->vertices[3][0] =
VORONOI2D_BOX_ANCHOR_X + VORONOI2D_BOX_SIDE_X - cell->x[0];
cell->vertices[3][1] = VORONOI2D_BOX_ANCHOR_Y - cell->x[1];
cell->vertices[3][0] = anchor[0] + side[0] - cell->x[0];
cell->vertices[3][1] = anchor[1] - cell->x[1];
/* The neighbours are ordered such that neighbour i shares the face in between
vertices i and i+1 (with last vertex + 1 = first vertex)
......
This diff is collapsed.
......@@ -48,23 +48,29 @@ struct voronoi_cell {
float vertices[3 * VORONOI3D_MAXNUMVERT];
/* Number of edges for every vertex. */
int orders[VORONOI3D_MAXNUMVERT];
char orders[VORONOI3D_MAXNUMVERT];
/* Offsets of the edges, edgeindices and neighbours corresponding to a
particular vertex in the internal arrays */
int offsets[VORONOI3D_MAXNUMVERT];
/* Edge information. */
/* Edge information. Edges are ordered counterclockwise w.r.t. a vector
pointing from the cell generator to the vertex. */
int edges[VORONOI3D_MAXNUMEDGE];
/* Additional edge information. */
int edgeindices[VORONOI3D_MAXNUMEDGE];
/* Neighbour information. */
char edgeindices[VORONOI3D_MAXNUMEDGE];
/* Neighbour information. This field is used differently depending on where we
are in the algorithm. During cell construction, it contains, for every edge
of every vertex, the index of the neighbour that generates the face
counterclockwise of the edge w.r.t. a vector pointing from the vertex along
the edge. After cell finalization, it contains a neighbour for every face,
in the same order as the face_areas and face_midpoints arrays. */
unsigned long long ngbs[VORONOI3D_MAXNUMEDGE];
/* Number of faces of the cell. */
int nface;
unsigned char nface;
/* Surface areas of the cell faces. */
float face_areas[VORONOI3D_MAXFACE];
......
......@@ -44,6 +44,12 @@ void hydro_props_init(struct hydro_props *p,
p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
p->delta_neighbours = parser_get_param_float(params, "SPH:delta_neighbours");
#ifdef SHADOWFAX_SPH
/* change the meaning of target_neighbours and delta_neighbours */
p->target_neighbours = 1.0f;
p->delta_neighbours = 0.0f;
#endif
/* Maximal smoothing length */
p->h_max = parser_get_opt_param_float(params, "SPH:h_max",
hydro_props_default_h_max);
......
......@@ -489,6 +489,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
const int count = c->count;
const int gcount = c->gcount;
const struct engine *e = r->e;
const struct space *s = e->s;
TIMER_TIC;
......@@ -514,7 +515,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
if (part_is_active(p, e)) {
/* Get ready for a density calculation */
hydro_init_part(p);
hydro_init_part(p, &s->hs);
}
}
......@@ -596,6 +597,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
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 =
......@@ -676,7 +678,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
redo += 1;
/* Re-initialise everything */
hydro_init_part(p);
hydro_init_part(p, &s->hs);
/* Off we go ! */
continue;
......
......@@ -2692,6 +2692,8 @@ void space_init(struct space *s, const struct swift_params *params,
bzero(s->xparts, Npart * sizeof(struct xpart));
}
hydro_space_init(&s->hs, s);
/* Set the particles in a state where they are ready for a run */
space_init_parts(s);
space_init_xparts(s);
......
......@@ -31,6 +31,7 @@
/* Includes. */
#include "cell.h"
#include "hydro_space.h"
#include "lock.h"
#include "parser.h"
#include "part.h"
......@@ -69,6 +70,9 @@ struct space {
/*! Is the space periodic? */
int periodic;
/*! Extra space information needed for some hydro schemes. */
struct hydro_space hs;
/*! Are we doing gravity? */
int gravity;
......
......@@ -144,7 +144,7 @@ void pairs_single_density(double *dim, long long int pid,
p = parts[k];
printf("pairs_single: part[%i].id == %lli.\n", k, pid);
hydro_init_part(&p);
hydro_init_part(&p, NULL);
/* Loop over all particle pairs. */
for (k = 0; k < N; k++) {
......@@ -459,7 +459,7 @@ void engine_single_density(double *dim, long long int pid,
p = parts[k];
/* Clear accumulators. */
hydro_init_part(&p);
hydro_init_part(&p, NULL);
/* Loop over all particle pairs (force). */
for (k = 0; k < N; k++) {
......
......@@ -450,10 +450,6 @@ void test_interactions(struct part test_part, struct part *parts, size_t count,
#endif
}
#if defined(SHADOWFAX_SPH)
VORONOI_DECLARE_GLOBAL_VARIABLES()
#endif
/* And go... */
int main(int argc, char *argv[]) {
size_t runs = 10000;
......
......@@ -422,10 +422,6 @@ void runner_doself1_density(struct runner *r, struct cell *ci);
void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself2_force(struct runner *r, struct cell *ci);
#if defined(SHADOWFAX_SPH)
VORONOI_DECLARE_GLOBAL_VARIABLES()
#endif
/* And go... */
int main(int argc, char *argv[]) {
......
......@@ -119,7 +119,13 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
part->conserved.mass = density * volume / count;
voronoi_cell_init(&part->cell, part->x);
#ifdef SHADOWFAX_SPH
double anchor[3] = {0., 0., 0.};
double side[3] = {1., 1., 1.};
voronoi_cell_init(&part->cell, part->x, anchor, side);
#endif
#else
part->mass = density * volume / count;
#endif
......@@ -177,7 +183,7 @@ void clean_up(struct cell *ci) {
*/
void zero_particle_fields(struct cell *c) {
for (int pid = 0; pid < c->count; pid++) {
hydro_init_part(&c->parts[pid]);
hydro_init_part(&c->parts[pid], NULL);
}
}
......@@ -299,10 +305,6 @@ void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
void runner_doself1_density(struct runner *r, struct cell *ci);
void runner_doself1_density_vec(struct runner *r, struct cell *ci);
#if defined(SHADOWFAX_SPH)
VORONOI_DECLARE_GLOBAL_VARIABLES()
#endif
/* And go... */
int main(int argc, char *argv[]) {
......
......@@ -22,10 +22,7 @@
#include "const.h"
#include "dimension.h"
#include "error.h"
double random_uniform(double a, double b) {
return (rand() / (double)RAND_MAX) * (b - a) + a;
}
#include "tools.h"
void setup_matrix(float A[3][3]) {
A[0][0] = random_uniform(-1.0, 1.0);
......
......@@ -116,7 +116,7 @@ void clean_up(struct cell *ci) {
*/
void zero_particle_fields(struct cell *c) {
for (int pid = 0; pid < c->count; pid++) {
hydro_init_part(&c->parts[pid]);
hydro_init_part(&c->parts[pid], NULL);
}
}
......@@ -187,10 +187,6 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
/* Just a forward declaration... */
void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
#if defined(SHADOWFAX_SPH)
VORONOI_DECLARE_GLOBAL_VARIABLES()
#endif
int main(int argc, char *argv[]) {
size_t particles = 0, runs = 0, volume, type = 0;
double offset[3] = {0, 0, 0}, h = 1.1255, size = 1., rho = 1.;
......
......@@ -20,10 +20,7 @@
#include <string.h>
#include "error.h"
#include "riemann/riemann_exact.h"
double random_uniform(double a, double b) {
return (rand() / (double)RAND_MAX) * (b - a) + a;
}
#include "tools.h"
int opposite(float a, float b) {
if ((a - b)) {
......
......@@ -20,10 +20,7 @@
#include <string.h>
#include "error.h"
#include "riemann/riemann_hllc.h"
double random_uniform(double a, double b) {
return (rand() / (double)RAND_MAX) * (b - a) + a;
}
#include "tools.h"
int consistent_with_zero(float val) { return fabs(val) < 1.e-4; }
......
......@@ -20,10 +20,7 @@
#include <string.h>
#include "error.h"
#include "riemann/riemann_trrs.h"
double random_uniform(double a, double b) {
return (rand() / (double)RAND_MAX) * (b - a) + a;
}
#include "tools.h"