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

Tweaked entropy switch parameters for GizmoMFM. It works, but results are not great.

parent 332e9d44
Branches gizmo
No related tags found
No related merge requests found
......@@ -697,17 +697,20 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.momentum[2] += dt_grav * p->conserved.mass * a_grav[2];
/* apply the entropy switch */
const float psize = powf(p->geometry.volume / hydro_dimension_unit_sphere,
hydro_dimension_inv);
const float dEgrav = p->conserved.mass *
sqrtf(a_grav[0] * a_grav[0] + a_grav[1] * a_grav[1] +
a_grav[2] * a_grav[2]) *
cosmo->a_inv * p->h;
if (p->conserved.energy < const_entropy_switch_ekin_fac *
(p->force.Ekinmax + p->conserved.energy) ||
cosmo->a_inv * psize;
const float dEkin = 0.5f * p->conserved.mass * p->force.Ekinmax;
if (p->conserved.energy <
const_entropy_switch_ekin_fac * (dEkin + p->conserved.energy) ||
p->conserved.energy < const_entropy_switch_grav_fac * dEgrav) {
p->conserved.energy = hydro_one_over_gamma_minus_one *
p->conserved.entropy * pow_gamma_minus_one(p->rho);
if (p->conserved.energy < const_entropy_switch_ekin_fac *
(p->force.Ekinmax + p->conserved.energy)) {
if (p->conserved.energy <
const_entropy_switch_ekin_fac * (dEkin + p->conserved.energy)) {
hydro_add_flag(p, GIZMO_FLAG_EKIN_SWITCH);
}
if (p->conserved.energy < const_entropy_switch_grav_fac * dEgrav) {
......
......@@ -38,7 +38,7 @@ enum gizmo_flag_variables {
* @param p Particle.
*/
__attribute__((always_inline)) INLINE static void hydro_flag_variable_init(
struct part *p) {
struct part* p) {
p->flag_variable = 0;
}
......@@ -50,9 +50,10 @@ __attribute__((always_inline)) INLINE static void hydro_flag_variable_init(
* @param num_fields The number of i/o fields to write.
*/
__attribute__((always_inline)) INLINE static void hydro_write_flag_variable(
const struct part* parts, struct io_props* list, int* num_fields) {
const struct part* parts, struct io_props* list, int* num_fields) {
list[*num_fields] = io_make_output_field("FlagVariable", INT, 1, UNIT_CONV_NO_UNITS, parts, flag_variable);
list[*num_fields] = io_make_output_field(
"FlagVariable", INT, 1, UNIT_CONV_NO_UNITS, parts, flag_variable);
++(*num_fields);
}
......@@ -62,8 +63,8 @@ __attribute__((always_inline)) INLINE static void hydro_write_flag_variable(
* @param p Particle.
* @param value Flag value to add.
*/
__attribute__((always_inline)) INLINE static void hydro_add_flag(
struct part* p, int value) {
__attribute__((always_inline)) INLINE static void hydro_add_flag(struct part* p,
int value) {
p->flag_variable |= (1 << value);
}
......@@ -76,7 +77,7 @@ __attribute__((always_inline)) INLINE static void hydro_add_flag(
* @param p Particle.
*/
__attribute__((always_inline)) INLINE static void hydro_flag_variable_init(
struct part *p) {}
struct part* p) {}
/**
* @brief Add the flag variable to the snapshot output.
......@@ -86,7 +87,7 @@ __attribute__((always_inline)) INLINE static void hydro_flag_variable_init(
* @param num_fields The number of i/o fields to write.
*/
__attribute__((always_inline)) INLINE static void hydro_write_flag_variable(
const struct part* parts, struct io_props* list, int* num_fields) {}
const struct part* parts, struct io_props* list, int* num_fields) {}
/**
* @brief Add the given flag value to the flag variable.
......@@ -94,8 +95,8 @@ __attribute__((always_inline)) INLINE static void hydro_write_flag_variable(
* @param p Particle.
* @param value Flag value to add.
*/
__attribute__((always_inline)) INLINE static void hydro_add_flag(
struct part* p, int value) {}
__attribute__((always_inline)) INLINE static void hydro_add_flag(struct part* p,
int value) {}
#endif
......
......@@ -261,17 +261,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
vmax = 0.0f;
/* calculate flux limiter quantities */
const float Wjminvi[3] = {Wj[1] - vi[0], Wj[2] - vi[1], Wj[3] - vi[2]};
const float Ekinngbj = 0.5f * pj->conserved.mass *
(Wjminvi[0] * Wjminvi[0] + Wjminvi[1] * Wjminvi[1] +
Wjminvi[2] * Wjminvi[2]);
const float Wjminvi[3] = {Wj[1] - Wi[1], Wj[2] - Wi[2], Wj[3] - Wi[3]};
const float Ekinngbj = Wjminvi[0] * Wjminvi[0] + Wjminvi[1] * Wjminvi[1] +
Wjminvi[2] * Wjminvi[2];
pi->force.Ekinmax = max(pi->force.Ekinmax, Ekinngbj);
if (mode == 1) {
const float Wiminvj[3] = {Wi[1] - vj[0], Wi[2] - vj[1], Wi[3] - vj[2]};
const float Ekinngbi = 0.5f * pi->conserved.mass *
(Wiminvj[0] * Wiminvj[0] + Wiminvj[1] * Wiminvj[1] +
Wiminvj[2] * Wiminvj[2]);
pj->force.Ekinmax = max(pj->force.Ekinmax, Ekinngbi);
pj->force.Ekinmax = max(pj->force.Ekinmax, Ekinngbj);
}
/* Velocity on the axis linking the particles */
......
......@@ -53,10 +53,10 @@
#define const_viscosity_beta 3.0f
/* Prefactor for the kinetic energy condition for the entropy switch. */
#define const_entropy_switch_ekin_fac 0.0001f
#define const_entropy_switch_ekin_fac 0.000001f
/* Prefactor for the gravitational energy condition for the entropy switch. */
#define const_entropy_switch_grav_fac 0.0001f
#define const_entropy_switch_grav_fac 0.000001f
/*! Activate this to write a diagnostic flag variable to the snapshots. */
#define GIZMO_FLAG_VARIABLE
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment