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

Rewrote the Gadget-2 hydro implementation. Everything OK until the end of the density loop.

parent 626f555a
Expand all Hide whitespace changes
Inline Side-by-side
src/cell.c
View file @ a3596a2e
......@@ -564,12 +564,27 @@ void cell_init_parts(struct cell *c, void *data) {
xp[i].v_full[0] = p[i].v[0];
xp[i].v_full[1] = p[i].v[1];
xp[i].v_full[2] = p[i].v[2];
hydro_init_part(p);
hydro_reset_acceleration(p);
hydro_init_part(&p[i]);
hydro_reset_acceleration(&p[i]);
}
c->t_end_min = 0.;
}
/**
* @brief Converts hydro quantities to a valid state after the initial density calculation
*
* @param c Cell to act upon
* @param data Unused parameter
*/
void cell_convert_hydro(struct cell *c, void *data) {
struct part *p = c->parts;
for(int i=0; i<c->count; ++i) {
hydro_convert_quantities(&p[i]);
}
}
/**
* @brief Cleans the links in a given cell.
......
src/cell.h
View file @ a3596a2e
......@@ -178,6 +178,7 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s);
int cell_getsize(struct cell *c);
int cell_link(struct cell *c, struct part *parts);
void cell_init_parts(struct cell *c, void *data);
void cell_convert_hydro(struct cell *c, void *data);
void cell_clean_links(struct cell * c, void * data);
#endif /* SWIFT_CELL_H */
src/debug.c
View file @ a3596a2e
......@@ -46,17 +46,24 @@ void printParticle(struct part *parts, long long int id, int N) {
for (i = 0; i < N; i++)
if (parts[i].id == id) {
printf(
"## Particle[%d]: id=%lld, x=[%.16e,%.16e,%.16e], "
"v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], h=%.3e, h_dt=%.3e, "
"wcount=%.3e, m=%.3e, rho=%.3e, rho_dh=%.3e, div_v=%.3e, u=%.3e, "
"dudt=%.3e, bals=%.3e, POrho2=%.3e, v_sig=%.3e, t_begin=%.3e, "
"t_end=%.3e\n",
"## Particle[%d]:\n id=%lld, x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e],\n h=%.3e, "
"wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%3.e, S=%.3e,\n divV=%.3e, "
" curlV=%.3e rotV=[%.3e,%.3e,%.3e] \n "
"t_begin=%.3e, t_end=%.3e\n",
i, parts[i].id, parts[i].x[0], parts[i].x[1], parts[i].x[2],
parts[i].v[0], parts[i].v[1], parts[i].v[2], parts[i].a[0],
parts[i].a[1], parts[i].a[2], parts[i].h, parts[i].force.h_dt,
parts[i].density.wcount, parts[i].mass, parts[i].rho, parts[i].rho_dh,
parts[i].density.div_v, parts[i].u, parts[i].force.u_dt,
parts[i].force.balsara, parts[i].force.POrho2, parts[i].force.v_sig,
parts[i].a[1], parts[i].a[2], 2.*parts[i].h,
(int)parts[i].density.wcount, parts[i].mass,
parts[i].rho_dh,
parts[i].rho,
parts[i].pressure,
parts[i].entropy,
parts[i].div_v,
parts[i].curl_v,
parts[i].rot_v[0],
parts[i].rot_v[1],
parts[i].rot_v[2],
parts[i].t_begin, parts[i].t_end);
found = 1;
}
......@@ -94,15 +101,17 @@ void printgParticle(struct gpart *parts, long long int id, int N) {
void printParticle_single(struct part *p) {
printf(
"## Particle: id=%lld, x=[%e,%e,%e], v=[%.3e,%.3e,%.3e], "
"a=[%.3e,%.3e,%.3e], h=%.3e, h_dt=%.3e, wcount=%.3e, m=%.3e, rho=%.3e, "
"rho_dh=%.3e, div_v=%.3e, u=%.3e, dudt=%.3e, bals=%.3e, POrho2=%.3e, "
"v_sig=%.3e, t_begin=%.3e, t_end=%.3e\n",
p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a[0],
p->a[1], p->a[2], p->h, p->force.h_dt, p->density.wcount, p->mass, p->rho,
p->rho_dh, p->density.div_v, p->u, p->force.u_dt, p->force.balsara,
p->force.POrho2, p->force.v_sig, p->t_begin, p->t_end);
/* printf( */
/* "## Particle: id=%lld, x=[%e,%e,%e], v=[%.3e,%.3e,%.3e], " */
/* "a=[%.3e,%.3e,%.3e], h=%.3e, h_dt=%.3e, wcount=%.3e, m=%.3e, rho=%.3e, " */
/* "rho_dh=%.3e, div_v=%.3e, u=%.3e, dudt=%.3e, bals=%.3e, POrho2=%.3e, " */
/* "v_sig=%.3e, t_begin=%.3e, t_end=%.3e\n", */
/* p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a[0], */
/* p->a[1], p->a[2], p->h, p->force.h_dt, p->density.wcount, p->mass, p->rho, */
/* p->rho_dh, p->density.div_v, p->u, p->force.u_dt, p->force.balsara, */
/* p->force.POrho2, p->force.v_sig, p->t_begin, p->t_end); */
// MATTHIEU
}
#ifdef HAVE_METIS
......
src/engine.c
View file @ a3596a2e
......@@ -1714,8 +1714,8 @@ void engine_init_particles(struct engine *e) {
message("Initialising particles");
space_map_cells_pre(s, 1, cell_init_parts, NULL);
printParticle(e->s->parts, 1000, e->s->nr_parts);
printParticle(e->s->parts, 515050, e->s->nr_parts);
//printParticle(e->s->parts, 1000, e->s->nr_parts);
//printParticle(e->s->parts, 515050, e->s->nr_parts);
message("\n0th DENSITY CALC\n");
......@@ -1729,11 +1729,15 @@ void engine_init_particles(struct engine *e) {
1 << task_subtype_density);
TIMER_TOC(timer_runners);
message("\n0th ENTROPY CONVERSION\n");
space_map_cells_pre(s, 1, cell_convert_hydro, NULL);
printParticle(e->s->parts, 1000, e->s->nr_parts);
printParticle(e->s->parts, 515050, e->s->nr_parts);
abort();
exit(0);
/* Ready to go */
e->step = -1;
......@@ -1788,13 +1792,13 @@ void engine_step(struct engine *e) {
MPI_SUCCESS)
error("Failed to aggregate t_end_max.");
t_end_max = in[0];
out[0] = count;
out[0] = updates;
out[1] = ekin;
out[2] = epot;
if (MPI_Allreduce(out, in, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate energies.");
count = in[0];
updates = in[0];
ekin = in[1];
epot = in[2];
/* int nr_parts;
......@@ -1819,8 +1823,6 @@ if ( e->nodeID == 0 )
printParticle(e->s->parts, 1000, e->s->nr_parts);
printParticle(e->s->parts, 515050, e->s->nr_parts);
abort();
/* Move forward in time */
e->timeOld = e->time;
e->time = t_end_min;
......@@ -1830,6 +1832,7 @@ if ( e->nodeID == 0 )
printf("%d %f %f %d\n", e->step, e->time, e->timeStep, updates);
fflush(stdout);
message("\nACCELERATION AND KICK\n");
/* Re-distribute the particles amongst the nodes? */
if (e->forcerepart) engine_repartition(e);
......@@ -1845,11 +1848,17 @@ if ( e->nodeID == 0 )
(1 << task_type_init) | (1 << task_type_ghost) |
(1 << task_type_kick) | (1 << task_type_send) |
(1 << task_type_recv),
0);
(1 << task_subtype_density) | (1 << task_subtype_force));
TIMER_TOC(timer_runners);
TIMER_TOC2(timer_step);
printParticle(e->s->parts, 1000, e->s->nr_parts);
printParticle(e->s->parts, 515050, e->s->nr_parts);
exit(0);
}
/**
......
src/hydro/Default/hydro_part.h
View file @ a3596a2e
......@@ -73,7 +73,7 @@ struct part {
float alpha;
/* Store density/force specific stuff. */
union {
//union {
struct {
......@@ -112,7 +112,7 @@ struct part {
float c;
} force;
};
//};
/* Particle mass. */
float mass;
......
src/hydro/Gadget2/hydro.h
View file @ a3596a2e
......@@ -28,19 +28,20 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
struct part* p, struct xpart* xp) {
/* CFL condition */
float dt_cfl = const_cfl * p->h / p->force.v_sig;
const float dt_cfl = const_cfl * p->h / p->v_sig;
/* Limit change in h */
float dt_h_change = (p->force.h_dt != 0.0f)
? fabsf(const_ln_max_h_change * p->h / p->force.h_dt)
: FLT_MAX;
/* /\* Limit change in h *\/ */
/* float dt_h_change = (p->h_dt != 0.0f) */
/* ? fabsf(const_ln_max_h_change * p->h / p->h_dt) */
/* : FLT_MAX; */
/* Limit change in u */
float dt_u_change = (p->force.u_dt != 0.0f)
? fabsf(const_max_u_change * p->u / p->force.u_dt)
: FLT_MAX;
/* /\* Limit change in u *\/ */
/* float dt_u_change = (p->force.u_dt != 0.0f) */
/* ? fabsf(const_max_u_change * p->u / p->force.u_dt) */
/* : FLT_MAX; */
return fminf(dt_cfl, fminf(dt_h_change, dt_u_change));
// return fminf(dt_cfl, fminf(dt_h_change, dt_u_change));
return dt_cfl;
}
......@@ -58,10 +59,11 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(struct part* p
p->density.wcount_dh = 0.f;
p->rho = 0.f;
p->rho_dh = 0.f;
p->density.div_v = 0.f;
p->density.curl_v[0] = 0.f;
p->density.curl_v[1] = 0.f;
p->density.curl_v[2] = 0.f;
p->div_v = 0.f;
p->curl_v = 0.f;
p->rot_v[0] = 0.f;
p->rot_v[1] = 0.f;
p->rot_v[2] = 0.f;
}
/**
......@@ -71,8 +73,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(struct part* p
* and add the self-contribution term.
*
* @param p The particle to act upon
* @param time The current time
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(struct part* p) {
__attribute__((always_inline)) INLINE static void hydro_end_density(struct part* p, float time) {
/* Some smoothing length multiples. */
const float h = p->h;
......@@ -80,12 +83,31 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(struct part*
const float ih2 = ih * ih;
const float ih4 = ih2 * ih2;
/* Final operation on the density. */
p->rho = ih * ih2 * (p->rho + p->mass * kernel_root);
p->rho_dh = (p->rho_dh - 3.0f * p->mass * kernel_root) * ih4;
p->density.wcount = (p->density.wcount + kernel_root) *
(4.0f / 3.0 * M_PI * kernel_gamma3);
p->density.wcount_dh = p->density.wcount_dh * ih * (4.0f / 3.0 * M_PI * kernel_gamma3);
/* Final operation on the density (add self-contribution). */
p->rho += p->mass * kernel_root;
p->rho_dh -= 3.0f * p->mass * kernel_root * kernel_igamma;
p->density.wcount += kernel_root;
/* Finish the calculation by inserting the missing h-factors */
p->rho *= ih * ih2;
p->rho_dh *= ih4;
p->density.wcount *= (4.0f / 3.0 * M_PI * kernel_gamma3);
p->density.wcount_dh *= ih * (4.0f / 3.0 * M_PI * kernel_gamma3);
/* Compute the derivative term */
p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh / p->rho);
/* Finish calculation of the velocity curl */
p->curl_v = sqrtf(p->rot_v[0] * p->rot_v[0] +
p->rot_v[1] * p->rot_v[1] +
p->rot_v[2] * p->rot_v[2]) / p->rho;
/* Finish calculation of the velocity divergence */
p->div_v /= p->rho;
/* Compute the pressure */
const float dt = time - 0.5f*(p->t_begin + p->t_end);
p->pressure = (p->entropy + p->entropy_dt * dt) * pow(p->rho, const_hydro_gamma);
}
......@@ -102,31 +124,31 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(struct par
struct xpart* xp) {
/* Some smoothing length multiples. */
const float h = p->h;
const float ih = 1.0f / h;
const float ih2 = ih * ih;
const float ih4 = ih2 * ih2;
/* const float h = p->h; */
/* const float ih = 1.0f / h; */
/* const float ih2 = ih * ih; */
/* const float ih4 = ih2 * ih2; */
/* Pre-compute some stuff for the balsara switch. */
const float normDiv_v = fabs(p->density.div_v / p->rho * ih4);
const float normCurl_v = sqrtf(p->density.curl_v[0] * p->density.curl_v[0] +
p->density.curl_v[1] * p->density.curl_v[1] +
p->density.curl_v[2] * p->density.curl_v[2]) /
p->rho * ih4;
/* Compute this particle's sound speed. */
const float u = p->u;
const float fc = p->force.c =
sqrtf(const_hydro_gamma * (const_hydro_gamma - 1.0f) * u);
/* /\* Pre-compute some stuff for the balsara switch. *\/ */
/* const float normDiv_v = fabs(p->density.div_v / p->rho * ih4); */
/* const float normCurl_v = sqrtf(p->density.curl_v[0] * p->density.curl_v[0] + */
/* p->density.curl_v[1] * p->density.curl_v[1] + */
/* p->density.curl_v[2] * p->density.curl_v[2]) / */
/* p->rho * ih4; */
/* /\* Compute this particle's sound speed. *\/ */
/* const float u = p->u; */
/* const float fc = p->force.c = */
/* sqrtf(const_hydro_gamma * (const_hydro_gamma - 1.0f) * u); */
/* Compute the P/Omega/rho2. */
xp->omega = 1.0f + 0.3333333333f * h * p->rho_dh / p->rho;
p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (p->rho * xp->omega);
/* /\* Compute the P/Omega/rho2. *\/ */
/* xp->omega = 1.0f + 0.3333333333f * h * p->rho_dh / p->rho; */
/* p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (p->rho * xp->omega); */
/* Balsara switch */
p->force.balsara =
normDiv_v / (normDiv_v + normCurl_v + 0.0001f * fc * ih);
/* /\* Balsara switch *\/ */
/* p->force.balsara = */
/* normDiv_v / (normDiv_v + normCurl_v + 0.0001f * fc * ih); */
}
......@@ -147,9 +169,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(struc
p->a[2] = 0.0f;
/* Reset the time derivatives. */
p->force.u_dt = 0.0f;
p->force.h_dt = 0.0f;
p->force.v_sig = 0.0f;
p->v_sig = 0.0f;
}
......@@ -163,19 +183,6 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(struc
__attribute__((always_inline)) INLINE static void hydro_predict_extra(struct part* p,
struct xpart* xp,
float dt) {
float u, w;
/* Predict internal energy */
w = p->force.u_dt / p->u * dt;
if (fabsf(w) < 0.01f) /* 1st order expansion of exp(w) */
u = p->u *=
1.0f +
w * (1.0f + w * (0.5f + w * (1.0f / 6.0f + 1.0f / 24.0f * w)));
else
u = p->u *= expf(w);
/* Predict gradient term */
p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (p->rho * xp->omega);
}
......@@ -190,6 +197,17 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(struct par
*/
__attribute__((always_inline)) INLINE static void hydro_end_force(struct part* p) {
p->force.h_dt *= p->h * 0.333333333f;
}
/**
* @brief Converts hydro quantity of a particle
*
* Requires the density to be known
*
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(struct part* p) {
p->entropy = (const_hydro_gamma - 1.f) * p->entropy * pow(p->rho, -(const_hydro_gamma - 1.f));
}
src/hydro/Gadget2/hydro_iact.h
View file @ a3596a2e
This diff is collapsed.
src/hydro/Gadget2/hydro_part.h
View file @ a3596a2e
......@@ -59,8 +59,18 @@ struct part {
/* Particle time of end of time-step. */
float t_end;
/* Particle internal energy. */
float u;
struct {
/* Number of neighbours */
float wcount;
/* Number of neighbours spatial derivative */
float wcount_dh;
} density;
/* Particle entropy. */
float entropy;
/* Particle density. */
float rho;
......@@ -69,52 +79,30 @@ struct part {
*/
float rho_dh;
/* Store density/force specific stuff. */
union {
struct {
/* Particle velocity divergence. */
float div_v;
/* Derivative of particle number density. */
float wcount_dh;
/* Particle velocity curl. */
float curl_v[3];
/* Particle number density. */
float wcount;
} density;
struct {
/* Balsara switch */
float balsara;
/* Aggregate quantities. */
float POrho2;
/* Change in particle energy over time. */
float u_dt;
/* Change in smoothing length over time. */
float h_dt;
/* Signal velocity */
float v_sig;
/* Sound speed */
float c;
} force;
};
/* Particle mass. */
float mass;
/* Particle pressure */
float pressure;
/* Entropy time derivative */
float entropy_dt;
/* Velocity divergence */
float div_v;
/* Velocity curl */
float curl_v;
float rot_v[3];
/* Signal velocity */
float v_sig;
/* temp */
struct {
float h_dt; //MATTHIEU
} force;
/* Particle ID. */
unsigned long long id;
......
src/runner.c
View file @ a3596a2e
......@@ -36,6 +36,7 @@
/* Local headers. */
#include "atomic.h"
#include "const.h"
#include "debug.h"
#include "engine.h"
#include "error.h"
#include "scheduler.h"
......@@ -587,7 +588,7 @@ void runner_doghost(struct runner *r, struct cell *c) {
if (p->t_end <= t_end) {
/* Finish the density calculation */
hydro_end_density(p);
hydro_end_density(p, t_end);
/* If no derivative, double the smoothing length. */
if (p->density.wcount_dh == 0.0f) h_corr = p->h;
......@@ -620,7 +621,9 @@ void runner_doghost(struct runner *r, struct cell *c) {
}
/* We now have a particle whose smoothing length has converged */
//if(p->id == 1000)
// printParticle(parts, 1000, count);
/* As of here, particle force variables will be set. Do _NOT_
try to read any particle density variables! */
......
src/runner_doiact.h
View file @ a3596a2e
......@@ -1400,6 +1400,9 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
/* Get a pointer to the ith particle. */
pi = &parts[pid];
if(pi->id == 1000) message("oO 1000");
if(pi->id == 515050) message("oO 515050");
/* Get the particle position and radius. */
for (k = 0; k < 3; k++) pix[k] = pi->x[k];
hi = pi->h;
......@@ -1624,6 +1627,10 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
/* Get a pointer to the ith particle. */
pi = &parts[pid];
if(pi->id == 1000) message("oO 1000");
if(pi->id == 515050) message("oO 515050");
/* Get the particle position and radius. */
for (k = 0; k < 3; k++) pix[k] = pi->x[k];
hi = pi->h;
......
src/single_io.c
View file @ a3596a2e
......@@ -227,7 +227,7 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
readArray(h_grp, "Velocities", FLOAT, *N, 3, *parts, v, COMPULSORY);
readArray(h_grp, "Masses", FLOAT, *N, 1, *parts, mass, COMPULSORY);
readArray(h_grp, "SmoothingLength", FLOAT, *N, 1, *parts, h, COMPULSORY);
readArray(h_grp, "InternalEnergy", FLOAT, *N, 1, *parts, u, COMPULSORY);
readArray(h_grp, "InternalEnergy", FLOAT, *N, 1, *parts, entropy, COMPULSORY); //MATTHIEU
readArray(h_grp, "ParticleIDs", ULONGLONG, *N, 1, *parts, id, COMPULSORY);
// readArray(h_grp, "TimeStep", FLOAT, *N, 1, *parts, dt, OPTIONAL);
readArray(h_grp, "Acceleration", FLOAT, *N, 3, *parts, a, OPTIONAL);
......@@ -469,8 +469,8 @@ void write_output_single(struct engine* e, struct UnitSystem* us) {
UNIT_CONV_MASS);
writeArray(h_grp, fileName, xmfFile, "SmoothingLength", FLOAT, N, 1, parts, h,
us, UNIT_CONV_LENGTH);
writeArray(h_grp, fileName, xmfFile, "InternalEnergy", FLOAT, N, 1, parts, u,
us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
writeArray(h_grp, fileName, xmfFile, "InternalEnergy", FLOAT, N, 1, parts, entropy,
us, UNIT_CONV_ENERGY_PER_UNIT_MASS); //MATTHIEU
writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N, 1, parts,
id, us, UNIT_CONV_NO_UNITS);
/* writeArray(h_grp, fileName, xmfFile, "TimeStep", FLOAT, N, 1, parts, dt,
......
src/space.c
View file @ a3596a2e
......@@ -1198,7 +1198,7 @@ void space_init(struct space *s, double dim[3], struct part *parts, int N,
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
xp->u_hdt = p->u;
//xp->u_hdt = p->u; // MATTHIEU
}
/* For now, clone the parts to make gparts. */
......
src/tools.c
View file @ a3596a2e
......@@ -300,20 +300,22 @@ void density_dump(int N) {
int k;
float r2[4] = {0.0f, 0.0f, 0.0f, 0.0f}, hi[4], hj[4];
struct part *pi[4], *pj[4], Pi[4], Pj[4];
struct part /**pi[4], *pj[4],*/ Pi[4], Pj[4];
/* Init the interaction parameters. */
for (k = 0; k < 4; k++) {
Pi[k].mass = 1.0f;
Pi[k].rho = 0.0f;
Pi[k].density.wcount = 0.0f;
Pi[k].id = k;
Pj[k].mass = 1.0f;
Pj[k].rho = 0.0f;
Pj[k].density.wcount = 0.0f;
Pj[k].id = k+4;
hi[k] = 1.0;
hj[k] = 1.0;
pi[k] = &Pi[k];
pj[k] = &Pj[k];
//pi[k] = &Pi[k];
//pj[k] = &Pj[k];
}
for (k = 0; k <= N; k++) {
......@@ -325,17 +327,19 @@ void density_dump(int N) {
Pj[0].density.wcount = 0;
runner_iact_density(r2[0], NULL, hi[0], hj[0], &Pi[0], &Pj[0]);
printf(" %e %e %e", r2[0], Pi[0].density.wcount, Pj[0].density.wcount);
Pi[0].density.wcount = 0;
Pj[0].density.wcount = 0;
Pi[1].density.wcount = 0;
Pj[1].density.wcount = 0;
Pi[2].density.wcount = 0;
Pj[2].density.wcount = 0;
Pi[3].density.wcount = 0;