Skip to content
Snippets Groups Projects
Commit 3fc7d30c authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

1D SHADOWFAX_SPH now compiles again.

parent 1c8b3ec6
Branches
Tags
1 merge request!3211D and 2D moving mesh algorithm
......@@ -49,7 +49,7 @@
#include "./hydro/Default/hydro_debug.h"
#elif defined(GIZMO_SPH)
#include "./hydro/Gizmo/hydro_debug.h"
#elif defined(SHADOWSWIFT)
#elif defined(SHADOWFAX_SPH)
#include "./hydro/Shadowswift/hydro_debug.h"
#else
#error "Invalid choice of SPH variant"
......
......@@ -42,6 +42,23 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return CFL_condition * R / fabsf(p->timestepvars.vmax);
}
/**
* @brief Does some extra hydro operations once the actual physical time step
* for the particle is known.
*
* We use this to store the physical time step, since it is used for the flux
* exchange during the force loop.
*
* We also set the active flag of the particle to inactive. It will be set to
* active in hydro_init_part, which is called the next time the particle becomes
* active.
*
* @param p The particle to act upon.
* @param dt Physical time step of the particle during the next step.
*/
__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
struct part* p, float dt) {}
/**
* @brief Initialises the particles for the first time
*
......@@ -99,10 +116,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
* This method also initializes the gradient variables (if gradients are used).
*
* @param p The particle to act upon.
* @param The current physical time.
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
struct part* restrict p, float time) {
struct part* restrict p) {
float volume;
float m, momentum[3], energy;
......@@ -164,11 +180,10 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
* @param timeBase Conversion factor between integer time and physical time.
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* restrict p, struct xpart* restrict xp, int ti_current,
double timeBase) {
struct part* restrict p, struct xpart* restrict xp) {
/* Set the physical time step */
p->force.dt = (p->ti_end - p->ti_begin) * timeBase;
// p->force.dt = (p->ti_end - p->ti_begin) * timeBase;
/* Initialize time step criterion variables */
p->timestepvars.vmax = 0.0f;
......@@ -213,6 +228,16 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
p->force.h_dt = 0.0f;
}
/**
* @brief Sets the values to be predicted in the drifts to their values at a
* kick time
*
* @param p The particle.
* @param xp The extended data of this particle.
*/
__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
struct part* restrict p, const struct xpart* restrict xp) {}
/**
* @brief Converts the hydrodynamic variables from the initial condition file to
* conserved variables that can be used during the integration
......@@ -228,9 +253,10 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
* velocity of the particle.
*
* @param p The particle to act upon.
* @param xp The extended particle data to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part* p) {
struct part* p, struct xpart* xp) {
float volume;
float m;
......@@ -258,13 +284,9 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
* @param p Particle to act upon.
* @param xp The extended particle data to act upon.
* @param dt The drift time-step.
* @param t0 Integer start time of the drift interval.
* @param t1 Integer end time of the drift interval.
* @param timeBase Conversion factor between integer and physical time.
*/
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part* p, struct xpart* xp, float dt, int t0, int t1,
double timeBase) {}
struct part* p, struct xpart* xp, float dt) {}
/**
* @brief Set the particle acceleration after the flux loop
......@@ -361,20 +383,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param p Particle to act upon.
* @param xp Extended particle data to act upon.
* @param dt Physical time step.
* @param half_dt Half the physical time step.
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt, float half_dt) {}
struct part* p, struct xpart* xp, float dt) {}
/**
* @brief Returns the internal energy of a particle
*
* @param p The particle of interest.
* @param dt Time since the last kick.
* @return Internal energy of the particle.
*/
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const struct part* restrict p, float dt) {
const struct part* restrict p) {
return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
}
......@@ -383,11 +403,10 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
* @brief Returns the entropy of a particle
*
* @param p The particle of interest.
* @param dt Time since the last kick.
* @return Entropy of the particle.
*/
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
const struct part* restrict p, float dt) {
const struct part* restrict p) {
return p->primitives.P / pow_gamma(p->primitives.rho);
}
......@@ -396,11 +415,10 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
* @brief Returns the sound speed of a particle
*
* @param p The particle of interest.
* @param dt Time since the last kick.
* @param Sound speed of the particle.
*/
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part* restrict p, float dt) {
const struct part* restrict p) {
return sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
}
......@@ -409,11 +427,10 @@ __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
* @brief Returns the pressure of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
* @param Pressure of the particle.
*/
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
const struct part* restrict p, float dt) {
const struct part* restrict p) {
return p->primitives.P;
}
......
......@@ -143,18 +143,18 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
*/
void writeSPHflavour(hid_t h_grpsph) {
/* Gradient information */
writeAttribute_s(h_grpsph, "Gradient reconstruction model",
HYDRO_GRADIENT_IMPLEMENTATION);
io_write_attribute_s(h_grpsph, "Gradient reconstruction model",
HYDRO_GRADIENT_IMPLEMENTATION);
/* Slope limiter information */
writeAttribute_s(h_grpsph, "Cell wide slope limiter model",
HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION);
writeAttribute_s(h_grpsph, "Piecewise slope limiter model",
HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION);
io_write_attribute_s(h_grpsph, "Cell wide slope limiter model",
HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION);
io_write_attribute_s(h_grpsph, "Piecewise slope limiter model",
HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION);
/* Riemann solver information */
writeAttribute_s(h_grpsph, "Riemann solver type",
RIEMANN_SOLVER_IMPLEMENTATION);
io_write_attribute_s(h_grpsph, "Riemann solver type",
RIEMANN_SOLVER_IMPLEMENTATION);
}
/**
......
......@@ -32,7 +32,7 @@
#include "./hydro/Default/hydro_io.h"
#elif defined(GIZMO_SPH)
#include "./hydro/Gizmo/hydro_io.h"
#elif defined(SHADOWSWIFT)
#elif defined(SHADOWFAX_SPH)
#include "./hydro/Shadowswift/hydro_io.h"
#else
#error "Invalid choice of SPH variant"
......
......@@ -92,7 +92,7 @@ struct part *make_particles(size_t count, double *offset, double spacing,
p->h = h;
p->id = ++(*partId);
#if !defined(GIZMO_SPH)
#if !defined(GIZMO_SPH) && !defined(SHADOWFAX_SPH)
p->mass = 1.0f;
#endif
......@@ -113,7 +113,7 @@ struct part *make_particles(size_t count, double *offset, double spacing,
p->h = h;
p->id = ++(*partId);
#if !defined(GIZMO_SPH)
#if !defined(GIZMO_SPH) && !defined(SHADOWFAX_SPH)
p->mass = 1.0f;
#endif
}
......@@ -125,7 +125,7 @@ struct part *make_particles(size_t count, double *offset, double spacing,
*/
void prepare_force(struct part *parts, size_t count) {
#if !defined(GIZMO_SPH)
#if !defined(GIZMO_SPH) && !defined(SHADOWFAX_SPH)
struct part *p;
for (size_t i = 0; i < count; ++i) {
p = &parts[i];
......@@ -152,18 +152,18 @@ 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->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
#if defined(GIZMO_SPH)
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
0., 0.,
#else
p->rho, p->density.rho_dh,
#endif
p->density.wcount, p->density.wcount_dh, p->force.h_dt,
#if defined(GIZMO_SPH)
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
0.,
#else
p->force.v_sig,
#endif
#if defined(MINIMAL_SPH) || defined(GIZMO_SPH)
#if defined(MINIMAL_SPH) || defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
0., 0., 0., 0.
#else
p->density.div_v, p->density.rot_v[0], p->density.rot_v[1],
......
......@@ -101,7 +101,7 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(MINIMAL_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(GIZMO_SPH) || defined(SHADOWSWIFT)
#elif defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
part->primitives.P = pressure;
#else
error("Need to define pressure here !");
......@@ -198,7 +198,7 @@ void reset_particles(struct cell *c, enum velocity_field vel,
set_velocity(p, vel, size);
set_energy_state(p, press, size, density);
#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
float volume = p->conserved.mass / density;
#if defined(GIZMO_SPH)
p->geometry.volume = volume;
......@@ -265,7 +265,7 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
part->x[2] = offset[2] + size * (z + 0.5) / (float)n;
part->h = size * h / (float)n;
#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
part->conserved.mass = density * volume / count;
#else
part->mass = density * volume / count;
......@@ -368,7 +368,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,
hydro_get_density(&main_cell->parts[pid]),
#if defined(MINIMAL_SPH) || defined(SHADOWSWIFT)
#if defined(MINIMAL_SPH) || defined(SHADOWFAX_SPH)
0.f,
#else
main_cell->parts[pid].density.div_v,
......@@ -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(SHADOWSWIFT) && defined(HYDRO_DIMENSION_3D)
VORONOI3D_DECLARE_GLOBAL_VARIABLES()
#endif
/* And go... */
int main(int argc, char *argv[]) {
......@@ -446,12 +442,6 @@ int main(int argc, char *argv[]) {
/* Get some randomness going */
srand(0);
#if defined(SHADOWSWIFT) && defined(HYDRO_DIMENSION_3D)
float box_anchor[3] = {-2.0f, -2.0f, -2.0f};
float box_side[3] = {8.0f, 8.0f, 8.0f};
voronoi_set_box(box_anchor, box_side);
#endif
char c;
while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:p:")) != -1) {
switch (c) {
......
......@@ -216,7 +216,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],
hydro_get_density(&main_cell->parts[pid]),
#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
0.f,
#else
main_cell->parts[pid].density.rho_dh,
......@@ -253,7 +253,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
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], hydro_get_density(&cj->parts[pjd]),
#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
0.f,
#else
main_cell->parts[pjd].density.rho_dh,
......@@ -299,10 +299,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(SHADOWSWIFT) && defined(HYDRO_DIMENSION_3D)
VORONOI3D_DECLARE_GLOBAL_VARIABLES()
#endif
/* And go... */
int main(int argc, char *argv[]) {
......@@ -325,12 +321,6 @@ int main(int argc, char *argv[]) {
/* Get some randomness going */
srand(0);
#if defined(SHADOWSWIFT) && defined(HYDRO_DIMENSION_3D)
float box_anchor[3] = {-2.0f, -2.0f, -2.0f};
float box_side[3] = {6.0f, 6.0f, 6.0f};
voronoi_set_box(box_anchor, box_side);
#endif
char c;
while ((c = getopt(argc, argv, "m:s:h:n:r:t:d:f:v:a:")) != -1) {
switch (c) {
......
......@@ -142,7 +142,7 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
ci->parts[pid].id, ci->parts[pid].x[0], ci->parts[pid].x[1],
ci->parts[pid].x[2], ci->parts[pid].v[0], ci->parts[pid].v[1],
ci->parts[pid].v[2], hydro_get_density(&ci->parts[pid]),
#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
0.f,
#else
ci->parts[pid].density.rho_dh,
......@@ -166,7 +166,7 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
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], hydro_get_density(&cj->parts[pjd]),
#if defined(GIZMO_SPH) || defined(SHADOWSWIFT)
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
0.f,
#else
cj->parts[pjd].density.rho_dh,
......@@ -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(SHADOWSWIFT) && defined(HYDRO_DIMENSION_3D)
VORONOI3D_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.;
......@@ -211,12 +207,6 @@ int main(int argc, char *argv[]) {
srand(0);
#if defined(SHADOWSWIFT) && defined(HYDRO_DIMENSION_3D)
float box_anchor[3] = {-2.0f, -2.0f, -2.0f};
float box_side[3] = {6.0f, 6.0f, 6.0f};
voronoi_set_box(box_anchor, box_side);
#endif
while ((c = getopt(argc, argv, "h:p:r:t:d:f:")) != -1) {
switch (c) {
case 'h':
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment