Commit 133171b2 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

And OK for the MINIMAL_SPH

parent 08c53168
......@@ -56,8 +56,8 @@
/* SPH variant to use */
//#define MINIMAL_SPH
//#define GADGET2_SPH
#define DEFAULT_SPH
#define GADGET2_SPH
//#define DEFAULT_SPH
/* External potential properties */
#define EXTERNAL_POTENTIAL_POINTMASS
......
......@@ -246,3 +246,39 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
return p->u;
}
/**
* @brief Returns the pressure of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
const struct part *restrict p, float dt) {
return p->u * p->rho * hydro_gamma_minus_one;
}
/**
* @brief Returns the entropy of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
const struct part *restrict p, float dt) {
return hydro_gamma_minus_one * p->u * pow_minus_gamma_minus_one(p->rho);
}
/**
* @brief Returns the sound speed of a particle
*
* @param p The particle of interest
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part *restrict p, float dt) {
return sqrt(p->u * hydro_gamma * hydro_gamma_minus_one);
}
......@@ -97,6 +97,8 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
part->entropy = pressure / pow_gamma(density);
#elif defined(DEFAULT_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(MINIMAL_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
#else
error("Need to define pressure here !");
#endif
......@@ -298,7 +300,12 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->parts[pid].x[1], main_cell->parts[pid].x[2],
main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
main_cell->parts[pid].v[2], main_cell->parts[pid].h,
main_cell->parts[pid].rho, main_cell->parts[pid].density.div_v,
main_cell->parts[pid].rho,
#ifdef MINIMAL_SPH
0.f,
#else
main_cell->parts[pid].density.div_v,
#endif
hydro_get_entropy(&main_cell->parts[pid], 0.f),
hydro_get_internal_energy(&main_cell->parts[pid], 0.f),
hydro_get_pressure(&main_cell->parts[pid], 0.f),
......@@ -310,6 +317,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->parts[pid].entropy_dt, 0.f
#elif defined(DEFAULT_SPH)
0.f, main_cell->parts[pid].force.u_dt
#elif defined(MINIMAL_SPH)
0.f, main_cell->parts[pid].u_dt
#else
0.f, 0.f
#endif
......
......@@ -78,8 +78,14 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
"%6llu %10f %10f %10f %10f %10f %10f %13e %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->rho,
p->rho_dh, p->density.wcount, p->density.wcount_dh, p->density.div_v,
p->density.rot_v[0], p->density.rot_v[1], p->density.rot_v[2]);
p->rho_dh, p->density.wcount, p->density.wcount_dh,
#ifdef MINIMAL_SPH
0.f, 0.f, 0.f, 0.f
#else
p->density.div_v,
p->density.rot_v[0], p->density.rot_v[1], p->density.rot_v[2]
#endif
);
fclose(file);
}
......
......@@ -75,8 +75,14 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
"%6llu %10f %10f %10f %10f %10f %10f %13e %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->rho,
p->rho_dh, p->density.wcount, p->density.wcount_dh, p->density.div_v,
p->density.rot_v[0], p->density.rot_v[1], p->density.rot_v[2]);
p->rho_dh, p->density.wcount, p->density.wcount_dh,
#ifdef MINIMAL_SPH
0.f, 0.f, 0.f, 0.f
#else
p->density.div_v,
p->density.rot_v[0], p->density.rot_v[1], p->density.rot_v[2]
#endif
);
fclose(file);
}
......
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