Commit c89b6aab authored by Jacob Kegerreis's avatar Jacob Kegerreis
Browse files

Set Balsara as the default for Planetary hydro

parent cfa8ac1b
......@@ -26,8 +26,8 @@
* equations) with multiple materials.
*
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
* viscosity term with the Balsara (1995) switch (optional).
* No thermal conduction term is implemented.
*
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with
* \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational
......@@ -45,9 +45,9 @@
#include "minmax.h"
/*
* Note: Define PLANETARY_SPH_BALSARA to use the Balsara (1995) switch for
* the artificial viscosity, instead of the default Monaghan (1992).
* i.e. compile with: make CFLAGS=-DPLANETARY_SPH_BALSARA to use.
* Note: Define PLANETARY_SPH_NO_BALSARA to disable the Balsara (1995) switch
* for the artificial viscosity and use the vanilla Monaghan (1992) instead.
* i.e. compile with: make CFLAGS=-DPLANETARY_SPH_NO_BALSARA
*/
/**
......@@ -495,7 +495,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo, const float dt_alpha) {
#ifdef PLANETARY_SPH_BALSARA
#ifndef PLANETARY_SPH_NO_BALSARA
const float fac_mu = cosmo->a_factor_mu;
/* Compute the norm of the curl */
......@@ -505,7 +505,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the norm of div v */
const float abs_div_v = fabsf(p->density.div_v);
#endif // PLANETARY_SPH_BALSARA
#endif
/* Compute the pressure */
const float pressure =
......@@ -527,20 +527,20 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
grad_h_term = 0.f;
}
#ifdef PLANETARY_SPH_BALSARA
#ifndef PLANETARY_SPH_NO_BALSARA
/* Compute the Balsara switch */
const float balsara =
abs_div_v / (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
#endif // PLANETARY_SPH_BALSARA
#endif
/* Update variables. */
p->force.f = grad_h_term;
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
#ifdef PLANETARY_SPH_BALSARA
#ifndef PLANETARY_SPH_NO_BALSARA
p->force.balsara = balsara;
#endif // PLANETARY_SPH_BALSARA
#endif
}
/**
......
......@@ -25,8 +25,8 @@
* @brief Minimal conservative implementation of SPH (Neighbour loop equations)
*
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
* viscosity term with the Balsara (1995) switch (optional).
* No thermal conduction term is implemented.
*
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with
* \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational
......@@ -176,11 +176,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
(pi->v[1] - pj->v[1]) * dx[1] +
(pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
#ifdef PLANETARY_SPH_BALSARA
#ifndef PLANETARY_SPH_NO_BALSARA
/* Balsara term */
const float balsara_i = pi->force.balsara;
const float balsara_j = pj->force.balsara;
#endif // PLANETARY_SPH_BALSARA
#endif
/* Are the particles moving towards each other? */
const float omega_ij = min(dvdr, 0.f);
......@@ -193,12 +193,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Now construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
#ifdef PLANETARY_SPH_BALSARA
#ifdef PLANETARY_SPH_NO_BALSARA
const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
#else
const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
(balsara_i + balsara_j) / rho_ij;
#else
const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
#endif // PLANETARY_SPH_BALSARA
#endif
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
......@@ -300,11 +300,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
(pi->v[1] - pj->v[1]) * dx[1] +
(pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
#ifdef PLANETARY_SPH_BALSARA
#ifndef PLANETARY_SPH_NO_BALSARA
/* Balsara term */
const float balsara_i = pi->force.balsara;
const float balsara_j = pj->force.balsara;
#endif // PLANETARY_SPH_BALSARA
#endif
/* Are the particles moving towards each other? */
const float omega_ij = min(dvdr, 0.f);
......@@ -319,12 +319,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
#ifdef PLANETARY_SPH_BALSARA
#ifdef PLANETARY_SPH_NO_BALSARA
const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
#else
const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
(balsara_i + balsara_j) / rho_ij;
#else
const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
#endif // PLANETARY_SPH_BALSARA
#endif
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
......
......@@ -25,8 +25,8 @@
* @brief Minimal conservative implementation of SPH (i/o routines)
*
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
* viscosity term with the Balsara (1995) switch (optional).
* No thermal conduction term is implemented.
*
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with
* \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
......@@ -197,14 +197,14 @@ INLINE static void hydro_write_flavour(hid_t h_grpsph) {
/* Viscosity and thermal conduction */
/* Nothing in this minimal model... */
io_write_attribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
#ifdef PLANETARY_SPH_BALSARA
#ifdef PLANETARY_SPH_NO_BALSARA
io_write_attribute_s(h_grpsph, "Viscosity Model",
"Minimal treatment as in Monaghan (1992)");
#else
io_write_attribute_s(
h_grpsph, "Viscosity Model",
"as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
#else
io_write_attribute_s(h_grpsph, "Viscosity Model",
"Minimal treatment as in Monaghan (1992)");
#endif // PLANETARY_SPH_BALSARA
#endif
/* Time integration properties */
io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
......
......@@ -25,8 +25,8 @@
* @brief Minimal conservative implementation of SPH (Particle definition)
*
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
* viscosity term with the Balsara (1995) switch (optional).
* No thermal conduction term is implemented.
*
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with
* \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational
......@@ -126,13 +126,13 @@ struct part {
/*! Derivative of density with respect to h */
float rho_dh;
#ifdef PLANETARY_SPH_BALSARA
#ifndef PLANETARY_SPH_NO_BALSARA
/*! Velocity divergence. */
float div_v;
/*! Velocity curl. */
float rot_v[3];
#endif // PLANETARY_SPH_BALSARA
#endif
} density;
......@@ -160,10 +160,10 @@ struct part {
/*! Time derivative of smoothing length */
float h_dt;
#ifdef PLANETARY_SPH_BALSARA
#ifndef PLANETARY_SPH_NO_BALSARA
/*! Balsara switch */
float balsara;
#endif // PLANETARY_SPH_BALSARA
#endif
} force;
};
......
......@@ -183,13 +183,14 @@ void hydro_props_print(const struct hydro_props *p) {
message("Minimal gas temperature set to %f", p->minimal_temperature);
// Matthieu: Temporary location for this i/o business.
#ifdef PLANETARY_SPH
#ifdef PLANETARY_SPH_BALSARA
message("Planetary SPH: Balsara switch enabled");
#ifdef PLANETARY_SPH_NO_BALSARA
message("Planetary SPH: Balsara switch DISABLED");
#else
message("Planetary SPH: Balsara switch disabled");
#endif // PLANETARY_SPH_BALSARA
#endif // PLANETARY_SPH
message("Planetary SPH: Balsara switch ENABLED");
#endif
#endif
}
#if defined(HAVE_HDF5)
......
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