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

Cleaned up Gizmo/hydro.h and added some extra documentation.

parent 911dbb35
Branches
Tags
1 merge request!317Cleaned up GIZMO code, added SineWavePotential tests.
Hide whitespace changes
Inline Side-by-side
Showing
with 30 additions and 42 deletions
src/hydro/Gizmo/hydro.h
+ 30
− 42
View file @ cf57807a
...@@ -111,27 +111,32 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part( ...@@ -111,27 +111,32 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->conserved.momentum[2] = mass * p->primitives.v[2]; p->conserved.momentum[2] = mass * p->primitives.v[2];
/* and the thermal energy */ /* and the thermal energy */
/* remember that we store the total thermal energy, not the specific thermal
energy (as in Gadget) */
#if defined(EOS_ISOTHERMAL_GAS) #if defined(EOS_ISOTHERMAL_GAS)
/* this overwrites the internal energy from the initial condition file */
p->conserved.energy = mass * const_isothermal_internal_energy; p->conserved.energy = mass * const_isothermal_internal_energy;
#else #else
p->conserved.energy *= mass; p->conserved.energy *= mass;
#endif #endif
#ifdef GIZMO_TOTAL_ENERGY #ifdef GIZMO_TOTAL_ENERGY
/* add the total kinetic energy */
p->conserved.energy += 0.5f * (p->conserved.momentum[0] * p->primitives.v[0] + p->conserved.energy += 0.5f * (p->conserved.momentum[0] * p->primitives.v[0] +
p->conserved.momentum[1] * p->primitives.v[1] + p->conserved.momentum[1] * p->primitives.v[1] +
p->conserved.momentum[2] * p->primitives.v[2]); p->conserved.momentum[2] * p->primitives.v[2]);
#endif #endif
#if defined(GIZMO_FIX_PARTICLES) #if defined(GIZMO_FIX_PARTICLES)
/* make sure the particles are initially at rest */
p->v[0] = 0.; p->v[0] = 0.;
p->v[1] = 0.; p->v[1] = 0.;
p->v[2] = 0.; p->v[2] = 0.;
#else #endif
xp->v_full[0] = p->v[0]; xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1]; xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2]; xp->v_full[2] = p->v[2];
#endif
} }
/** /**
...@@ -240,6 +245,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( ...@@ -240,6 +245,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->primitives.v[2] = momentum[2] / m; p->primitives.v[2] = momentum[2] / m;
#ifdef EOS_ISOTHERMAL_GAS #ifdef EOS_ISOTHERMAL_GAS
/* although the pressure is not formally used anywhere if an isothermal eos
has been selected, we still make sure it is set to the correct value */
p->primitives.P = const_isothermal_soundspeed * const_isothermal_soundspeed * p->primitives.P = const_isothermal_soundspeed * const_isothermal_soundspeed *
p->primitives.rho; p->primitives.rho;
#else #else
...@@ -247,15 +254,20 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( ...@@ -247,15 +254,20 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
float energy = p->conserved.energy; float energy = p->conserved.energy;
#ifdef GIZMO_TOTAL_ENERGY #ifdef GIZMO_TOTAL_ENERGY
/* subtract the kinetic energy; we want the thermal energy */
energy -= 0.5f * (momentum[0] * p->primitives.v[0] + energy -= 0.5f * (momentum[0] * p->primitives.v[0] +
momentum[1] * p->primitives.v[1] + momentum[1] * p->primitives.v[1] +
momentum[2] * p->primitives.v[2]); momentum[2] * p->primitives.v[2]);
#endif #endif
/* energy contains the total thermal energy, we want the specific energy.
this is why we divide by the volume, and not by the density */
p->primitives.P = hydro_gamma_minus_one * energy / volume; p->primitives.P = hydro_gamma_minus_one * energy / volume;
#endif #endif
/* sanity checks */ /* sanity checks */
/* it would probably be safer to throw a warning if netive densities or
pressures occur */
if (p->primitives.rho < 0.0f || p->primitives.P < 0.0f) { if (p->primitives.rho < 0.0f || p->primitives.P < 0.0f) {
p->primitives.rho = 0.0f; p->primitives.rho = 0.0f;
p->primitives.P = 0.0f; p->primitives.P = 0.0f;
...@@ -285,6 +297,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -285,6 +297,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->timestepvars.vmax = 0.0f; p->timestepvars.vmax = 0.0f;
/* Set the actual velocity of the particle */ /* Set the actual velocity of the particle */
/* if GIZMO_FIX_PARTICLES has been selected, v_full will always be zero */
p->force.v_full[0] = xp->v_full[0]; p->force.v_full[0] = xp->v_full[0];
p->force.v_full[1] = xp->v_full[1]; p->force.v_full[1] = xp->v_full[1];
p->force.v_full[2] = xp->v_full[2]; p->force.v_full[2] = xp->v_full[2];
...@@ -375,6 +388,9 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( ...@@ -375,6 +388,9 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
else else
p->h *= expf(w1); p->h *= expf(w1);
/* we temporarily disabled the primitive variable drift.
This should be reenabled once gravity works, and we have time to check that
drifting works properly. */
// const float w2 = -hydro_dimension * w1; // const float w2 = -hydro_dimension * w1;
// if (fabsf(w2) < 0.2f) { // if (fabsf(w2) < 0.2f) {
// p->primitives.rho *= approx_expf(w2); // p->primitives.rho *= approx_expf(w2);
...@@ -416,56 +432,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( ...@@ -416,56 +432,18 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
__attribute__((always_inline)) INLINE static void hydro_end_force( __attribute__((always_inline)) INLINE static void hydro_end_force(
struct part* p) { struct part* p) {
/* set the variables that are used to drift the primitive variables */
/* Add normalization to h_dt. */ /* Add normalization to h_dt. */
p->force.h_dt *= p->h * hydro_dimension_inv; p->force.h_dt *= p->h * hydro_dimension_inv;
/* Set the hydro acceleration, based on the new momentum and mass */
/* NOTE: the momentum and mass are only correct for active particles, since
only active particles have received flux contributions from all their
neighbours. Since this method is only called for active particles,
this is indeed the case. */
if (p->force.dt) { if (p->force.dt) {
// float mnew;
// float vnew[3];
// mnew = p->conserved.mass + p->conserved.flux.mass;
// if (mnew > 0.) {
// vnew[0] =
// (p->conserved.momentum[0] + p->conserved.flux.momentum[0]) /
// mnew;
// vnew[1] =
// (p->conserved.momentum[1] + p->conserved.flux.momentum[1]) /
// mnew;
// vnew[2] =
// (p->conserved.momentum[2] + p->conserved.flux.momentum[2]) /
// mnew;
// } else {
// vnew[0] = 0.;
// vnew[1] = 0.;
// vnew[2] = 0.;
// }
// p->a_hydro[0] = (vnew[0] - p->force.v_full[0]) / p->force.dt;
// p->a_hydro[1] = (vnew[1] - p->force.v_full[1]) / p->force.dt;
// p->a_hydro[2] = (vnew[2] - p->force.v_full[2]) / p->force.dt;
p->a_hydro[0] = 0.;
p->a_hydro[1] = 0.;
p->a_hydro[2] = 0.;
p->du_dt = p->conserved.flux.energy / p->force.dt; p->du_dt = p->conserved.flux.energy / p->force.dt;
} else { } else {
p->a_hydro[0] = 0.0f;
p->a_hydro[1] = 0.0f;
p->a_hydro[2] = 0.0f;
p->du_dt = 0.0f; p->du_dt = 0.0f;
} }
#if defined(GIZMO_FIX_PARTICLES) #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; p->du_dt = 0.0f;
/* disable the smoothing length update, since the smoothing lengths should /* disable the smoothing length update, since the smoothing lengths should
...@@ -510,6 +488,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -510,6 +488,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
a_grav[2] = p->gpart->a_grav[2]; a_grav[2] = p->gpart->a_grav[2];
/* Store the gravitational acceleration for later use. */ /* Store the gravitational acceleration for later use. */
/* This is currently only used for output purposes. */
p->gravity.old_a[0] = a_grav[0]; p->gravity.old_a[0] = a_grav[0];
p->gravity.old_a[1] = a_grav[1]; p->gravity.old_a[1] = a_grav[1];
p->gravity.old_a[2] = a_grav[2]; p->gravity.old_a[2] = a_grav[2];
...@@ -518,11 +497,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -518,11 +497,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->gpart->mass = p->conserved.mass; p->gpart->mass = p->conserved.mass;
/* Kick the momentum for half a time step */ /* Kick the momentum for half a time step */
/* Note that this also affects the particle movement, as the velocity for
the particles is set after this. */
p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0]; p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0];
p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1]; p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2]; p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
#if !defined(EOS_ISOTHERMAL_GAS) && defined(GIZMO_TOTAL_ENERGY) #if !defined(EOS_ISOTHERMAL_GAS) && defined(GIZMO_TOTAL_ENERGY)
/* This part still needs to be tested! */
p->conserved.energy += dt * (p->conserved.momentum[0] * a_grav[0] + p->conserved.energy += dt * (p->conserved.momentum[0] * a_grav[0] +
p->conserved.momentum[1] * a_grav[1] + p->conserved.momentum[1] * a_grav[1] +
p->conserved.momentum[2] * a_grav[2]); p->conserved.momentum[2] * a_grav[2]);
...@@ -563,6 +545,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -563,6 +545,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
xp->v_full[1] = p->conserved.momentum[1] / p->conserved.mass; xp->v_full[1] = p->conserved.momentum[1] / p->conserved.mass;
xp->v_full[2] = p->conserved.momentum[2] / p->conserved.mass; xp->v_full[2] = p->conserved.momentum[2] / p->conserved.mass;
} else { } else {
/* vacuum particles don't move */
xp->v_full[0] = 0.; xp->v_full[0] = 0.;
xp->v_full[1] = 0.; xp->v_full[1] = 0.;
xp->v_full[2] = 0.; xp->v_full[2] = 0.;
...@@ -572,6 +555,9 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -572,6 +555,9 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->v[2] = xp->v_full[2]; p->v[2] = xp->v_full[2];
/* Update gpart! */ /* Update gpart! */
/* This is essential, as the gpart drift is done independently from the part
drift, and we don't want the gpart and the part to have different
positions! */
if (p->gpart) { if (p->gpart) {
p->gpart->v_full[0] = xp->v_full[0]; p->gpart->v_full[0] = xp->v_full[0];
p->gpart->v_full[1] = xp->v_full[1]; p->gpart->v_full[1] = xp->v_full[1];
...@@ -684,6 +670,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy( ...@@ -684,6 +670,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
energy (u*m) */ energy (u*m) */
p->conserved.energy = u * p->conserved.mass; p->conserved.energy = u * p->conserved.mass;
#ifdef GIZMO_TOTAL_ENERGY #ifdef GIZMO_TOTAL_ENERGY
/* add the kinetic energy */
p->conserved.energy += 0.5f * p->conserved.mass * p->conserved.energy += 0.5f * p->conserved.mass *
(p->conserved.momentum[0] * p->primitives.v[0] + (p->conserved.momentum[0] * p->primitives.v[0] +
p->conserved.momentum[1] * p->primitives.v[1] + p->conserved.momentum[1] * p->primitives.v[1] +
...@@ -707,6 +694,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy( ...@@ -707,6 +694,7 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
p->conserved.energy = S * pow_gamma_minus_one(p->primitives.rho) * p->conserved.energy = S * pow_gamma_minus_one(p->primitives.rho) *
hydro_one_over_gamma_minus_one * p->conserved.mass; hydro_one_over_gamma_minus_one * p->conserved.mass;
#ifdef GIZMO_TOTAL_ENERGY #ifdef GIZMO_TOTAL_ENERGY
/* add the kinetic energy */
p->conserved.energy += 0.5f * p->conserved.mass * p->conserved.energy += 0.5f * p->conserved.mass *
(p->conserved.momentum[0] * p->primitives.v[0] + (p->conserved.momentum[0] * p->primitives.v[0] +
p->conserved.momentum[1] * p->primitives.v[1] + p->conserved.momentum[1] * p->primitives.v[1] +
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment