Commit 35fd53ea authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Made a small but important change in computation of hydro acceleration. Added...

Made a small but important change in computation of hydro acceleration. Added total energy to snapshots instead of thermal energy.
parent a79c0e06
......@@ -343,13 +343,13 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
this is indeed the case. */
if (p->force.dt) {
p->a_hydro[0] =
(p->conserved.momentum[0] / p->conserved.mass - p->primitives.v[0]) /
(p->conserved.momentum[0] / p->conserved.mass - p->force.v_full[0]) /
p->force.dt;
p->a_hydro[1] =
(p->conserved.momentum[1] / p->conserved.mass - p->primitives.v[1]) /
(p->conserved.momentum[1] / p->conserved.mass - p->force.v_full[1]) /
p->force.dt;
p->a_hydro[2] =
(p->conserved.momentum[2] / p->conserved.mass - p->primitives.v[2]) /
(p->conserved.momentum[2] / p->conserved.mass - p->force.v_full[2]) /
p->force.dt;
p->du_dt = p->conserved.flux.energy / p->force.dt;
......
......@@ -55,14 +55,45 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, primitives.rho);
}
/**
* @brief Get the internal energy of a particle
*
* @param e #engine.
* @param p Particle.
* @return Internal energy of the particle
*/
float convert_u(struct engine* e, struct part* p) {
return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
}
/**
* @brief Get the entropic function of a particle
*
* @param e #engine.
* @param p Particle.
* @return Entropic function of the particle
*/
float convert_A(struct engine* e, struct part* p) {
return p->primitives.P / pow_gamma(p->primitives.rho);
}
/**
* @brief Get the total energy of a particle
*
* @param e #engine.
* @param p Particle.
* @return Total energy of the particle
*/
float convert_Etot(struct engine* e, struct part* p) {
float momentum2;
momentum2 = p->conserved.momentum[0] * p->conserved.momentum[0] +
p->conserved.momentum[1] * p->conserved.momentum[1] +
p->conserved.momentum[2] * p->conserved.momentum[2];
return p->conserved.energy + 0.5f * momentum2 / p->conserved.mass;
}
/**
* @brief Specifies which particle fields to write to a dataset
*
......@@ -101,8 +132,9 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
"Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, primitives.P, convert_A);
list[11] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
parts, primitives.P);
list[12] = io_make_output_field("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
parts, conserved.energy);
list[12] =
io_make_output_field_convert_part("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
parts, conserved.energy, convert_Etot);
}
/**
......
......@@ -48,10 +48,10 @@
/* Task type names. */
const char *taskID_names[task_type_count] = {
"none", "sort", "self", "pair", "sub_self",
"sub_pair", "init", "ghost", "extra_ghost",
"kick", "kick_fixdt", "send", "recv", "grav_gather_m",
"grav_fft", "grav_mm", "grav_up", "grav_external"};
"none", "sort", "self", "pair", "sub_self",
"sub_pair", "init", "ghost", "extra_ghost", "kick",
"kick_fixdt", "send", "recv", "grav_gather_m", "grav_fft",
"grav_mm", "grav_up", "grav_external"};
const char *subtaskID_names[task_subtype_count] = {
"none", "density", "gradient", "force", "grav", "tend"};
......
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