Commit ac9208ba authored by Josh Borrow's avatar Josh Borrow Committed by Matthieu Schaller

Minor Updates to Hydro

parent 5fab1cdb
......@@ -1636,10 +1636,10 @@ esac
# Hydro scheme.
AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>],
[Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, sphenix, anarchy-pu default: gadget2@:>@]
[Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, phantom, gizmo-mfv, gizmo-mfm, shadowfax, planetary, sphenix, anarchy-pu default: sphenix@:>@]
)],
[with_hydro="$withval"],
[with_hydro="gadget2"]
[with_hydro="sphenix"]
)
case "$with_hydro" in
......@@ -1658,8 +1658,8 @@ case "$with_hydro" in
pressure-energy-monaghan)
AC_DEFINE([HOPKINS_PU_SPH_MONAGHAN], [1], [Pressure-Energy SPH with M&M Variable A.V.])
;;
default)
AC_DEFINE([DEFAULT_SPH], [1], [Default SPH])
phantom)
AC_DEFINE([PHANTOM_SPH], [1], [Phantom SPH])
;;
gizmo-mfv)
AC_DEFINE([GIZMO_MFV_SPH], [1], [GIZMO MFV SPH])
......
......@@ -19,6 +19,7 @@ schemes available in SWIFT, as well as how to implement your own.
hopkins_sph
anarchy_sph
sphenix_sph
phantom_sph
gizmo
adding_your_own
.. PHANTOM SPH
Josh Borrow 13th October 2020
Phantom
=======
This scheme is a reference implementation similar to the one presented in
Price (2018), the PHANTOM paper (not including MHD). It uses:
+ A simplified Cullen & Dehnen AV limiter (note that this is different to
PHANTOM as we do not explicitly include the matrix calculation).
+ A fixed alpha artificial conduction scheme used for hydro-only problems
as presented in the PHANTOM paper (i.e. we use the 'hydro-only' conduction
velocity, rather than the one used for gravitational problems).
+ Base Density-Energy SPH
The simplified version of the 'Inviscid SPH' artificial viscosity calculates
the time differential of the velocity divergence explicitly, using the value
from the previous step. We also use the Balsara switch instead of the improved
neighbour-based limiter from Cullen & Dehnen 2010, to avoid matrix
calculations. We also use a different value for the 'h-factors' due to SWIFT
using neighbour finding based on particle number density, rather than local
mass density.
To configure with this scheme, use
.. code-block:: bash
./configure --with-hydro=phantom --disable-hand-vec
The parameters available for this scheme, and their defaults, are:
.. code-block:: yaml
SPH:
viscosity_alpha: 0.1 # Initial value for the alpha viscosity
viscosity_length: 0.25 # Viscosity decay length (in terms of sound-crossing time)
# These are enforced each time-step
viscosity_alpha_max: 2.0 # Maximal allowed value for the viscosity alpha
viscosity_alpha_min: 0.0 # Minimal allowed value for the viscosity alpha
diffusion_alpha: 1.0 # Fixed value for the diffusion alpha
There is also a compile-time parameter, ``viscosity_beta`` that we set to
3.0. During feedback events, the viscosity is set to the compile-time
``hydro_props_default_viscosity_alpha_feedback_reset = 2.0`` and the
diffusion is set to ``hydro_props_default_diffusion_alpha_feedback_reset =
0.0``. These can be changed in ``src/hydro/Phantom/hydro_parameters.h``.
......@@ -37,9 +37,9 @@
#elif defined(HOPKINS_PU_SPH_MONAGHAN)
#error TODO
#include "./hydro/PressureEnergyMorrisMonaghanAV/logger_hydro.h"
#elif defined(DEFAULT_SPH)
#elif defined(PHANTOM_SPH)
#error TODO
#include "./hydro/Default/logger_hydro.h"
#include "./hydro/Phantom/logger_hydro.h"
#elif defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH)
#error TODO
#include "./hydro/Gizmo/logger_hydro.h"
......
......@@ -145,9 +145,9 @@ nobase_noinst_HEADERS += hydro.h hydro_io.h hydro_logger.h hydro_parameters.h
nobase_noinst_HEADERS += hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h
nobase_noinst_HEADERS += hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h
nobase_noinst_HEADERS += hydro/Minimal/hydro_parameters.h
nobase_noinst_HEADERS += hydro/Default/hydro.h hydro/Default/hydro_iact.h hydro/Default/hydro_io.h
nobase_noinst_HEADERS += hydro/Default/hydro_debug.h hydro/Default/hydro_part.h
nobase_noinst_HEADERS += hydro/Default/hydro_parameters.h
nobase_noinst_HEADERS += hydro/Phantom/hydro.h hydro/Phantom/hydro_iact.h hydro/Phantom/hydro_io.h
nobase_noinst_HEADERS += hydro/Phantom/hydro_debug.h hydro/Phantom/hydro_part.h
nobase_noinst_HEADERS += hydro/Phantom/hydro_parameters.h
nobase_noinst_HEADERS += hydro/Gadget2/hydro.h hydro/Gadget2/hydro_iact.h hydro/Gadget2/hydro_io.h
nobase_noinst_HEADERS += hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h
nobase_noinst_HEADERS += hydro/Gadget2/hydro_parameters.h hydro/Gadget2/hydro_logger.h
......
......@@ -50,8 +50,8 @@
#include "./hydro/PressureEnergy/hydro_debug.h"
#elif defined(HOPKINS_PU_SPH_MONAGHAN)
#include "./hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h"
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_debug.h"
#elif defined(PHANTOM_SPH)
#include "./hydro/Phantom/hydro_debug.h"
#elif defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH)
#include "./hydro/Gizmo/hydro_debug.h"
#elif defined(SHADOWFAX_SPH)
......
......@@ -51,10 +51,10 @@
#define SPH_IMPLEMENTATION \
"Pressure-Energy SPH (Hopkins 2013) with a Morris and Monaghan (1997) " \
"variable artificial viscosity."
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro.h"
#include "./hydro/Default/hydro_iact.h"
#define SPH_IMPLEMENTATION "Default version of SPH"
#elif defined(PHANTOM_SPH)
#include "./hydro/Phantom/hydro.h"
#include "./hydro/Phantom/hydro_iact.h"
#define SPH_IMPLEMENTATION "PHANTOM SPH reference implementation (Price 2018)"
#elif defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH)
#include "./hydro/Gizmo/hydro.h"
#include "./hydro/Gizmo/hydro_iact.h"
......
......@@ -598,17 +598,21 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
/* Compute the "grad h" term */
const float rho_inv = 1.f / p->rho;
float rho_dh = p->density.rho_dh;
/* Compute the "grad h" term - Note here that we have \tilde{x}
* as 1 as we use the local number density to find neighbours. This
* introduces a j-component that is considered in the force loop,
* meaning that this cached grad_h_term gives:
*
* f_ij = 1.f - grad_h_term_i / m_j */
const float common_factor = p->h * hydro_dimension_inv / p->density.wcount;
float grad_h_term = common_factor * p->density.rho_dh /
(1.f + common_factor * p->density.wcount_dh);
/* Ignore changing-kernel effects when h ~= h_max */
if (p->h > 0.9999f * hydro_props->h_max) {
rho_dh = 0.f;
grad_h_term = 0.f;
}
const float grad_h_term =
1.f / (1.f + hydro_dimension_inv * p->h * rho_dh * rho_inv);
/* Compute the Balsara switch */
/* Pre-multiply in the AV factor; hydro_props are not passed to the iact
* functions */
......
......@@ -235,9 +235,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
kernel_deval(xj, &wj, &wj_dx);
const float wj_dr = hjd_inv * wj_dx;
/* Variable smoothing length term */
const float f_ij = 1.f - pi->force.f / mj;
const float f_ji = 1.f - pj->force.f / mi;
/* Compute gradient terms */
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * f_ij;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * f_ji;
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
......@@ -265,7 +269,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float visc = -0.25f * v_sig * (balsara_i + balsara_j) * mu_ij / rho_ij;
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
const float visc_acc_term =
0.5f * visc * (wi_dr * f_ij + wj_dr * f_ji) * r_inv;
/* SPH acceleration term */
const float sph_acc_term =
......@@ -299,8 +304,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->u_dt += du_dt_j * mi;
/* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
pj->force.h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr * f_ij;
pj->force.h_dt -= mi * dvdr * r_inv / rhoi * wj_dr * f_ji;
/* Update the signal velocity. */
pi->force.v_sig = max(pi->force.v_sig, v_sig);
......@@ -339,7 +344,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float r = r2 * r_inv;
/* Recover some data */
// const float mi = pi->mass;
const float mi = pi->mass;
const float mj = pj->mass;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
......@@ -362,9 +367,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
kernel_deval(xj, &wj, &wj_dx);
const float wj_dr = hjd_inv * wj_dx;
/* Variable smoothing length term */
const float f_ij = 1.f - pi->force.f / mj;
const float f_ji = 1.f - pj->force.f / mi;
/* Compute gradient terms */
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * f_ij;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * f_ji;
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
......@@ -392,7 +401,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float visc = -0.25f * v_sig * (balsara_i + balsara_j) * mu_ij / rho_ij;
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
const float visc_acc_term =
0.5f * visc * (wi_dr * f_ij + wj_dr * f_ji) * r_inv;
/* SPH acceleration term */
const float sph_acc_term =
......@@ -419,7 +429,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->u_dt += du_dt_i * mj;
/* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr * f_ij;
/* Update the signal velocity. */
pi->force.v_sig = max(pi->force.v_sig, v_sig);
......
......@@ -17,15 +17,18 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_DEFAULT_HYDRO_H
#define SWIFT_DEFAULT_HYDRO_H
#ifndef SWIFT_PHANTOM_HYDRO_H
#define SWIFT_PHANTOM_HYDRO_H
/**
* @file Default/hydro.h
* @file Phantom/hydro.h
* @brief Density-Energy conservative implementation of SPH,
* with added diffusive physics (Cullen & Denhen 2011 AV,
* Price 2017 (PHANTOM) diffusion) (Non-neighbour loop
* equations)
* equations).
*
* This is a base reference implementation
* similar to the one presented in Price 2018.
*/
#include "adiabatic_index.h"
......@@ -589,17 +592,21 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
const float balsara =
abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed * fac_B / p->h);
/* Compute the "grad h" term */
const float rho_inv = 1.f / p->rho;
float rho_dh = p->density.rho_dh;
/* Compute the "grad h" term - Note here that we have \tilde{x}
* as 1 as we use the local number density to find neighbours. This
* introduces a j-component that is considered in the force loop,
* meaning that this cached grad_h_term gives:
*
* f_ij = 1.f - grad_h_term_i / m_j */
const float common_factor = p->h * hydro_dimension_inv / p->density.wcount;
float grad_h_term = common_factor * p->density.rho_dh /
(1.f + common_factor * p->density.wcount_dh);
/* Ignore changing-kernel effects when h ~= h_max */
if (p->h > 0.9999f * hydro_props->h_max) {
rho_dh = 0.f;
grad_h_term = 0.f;
}
const float grad_h_term =
1.f / (1.f + hydro_dimension_inv * p->h * rho_dh * rho_inv);
/* Update variables. */
p->force.f = grad_h_term;
p->force.pressure = pressure;
......@@ -1021,4 +1028,4 @@ hydro_set_init_internal_energy(struct part *p, float u_init) {
p->u = u_init;
}
#endif /* SWIFT_DEFAULT_HYDRO_H */
#endif /* SWIFT_PHANTOM_HYDRO_H */
......@@ -18,14 +18,17 @@
*
******************************************************************************/
#ifndef SWIFT_DEFAULT_HYDRO_DEBUG_H
#define SWIFT_DEFAULT_HYDRO_DEBUG_H
#ifndef SWIFT_PHANTOM_HYDRO_DEBUG_H
#define SWIFT_PHANTOM_HYDRO_DEBUG_H
/**
* @file Default/hydro_debug.h
* @file Phantom/hydro_debug.h
* @brief Density-Energy conservative implementation of SPH,
* with added diffusive physics (Cullen & Denhen 2011 AV,
* Price 2017 (PHANTOM) diffusion) (Debugging routines)
*
* This is a base reference implementation
* similar to the one presented in Price 2018.
*/
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
......@@ -44,4 +47,4 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
p->viscosity.alpha, p->time_bin);
}
#endif /* SWIFT_DEFAULT_HYDRO_DEBUG_H */
#endif /* SWIFT_PHANTOM_HYDRO_DEBUG_H */
......@@ -17,14 +17,17 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_DEFAULT_HYDRO_IACT_H
#define SWIFT_DEFAULT_HYDRO_IACT_H
#ifndef SWIFT_PHANTOM_HYDRO_IACT_H
#define SWIFT_PHANTOM_HYDRO_IACT_H
/**
* @file Default/hydro_iact.h
* @file Phantom/hydro_iact.h
* @brief Density-Energy conservative implementation of SPH,
* with added diffusive physics (Cullen & Denhen 2011 AV,
* Price 2017 (PHANTOM) diffusion) (interaction routines)
*
* This is a base reference implementation
* similar to the one presented in Price 2018.
*/
#include "adiabatic_index.h"
......@@ -330,6 +333,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float v_sig = pi->force.soundspeed + pj->force.soundspeed -
const_viscosity_beta * mu_ij;
/* Variable smoothing length term */
const float f_ij = 1.f - pi->force.f / mj;
const float f_ji = 1.f - pj->force.f / mi;
/* Balsara term */
const float balsara_i = pi->force.balsara;
const float balsara_j = pj->force.balsara;
......@@ -340,12 +347,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Convolve with the kernel */
const float visc_acc_term =
0.5f * visc * (wi_dr * pi->force.f / rhoi + wj_dr * pj->force.f / rhoj) *
r_inv;
0.5f * visc * (wi_dr * f_ij / rhoi + wj_dr * f_ji / rhoj) * r_inv;
/* Compute gradient terms */
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * f_ij;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * f_ji;
/* SPH acceleration term */
const float sph_acc_term =
......@@ -378,9 +384,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
/* wi_dx + wj_dx / 2 is F_ij */
const float diff_du_term =
alpha_diff * v_diff * (pi->u - pj->u) * 0.5f *
(wi_dr * pi->force.f / pi->rho + wj_dr * pi->force.f / pi->rho);
const float diff_du_term = alpha_diff * v_diff * (pi->u - pj->u) * 0.5f *
(wi_dr * f_ij / pi->rho + wj_dr * f_ji / pi->rho);
/* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;
......@@ -391,8 +396,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->u_dt += du_dt_j * mi;
/* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * pi->force.f * r_inv / rhoj * wi_dr;
pj->force.h_dt -= mi * dvdr * pj->force.f * r_inv / rhoi * wj_dr;
pi->force.h_dt -= mj * dvdr * f_ij * r_inv / rhoj * wi_dr;
pj->force.h_dt -= mi * dvdr * f_ji * r_inv / rhoi * wj_dr;
}
/**
......@@ -419,7 +424,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float r_inv = 1.0f / r;
/* Recover some data */
// const float mi = pi->mass;
const float mi = pi->mass;
const float mj = pj->mass;
const float rhoi = pi->rho;
......@@ -460,6 +465,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float v_sig = pi->force.soundspeed + pj->force.soundspeed -
const_viscosity_beta * mu_ij;
/* Variable smoothing length term */
const float f_ij = 1.f - pi->force.f / mj;
const float f_ji = 1.f - pj->force.f / mi;
/* Balsara term */
const float balsara_i = pi->force.balsara;
const float balsara_j = pj->force.balsara;
......@@ -470,12 +479,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Convolve with the kernel */
const float visc_acc_term =
0.5f * visc * (wi_dr * pi->force.f / rhoi + wj_dr * pj->force.f / rhoj) *
r_inv;
0.5f * visc * (wi_dr * f_ij / rhoi + wj_dr * f_ji / rhoj) * r_inv;
/* Compute gradient terms */
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * f_ij;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * f_ji;
/* SPH acceleration term */
const float sph_acc_term =
......@@ -502,9 +510,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
/* wi_dx + wj_dx / 2 is F_ij */
const float diff_du_term =
alpha_diff * v_diff * (pi->u - pj->u) * 0.5f *
(wi_dr * pi->force.f / pi->rho + wj_dr * pi->force.f / pi->rho);
const float diff_du_term = alpha_diff * v_diff * (pi->u - pj->u) * 0.5f *
(wi_dr * f_ij / pi->rho + wj_dr * f_ji / pi->rho);
/* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;
......@@ -513,7 +520,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->u_dt += du_dt_i * mj;
/* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * pi->force.f * r_inv / rhoj * wi_dr;
pi->force.h_dt -= mj * dvdr * f_ij * r_inv / rhoj * wi_dr;
}
#endif /* SWIFT_DEFAULT_HYDRO_IACT_H */
#endif /* SWIFT_PHANTOM_HYDRO_IACT_H */
......@@ -17,14 +17,16 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_DEFAULT_HYDRO_IO_H
#define SWIFT_DEFAULT_HYDRO_IO_H
#ifndef SWIFT_PHANTOM_HYDRO_IO_H
#define SWIFT_PHANTOM_HYDRO_IO_H
/**
* @file Default/hydro_io.h
* @file Phantom/hydro_io.h
* @brief Density-Energy conservative implementation of SPH,
* with added diffusive physics (Cullen & Denhen 2011 AV,
* Price 2017 (PHANTOM) diffusion) (i/o routines)
* This is a base reference implementation
* similar to the one presented in Price 2018.
*/
#include "adiabatic_index.h"
......@@ -237,4 +239,4 @@ INLINE static void hydro_write_flavour(hid_t h_grpsph) {
*/
INLINE static int writeEntropyFlag(void) { return 0; }
#endif /* SWIFT_DEFAULT_HYDRO_IO_H */
#endif /* SWIFT_PHANTOM_HYDRO_IO_H */
......@@ -18,8 +18,8 @@
*
******************************************************************************/
#ifndef SWIFT_DEFAULT_HYDRO_PARAMETERS_H
#define SWIFT_DEFAULT_HYDRO_PARAMETERS_H
#ifndef SWIFT_PHANTOM_HYDRO_PARAMETERS_H
#define SWIFT_PHANTOM_HYDRO_PARAMETERS_H
/* Configuration file */
#include "config.h"
......@@ -35,12 +35,15 @@
#include "inline.h"
/**
* @file Default/hydro_parameters.h
* @file Phantom/hydro_parameters.h
* @brief Density-Energy conservative implementation of SPH,
* with added diffusive physics (Cullen & Denhen 2011 AV,
* Price 2017 (PHANTOM) diffusion) (default compile-time
* parameters).
*
* This is a base reference implementation
* similar to the one presented in Price 2018.
*
* This file defines a number of things that are used in
* hydro_properties.c as defaults for run-time parameters
* as well as a number of compile-time parameters.
......@@ -253,4 +256,4 @@ static INLINE void diffusion_print_snapshot(
}
#endif
#endif /* SWIFT_DEFAULT_HYDRO_PARAMETERS_H */
#endif /* SWIFT_PHANTOM_HYDRO_PARAMETERS_H */
......@@ -17,14 +17,17 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_DEFAULT_HYDRO_PART_H
#define SWIFT_DEFAULT_HYDRO_PART_H
#ifndef SWIFT_PHANTOM_HYDRO_PART_H
#define SWIFT_PHANTOM_HYDRO_PART_H
/**
* @file Default/hydro_part.h
* @file Phantom/hydro_part.h
* @brief Density-Energy conservative implementation of SPH,
* with added diffusive physics (Cullen & Denhen 2011 AV,
* Price 2017 (PHANTOM) diffusion) (particle definition)
*
* This is a base reference implementation
* similar to the one presented in Price 2018.
*/
#include "black_holes_struct.h"
......@@ -220,4 +223,4 @@ struct part {
} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_DEFAULT_HYDRO_PART_H */
#endif /* SWIFT_PHANTOM_HYDRO_PART_H */
......@@ -34,8 +34,8 @@
#include "./hydro/PressureEnergy/hydro_io.h"
#elif defined(HOPKINS_PU_SPH_MONAGHAN)
#include "./hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h"
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_io.h"
#elif defined(PHANTOM_SPH)
#include "./hydro/Phantom/hydro_io.h"
#elif defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH)
#include "./hydro/Gizmo/hydro_io.h"
#elif defined(SHADOWFAX_SPH)
......
......@@ -39,7 +39,7 @@
#error TODO
#elif defined(HOPKINS_PU_SPH_MONAGHAN)
#error TODO
#elif defined(DEFAULT_SPH)
#elif defined(PHANTOM_SPH)
#error TODO
#elif defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH)
#error TODO
......
......@@ -42,8 +42,8 @@
#include "./hydro/PressureEnergy/hydro_parameters.h"
#elif defined(HOPKINS_PU_SPH_MONAGHAN)
#include "./hydro/PressureEnergyMorrisMonaghanAV/hydro_parameters.h"
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_parameters.h"
#elif defined(PHANTOM_SPH)
#include "./hydro/Phantom/hydro_parameters.h"
#elif defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH)
#include "./hydro/Gizmo/hydro_parameters.h"
#elif defined(SHADOWFAX_SPH)
......
......@@ -61,8 +61,8 @@ struct threadpool;
#elif defined(HOPKINS_PU_SPH_MONAGHAN)
#include "./hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h"
#define hydro_need_extra_init_loop 0
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_part.h"
#elif defined(PHANTOM_SPH)
#include "./hydro/Phantom/hydro_part.h"
#define EXTRA_HYDRO_LOOP
#define hydro_need_extra_init_loop 0
#elif defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH)
......
......@@ -109,11 +109,11 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
part->entropy = pressure / pow_gamma(density);
#elif defined(HOPKINS_PE_SPH)
part->entropy = pressure / pow_gamma(density);
#elif defined(DEFAULT_SPH)
#elif defined(PHANTOM_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH) || \
defined(SPHENIX_SPH) || defined(DEFAULT_SPH)
defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(PLANETARY_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
......@@ -405,7 +405,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH) || \
defined(HOPKINS_PU_SPH) || defined(HOPKINS_PU_SPH_MONAGHAN)
0.f,
#elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(DEFAULT_SPH)
#elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
main_cell->hydro.parts[pid].viscosity.div_v,
#else
main_cell->hydro.parts[pid].density.div_v,
......@@ -426,7 +426,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
defined(HOPKINS_PU_SPH_MONAGHAN)
main_cell->hydro.parts[pid].force.v_sig, 0.f,
main_cell->hydro.parts[pid].u_dt
#elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(DEFAULT_SPH)
#elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
main_cell->hydro.parts[pid].viscosity.v_sig, 0.f,
main_cell->hydro.parts[pid].u_dt
#else
......
......@@ -288,7 +288,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->hydro.parts[pid].density.rot_v[0],
main_cell->hydro.parts[pid].density.rot_v[1],
main_cell->hydro.parts[pid].density.rot_v[2]
#elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(DEFAULT_SPH)
#elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
/* this is required because of the variable AV scheme */
main_cell->hydro.parts[pid].viscosity.div_v,
main_cell->hydro.parts[pid].density.rot_v[0],
......@@ -333,7 +333,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
cj->hydro.parts[pjd].density.rot_v[0],
cj->hydro.parts[pjd].density.rot_v[1],
cj->hydro.parts[pjd].density.rot_v[2]
#elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(DEFAULT_SPH)
#elif defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
/* this is required because of the variable AV scheme */
cj->hydro.parts[pjd].viscosity.div_v,
cj->hydro.parts[pjd].density.rot_v[0],
......
......@@ -113,7 +113,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
part->entropy = 1.f;
#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH) || \
defined(SPHENIX_SPH) || defined(DEFAULT_SPH)
defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
part->u = 1.f;
#elif defined(HOPKINS_PE_SPH)
part->entropy = 1.f;
......@@ -197,11 +197,12 @@ void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo,
p->density.rot_v[2] = 0.f;
p->density.div_v = 0.f;
#endif /* GADGET-2 */
#if defined(MINIMAL_SPH) || defined(SPHENIX_SPH) || defined(DEFAULT_SPH)
#if defined(MINIMAL_SPH) || defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
p->rho = 1.f;
p->density.rho_dh = 0.f;
p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h));
p->density.wcount_dh = 0.f;
p->viscosity.v_sig = hydro_get_comoving_soundspeed(p);
#endif /* MINIMAL */
#ifdef HOPKINS_PE_SPH
p->rho = 1.f;
......
......@@ -114,7 +114,7 @@ void prepare_force(struct part *parts, size_t count) {
#if !defined(GIZMO_MFV_SPH) && !defined(SHADOWFAX_SPH) && \
!defined(MINIMAL_SPH) && !defined(PLANETARY_SPH) && \
!defined(HOPKINS_PU_SPH) && !defined(HOPKINS_PU_SPH_MONAGHAN) && \
!defined(ANARCHY_PU_SPH) && !defined(SPHENIX_SPH) && !defined(DEFAULT_SPH)
!defined(ANARCHY_PU_SPH) && !defined(SPHENIX_SPH) && !defined(PHANTOM_SPH)
struct part *p;
for (size_t i = 0; i < count; ++i) {
p = &parts[i];
......@@ -142,7 +142,7 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->h,
hydro_get_comoving_density(p),
#if defined(MINIMAL_SPH) || defined(PLANETARY_SPH) || \
defined(SHADOWFAX_SPH) || defined(DEFAULT_SPH)
defined(SHADOWFAX_SPH) || defined(PHANTOM_SPH)
0.f,
#else
p->density.div_v,
......@@ -153,11 +153,11 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
p->a_hydro[0], p->a_hydro[1], p->a_hydro[2], p->force.h_dt,
#if defined(GADGET2_SPH)
p->force.v_sig, p->entropy_dt, 0.f
#elif defined(DEFAULT_SPH)
#elif defined(PHANTOM_SPH)
p->force.v_sig, 0.f, p->force.u_dt
#elif defined(MINIMAL_SPH) || defined(HOPKINS_PU_SPH) || \
defined(HOPKINS_PU_SPH_MONAGHAN) || defined(ANARCHY_PU_SPH) || \
defined(SPHENIX_SPH) || defined(DEFAULT_SPH)
defined(SPHENIX_SPH) || defined(PHANTOM_SPH)
p->force.v_sig, 0.f, p->u_dt
#else
0.f, 0.f, 0.f
......
......@@ -251,7 +251,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell, int i, int j,
main_cell->hydro.parts[pid].density.rot_v[0],
main_cell->hydro.parts[pid].density.rot_v[1],
main_cell->hydro.parts[pid].density.rot_v[2]
#elif defined(DEFAULT_SPH) || defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH)
#elif defined(PHANTOM_SPH) || defined(ANARCHY_PU_SPH) || defined(SPHENIX_SPH)
main_cell->hydro.parts[pid].viscosity.div_v,
main_cell->hydro.parts[pid].density.rot_v[0],
main_cell->hydro.parts[pid].density.rot_v[1],
......
Markdown is supported
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