Skip to content
GitLab
Commit b04a14f4 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'master' of gitlab.cosma.dur.ac.uk:swift/swiftsim

parents cc293259 f7ae078d
Hide whitespace changes
Inline Side-by-side
src/const.h
View file @ b04a14f4
......@@ -55,6 +55,10 @@
#define SLOPE_LIMITER_PER_FACE
#define SLOPE_LIMITER_CELL_WIDE
/* Options to control the movement of particles for GIZMO_SPH. */
/* This option disables particle movement */
//#define GIZMO_FIX_PARTICLES
/* Source terms */
#define SOURCETERMS_NONE
//#define SOURCETERMS_SN_FEEDBACK
......
src/hydro/Default/hydro.h
View file @ b04a14f4
......@@ -148,6 +148,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return min(dt_cfl, dt_u_change);
}
/**
* @brief Does some extra hydro operations once the actual physical time step
* for the particle is known.
*
* @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 Prepares a particle for the density calculation.
*
......
src/hydro/Gadget2/hydro.h
View file @ b04a14f4
......@@ -152,6 +152,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return dt_cfl;
}
/**
* @brief Does some extra hydro operations once the actual physical time step
* for the particle is known.
*
* @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 Prepares a particle for the density calculation.
*
......
src/hydro/Gizmo/hydro.h
View file @ b04a14f4
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
......@@ -40,6 +41,27 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return CFL_condition * p->h / 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) {
p->force.dt = dt;
p->force.active = 0;
}
/**
* @brief Initialises the particles for the first time
*
......@@ -58,13 +80,29 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part* p, struct xpart* xp) {
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
const float mass = p->conserved.mass;
p->primitives.v[0] = p->v[0];
p->primitives.v[1] = p->v[1];
p->primitives.v[2] = p->v[2];
/* we can already initialize the momentum */
p->conserved.momentum[0] = mass * p->primitives.v[0];
p->conserved.momentum[1] = mass * p->primitives.v[1];
p->conserved.momentum[2] = mass * p->primitives.v[2];
/* and the thermal energy */
p->conserved.energy *= mass;
#if defined(GIZMO_FIX_PARTICLES)
p->v[0] = 0.;
p->v[1] = 0.;
p->v[2] = 0.;
#else
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
#endif
}
/**
......@@ -89,6 +127,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->geometry.matrix_E[2][0] = 0.0f;
p->geometry.matrix_E[2][1] = 0.0f;
p->geometry.matrix_E[2][2] = 0.0f;
/* Set the active flag to active. */
p->force.active = 1;
}
/**
......@@ -159,6 +200,12 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->primitives.v[2] = momentum[2] / m;
const float energy = p->conserved.energy;
p->primitives.P = hydro_gamma_minus_one * energy / volume;
/* sanity checks */
if (p->primitives.rho < 0.0f || p->primitives.P < 0.0f) {
p->primitives.rho = 0.0f;
p->primitives.P = 0.0f;
}
}
/**
......@@ -180,9 +227,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* restrict p, struct xpart* restrict xp) {
/* Set the physical time step */
p->force.dt = get_timestep(p->time_bin, 0.); // MATTHIEU 0
/* Initialize time step criterion variables */
p->timestepvars.vmax = 0.0f;
......@@ -246,40 +290,14 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* @brief Converts the hydrodynamic variables from the initial condition file to
* conserved variables that can be used during the integration
*
* Requires the volume to be known.
*
* The initial condition file contains a mixture of primitive and conserved
* variables. Mass is a conserved variable, and we just copy the particle
* mass into the corresponding conserved quantity. We need the volume to
* also derive a density, which is then used to convert the internal energy
* to a pressure. However, we do not actually use these variables anymore.
* We do need to initialize the linear momentum, based on the mass and the
* velocity of the particle.
* We no longer do this, as the mass needs to be provided in the initial
* condition file, and the mass alone is enough to initialize all conserved
* variables. This is now done in hydro_first_init_part.
*
* @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part* p, struct xpart* xp) {
const float volume = p->geometry.volume;
const float m = p->conserved.mass;
p->primitives.rho = m / volume;
/* first get the initial velocities, as they were overwritten in end_density
*/
p->primitives.v[0] = p->v[0];
p->primitives.v[1] = p->v[1];
p->primitives.v[2] = p->v[2];
p->conserved.momentum[0] = m * p->primitives.v[0];
p->conserved.momentum[1] = m * p->primitives.v[1];
p->conserved.momentum[2] = m * p->primitives.v[2];
p->primitives.P =
hydro_gamma_minus_one * p->conserved.energy * p->primitives.rho;
p->conserved.energy *= m;
}
struct part* p, struct xpart* xp) {}
/**
* @brief Extra operations to be done during the drift
......@@ -362,6 +380,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
p->du_dt = 0.0f;
}
#if defined(GIZMO_FIX_PARTICLES)
p->a_hydro[0] = 0.0f;
p->a_hydro[1] = 0.0f;
p->a_hydro[2] = 0.0f;
p->du_dt = 0.0f;
/* disable the smoothing length update, since the smoothing lengths should
stay the same for all steps (particles don't move) */
p->force.h_dt = 0.0f;
#endif
}
/**
......
src/hydro/Gizmo/hydro_gradients.h
View file @ b04a14f4
......@@ -140,57 +140,61 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
/* time */
dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
Wi[2] * pi->primitives.gradients.rho[1] +
Wi[3] * pi->primitives.gradients.rho[2] +
Wi[0] * (pi->primitives.gradients.v[0][0] +
pi->primitives.gradients.v[1][1] +
pi->primitives.gradients.v[2][2]));
dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
Wi[2] * pi->primitives.gradients.v[0][1] +
Wi[3] * pi->primitives.gradients.v[0][2] +
pi->primitives.gradients.P[0] / Wi[0]);
dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
Wi[2] * pi->primitives.gradients.v[1][1] +
Wi[3] * pi->primitives.gradients.v[1][2] +
pi->primitives.gradients.P[1] / Wi[0]);
dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
Wi[2] * pi->primitives.gradients.v[2][1] +
Wi[3] * pi->primitives.gradients.v[2][2] +
pi->primitives.gradients.P[2] / Wi[0]);
dWi[4] -=
0.5 * mindt * (Wi[1] * pi->primitives.gradients.P[0] +
Wi[2] * pi->primitives.gradients.P[1] +
Wi[3] * pi->primitives.gradients.P[2] +
hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
pi->primitives.gradients.v[1][1] +
pi->primitives.gradients.v[2][2]));
dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
Wj[2] * pj->primitives.gradients.rho[1] +
Wj[3] * pj->primitives.gradients.rho[2] +
Wj[0] * (pj->primitives.gradients.v[0][0] +
pj->primitives.gradients.v[1][1] +
pj->primitives.gradients.v[2][2]));
dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
Wj[2] * pj->primitives.gradients.v[0][1] +
Wj[3] * pj->primitives.gradients.v[0][2] +
pj->primitives.gradients.P[0] / Wj[0]);
dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
Wj[2] * pj->primitives.gradients.v[1][1] +
Wj[3] * pj->primitives.gradients.v[1][2] +
pj->primitives.gradients.P[1] / Wj[0]);
dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
Wj[2] * pj->primitives.gradients.v[2][1] +
Wj[3] * pj->primitives.gradients.v[2][2] +
pj->primitives.gradients.P[2] / Wj[0]);
dWj[4] -=
0.5 * mindt * (Wj[1] * pj->primitives.gradients.P[0] +
Wj[2] * pj->primitives.gradients.P[1] +
Wj[3] * pj->primitives.gradients.P[2] +
hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
pj->primitives.gradients.v[1][1] +
pj->primitives.gradients.v[2][2]));
if (Wi[0] > 0.0f) {
dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
Wi[2] * pi->primitives.gradients.rho[1] +
Wi[3] * pi->primitives.gradients.rho[2] +
Wi[0] * (pi->primitives.gradients.v[0][0] +
pi->primitives.gradients.v[1][1] +
pi->primitives.gradients.v[2][2]));
dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
Wi[2] * pi->primitives.gradients.v[0][1] +
Wi[3] * pi->primitives.gradients.v[0][2] +
pi->primitives.gradients.P[0] / Wi[0]);
dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
Wi[2] * pi->primitives.gradients.v[1][1] +
Wi[3] * pi->primitives.gradients.v[1][2] +
pi->primitives.gradients.P[1] / Wi[0]);
dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
Wi[2] * pi->primitives.gradients.v[2][1] +
Wi[3] * pi->primitives.gradients.v[2][2] +
pi->primitives.gradients.P[2] / Wi[0]);
dWi[4] -= 0.5 * mindt *
(Wi[1] * pi->primitives.gradients.P[0] +
Wi[2] * pi->primitives.gradients.P[1] +
Wi[3] * pi->primitives.gradients.P[2] +
hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
pi->primitives.gradients.v[1][1] +
pi->primitives.gradients.v[2][2]));
}
if (Wj[0] > 0.0f) {
dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
Wj[2] * pj->primitives.gradients.rho[1] +
Wj[3] * pj->primitives.gradients.rho[2] +
Wj[0] * (pj->primitives.gradients.v[0][0] +
pj->primitives.gradients.v[1][1] +
pj->primitives.gradients.v[2][2]));
dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
Wj[2] * pj->primitives.gradients.v[0][1] +
Wj[3] * pj->primitives.gradients.v[0][2] +
pj->primitives.gradients.P[0] / Wj[0]);
dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
Wj[2] * pj->primitives.gradients.v[1][1] +
Wj[3] * pj->primitives.gradients.v[1][2] +
pj->primitives.gradients.P[1] / Wj[0]);
dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
Wj[2] * pj->primitives.gradients.v[2][1] +
Wj[3] * pj->primitives.gradients.v[2][2] +
pj->primitives.gradients.P[2] / Wj[0]);
dWj[4] -= 0.5 * mindt *
(Wj[1] * pj->primitives.gradients.P[0] +
Wj[2] * pj->primitives.gradients.P[1] +
Wj[3] * pj->primitives.gradients.P[2] +
hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
pj->primitives.gradients.v[1][1] +
pj->primitives.gradients.v[2][2]));
}
Wi[0] += dWi[0];
Wi[1] += dWi[1];
......
src/hydro/Gizmo/hydro_iact.h
View file @ b04a14f4
......@@ -231,7 +231,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
dtj = pj->force.dt;
/* calculate the maximal signal velocity */
if (Wi[0] && Wj[0]) {
if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
vmax =
sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
} else {
......@@ -412,11 +412,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
==> we update particle j if (MODE IS 1) OR (j IS INACTIVE)
*/
// MATTHIEU
const integertime_t pj_ti_end = 0; // get_integer_time_end(pj->time_bin);
const integertime_t pi_ti_end = 0; // get_integer_time_end(pi->time_bin);
if (mode == 1 || pj_ti_end > pi_ti_end) {
if (mode == 1 || pj->force.active == 0) {
/* Store mass flux */
mflux = dtj * Anorm * totflux[0];
pj->gravity.mflux[0] -= mflux * dx[0];
......
src/hydro/Gizmo/hydro_io.h
View file @ b04a14f4
......@@ -23,6 +23,13 @@
#include "io_properties.h"
#include "riemann.h"
/* Set the description of the particle movement. */
#if defined(GIZMO_FIX_PARTICLES)
#define GIZMO_PARTICLE_MOVEMENT "Fixed particles."
#else
#define GIZMO_PARTICLE_MOVEMENT "Particles move with flow velocity."
#endif
/**
* @brief Specifies which particle fields to read from a dataset
*
......@@ -155,6 +162,9 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Riemann solver information */
io_write_attribute_s(h_grpsph, "Riemann solver type",
RIEMANN_SOLVER_IMPLEMENTATION);
/* Particle movement information */
io_write_attribute_s(h_grpsph, "Particle movement", GIZMO_PARTICLE_MOVEMENT);
}
/**
......
src/hydro/Gizmo/hydro_part.h
View file @ b04a14f4
......@@ -178,6 +178,9 @@ struct part {
/* Physical time step of the particle. */
float dt;
/* Flag keeping track of whether this is an active or inactive particle. */
char active;
/* Actual velocity of the particle. */
float v_full[3];
......
src/hydro/Minimal/hydro.h
View file @ b04a14f4
......@@ -165,6 +165,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return dt_cfl;
}
/**
* @brief Does some extra hydro operations once the actual physical time step
* for the particle is known.
*
* @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 Prepares a particle for the density calculation.
*
......
src/hydro/PressureEntropy/hydro.h
View file @ b04a14f4
......@@ -152,6 +152,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return dt_cfl;
}
/**
* @brief Does some extra hydro operations once the actual physical time step
* for the particle is known.
*
* @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 Prepares a particle for the density calculation.
*
......
src/runner.c
View file @ b04a14f4
......@@ -1096,6 +1096,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
struct xpart *restrict xparts = c->xparts;
struct gpart *restrict gparts = c->gparts;
struct spart *restrict sparts = c->sparts;
const double timeBase = e->timeBase;
TIMER_TIC;
......@@ -1131,6 +1132,10 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
p->time_bin = get_time_bin(ti_new_step);
if (p->gpart != NULL) p->gpart->time_bin = get_time_bin(ti_new_step);
/* Tell the particle what the new physical time step is */
float dt = get_timestep(p->time_bin, timeBase);
hydro_timestep_extra(p, dt);
/* Number of updated particles */
updated++;
if (p->gpart != NULL) g_updated++;
......
tests/benchmarkInteractions.c
View file @ b04a14f4
......@@ -91,7 +91,10 @@ struct part *make_particles(size_t count, double *offset, double spacing,
p->h = h;
p->id = ++(*partId);
#if !defined(GIZMO_SPH)
p->mass = 1.0f;
#endif
/* Place rest of particles around the test particle
* with random position within a unit sphere. */
......@@ -110,7 +113,9 @@ struct part *make_particles(size_t count, double *offset, double spacing,
p->h = h;
p->id = ++(*partId);
#if !defined(GIZMO_SPH)
p->mass = 1.0f;
#endif
}
return particles;
}
......@@ -120,6 +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)
struct part *p;
for (size_t i = 0; i < count; ++i) {
p = &parts[i];
......@@ -130,6 +136,7 @@ void prepare_force(struct part *parts, size_t count) {
p->force.v_sig = 0.0f;
p->force.h_dt = 0.0f;
}
#endif
}
/**
......@@ -144,10 +151,19 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
"%13e %13e %13e %13e "
"%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], p->rho,
p->density.rho_dh, p->density.wcount, p->density.wcount_dh,
p->force.h_dt, p->force.v_sig,
#if defined(MINIMAL_SPH)
p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
#if defined(GIZMO_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)
0.,
#else
p->force.v_sig,
#endif
#if defined(MINIMAL_SPH) || defined(GIZMO_SPH)
0., 0., 0., 0.
#else
p->density.div_v, p->density.rot_v[0], p->density.rot_v[1],
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment