Skip to content
Snippets Groups Projects
Commit 0596dfcf authored by Josh Borrow's avatar Josh Borrow
Browse files

Finally got a version of P-U working...

parent a7c5d52d
No related branches found
No related tags found
1 merge request!540Add Pressure-Energy (P-U) SPH
Hide whitespace changes
Inline Side-by-side
src/hydro/PressureEnergy/hydro.h
+ 177
− 172
View file @ 0596dfcf
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
......@@ -20,15 +20,17 @@
#define SWIFT_PRESSURE_ENERGY_HYDRO_H
/**
* @file PressureEnergy/hydro.h
* @brief Pressure-Energy implementation of SPH (Non-neighbour loop
* @file Minimal/hydro.h
* @brief Minimal conservative implementation of SPH (Non-neighbour loop
* equations)
*
* The thermal variable is the energy (u) and the pressure is smoothed over
* contact discontinuities to prevent spurious surface tension.
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
*
* Follows equations (16), (17) and (18) of Hopkins, P., MNRAS, 2013,
* Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
* 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
* Physics, 2012, Volume 231, Issue 3, pp. 759-794.
*/
#include "adiabatic_index.h"
......@@ -37,6 +39,7 @@
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "hydro_space.h"
#include "kernel_hydro.h"
#include "minmax.h"
......@@ -45,10 +48,14 @@
/**
* @brief Returns the comoving internal energy of a particle
*
* For implementations where the main thermodynamic variable
* is not internal energy, this function computes the internal
* energy from the thermodynamic variable.
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_comoving_internal_energy(
const struct part *restrict p) {
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_internal_energy(const struct part *restrict p) {
return p->u;
}
......@@ -56,11 +63,17 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_internal_e
/**
* @brief Returns the physical internal energy of a particle
*
* @param p The particle of interest
* For implementations where the main thermodynamic variable
* is not internal energy, this function computes the internal
* energy from the thermodynamic variable and converts it to
* physical coordinates.
*
* @param p The particle of interest.
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float hydro_get_physical_internal_energy(
const struct part *restrict p,
const struct cosmology *cosmo) {
__attribute__((always_inline)) INLINE static float
hydro_get_physical_internal_energy(const struct part *restrict p,
const struct cosmology *cosmo) {
return p->u * cosmo->a_factor_internal_energy;
}
......@@ -68,6 +81,8 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_internal_e
/**
* @brief Returns the comoving pressure of a particle
*
* Computes the pressure based on the particle's properties.
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
......@@ -79,20 +94,25 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
/**
* @brief Returns the physical pressure of a particle
*
* Computes the pressure based on the particle's properties and
* convert it to physical coordinates.
*
* @param p The particle of interest
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
const struct part *restrict p,
const struct cosmology *cosmo) {
const struct part *restrict p, const struct cosmology *cosmo) {
return gas_pressure_from_internal_energy(p->rho, p->u) * cosmo->a_factor_pressure;
return cosmo->a_factor_pressure *
gas_pressure_from_internal_energy(p->rho, p->u);
}
/**
* @brief Returns the comoving entropy of a particle
*
* This function actually calculates the entropy from the internal energy of the
* particle, as in PressureEnergy (unsurprisingly) we follow U.
* For implementations where the main thermodynamic variable
* is not entropy, this function computes the entropy from
* the thermodynamic variable.
*
* @param p The particle of interest
*/
......@@ -105,18 +125,19 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
/**
* @brief Returns the physical entropy of a particle
*
* This function actually calculates the entropy from the internal energy of the
* particle, as in PressureEnergy (unsurprisingly) we follow U.
* For implementations where the main thermodynamic variable
* is not entropy, this function computes the entropy from
* the thermodynamic variable and converts it to
* physical coordinates.
*
* @param p The particle of interest
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
const struct part *restrict p,
const struct cosmology *cosmo) {
/* Note we choose co-ordinates such that this always works out to
* require no conversion */
const struct part *restrict p, const struct cosmology *cosmo) {
/* Note: no cosmological conversion required here with our choice of
* coordinates. */
return gas_entropy_from_internal_energy(p->rho, p->u);
}
......@@ -125,8 +146,8 @@ __attribute__((always_inline)) INLINE static float hydro_get_physical_entropy(
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_comoving_soundspeed(
const struct part *restrict p) {
__attribute__((always_inline)) INLINE static float
hydro_get_comoving_soundspeed(const struct part *restrict p) {
return p->force.soundspeed;
}
......@@ -135,12 +156,13 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_soundspeed
* @brief Returns the physical sound speed of a particle
*
* @param p The particle of interest
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float hydro_get_physical_soundspeed(
const struct part *restrict p,
const struct cosmology *cosmo) {
__attribute__((always_inline)) INLINE static float
hydro_get_physical_soundspeed(const struct part *restrict p,
const struct cosmology *cosmo) {
return p->force.soundspeed * cosmo->a_factor_sound_speed;
return cosmo->a_factor_sound_speed * p->force.soundspeed;
}
/**
......@@ -155,14 +177,15 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
}
/**
* @brief Returns the physical density of a particle
* @brief Returns the comoving density of a particle.
*
* @param p The particle of interest
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float hydro_get_physical_density(
const struct part *restrict p, const struct cosmology *restrict cosmo) {
const struct part *restrict p, const struct cosmology *cosmo) {
return p->rho * cosmo->a3_inv;
return cosmo->a3_inv * p->rho;
}
/**
......@@ -176,13 +199,13 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
return p->mass;
}
/**
* @brief Returns the velocities drifted to the current time of a particle.
*
* @param p The particle of interest
* @param xp The extended data of the particle.
* @param dt The time since the last kick.
* @param dt_kick_hydro The time (for hydro accelerations) since the last kick.
* @param dt_kick_grav The time (for gravity accelerations) since the last kick.
* @param v (return) The velocities at the current time.
*/
__attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
......@@ -190,11 +213,11 @@ __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
float dt_kick_grav, float v[3]) {
v[0] = xp->v_full[0] + p->a_hydro[0] * dt_kick_hydro +
xp->a_grav[0] * dt_kick_grav;
xp->a_grav[0] * dt_kick_grav;
v[1] = xp->v_full[1] + p->a_hydro[1] * dt_kick_hydro +
xp->a_grav[1] * dt_kick_grav;
xp->a_grav[1] * dt_kick_grav;
v[2] = xp->v_full[2] + p->a_hydro[2] * dt_kick_hydro +
xp->a_grav[2] * dt_kick_grav;
xp->a_grav[2] * dt_kick_grav;
}
/**
......@@ -204,52 +227,48 @@ __attribute__((always_inline)) INLINE static void hydro_get_drifted_velocities(
*
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static float hydro_get_energy_dt(
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
const struct part *restrict p) {
/* Since we track energy, we can forget about entropy conversions etc. */
return p->u_dt;
}
/**
* @brief Sets the time derivative of internal energy of a particle
* @brief Returns the time derivative of internal energy of a particle
*
* We assume a constant density.
*
* @param p The particle of interest.
* @param du_dt The new time derivative of the internal energy.
*/
__attribute__((always_inline)) INLINE static void hydro_set_energy_dt(
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
struct part *restrict p, float du_dt) {
/* Since we track energy, we can forget about entropy conversions etc. */
p->u_dt = du_dt;
}
/**
* @brief Computes the hydro time-step of a given particle
*
* This function returns the time-step of a particle given its hydro-dynamical
* state. A typical time-step calculation would be the use of the CFL condition.
*
* @param p Pointer to the particle data
* @param xp Pointer to the extended particle data
*
* @param hydro_properties The SPH parameters
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const struct part *restrict p, const struct xpart *restrict xp,
const struct hydro_props *restrict hydro_properties,
const struct cosmology *restrict cosmo) {
const float CFL_condition = hydro_properties->CFL_condition;
/* CFL condition */
const float dt_cfl =
2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
(p->force.v_sig * cosmo->a_factor_sound_speed);
const float dt_cfl = 2.f * kernel_gamma * CFL_condition * cosmo->a * p->h /
(cosmo->a_factor_sound_speed * p->force.v_sig);
/* U change limited */
/* TODO: Ensure this is consistent with cosmology */
const float dt_u_changes =
(p->u_dt != 0.0f) ? fabsf(const_max_u_change * p->u / p->u_dt)
: FLT_MAX;
return min(dt_u_changes, dt_cfl);
return dt_cfl;
}
/**
......@@ -266,23 +285,21 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
* @brief Prepares a particle for the density calculation.
*
* Zeroes all the relevant arrays in preparation for the sums taking place in
* the variaous density tasks
* the various density loop over neighbours. Typically, all fields of the
* density sub-structure of a particle get zeroed in here.
*
* @param p The particle to act upon
* @param hs #hydro_space containing hydro specific space information.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part *restrict p, const struct hydro_space *hs) {
p->rho = 0.f;
p->pressure_bar = 0.f;
p->pressure_bar_dh = 0.f;
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->density.div_v = 0.f;
p->density.rot_v[0] = 0.f;
p->density.rot_v[1] = 0.f;
p->density.rot_v[2] = 0.f;
p->rho = 0.f;
p->density.rho_dh = 0.f;
p->pressure_bar = 0.f;
p->density.pressure_bar_dh = 0.f;
}
/**
......@@ -290,8 +307,13 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
*
* Multiplies the density and number of neighbours by the appropiate constants
* and add the self-contribution term.
* Additional quantities such as velocity gradients will also get the final
* terms added to them here.
*
* Also adds/multiplies the cosmological terms if need be.
*
* @param p The particle to act upon
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
struct part *restrict p, const struct cosmology *cosmo) {
......@@ -301,40 +323,34 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
const float h_inv = 1.0f / h; /* 1/h */
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
/* Final operation on the density (add self-contribution). */
p->density.wcount += kernel_root;
p->density.wcount_dh -= hydro_dimension * kernel_root * h_inv;
p->rho += p->mass * kernel_root;
p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
p->pressure_bar += p->mass * p->u * kernel_root;
p->pressure_bar_dh -=
p->mass * p->u * h_inv * (hydro_dimension * kernel_root);
p->density.pressure_bar_dh -= hydro_dimension * p->mass * p->u * kernel_root;
p->density.wcount += kernel_root;
p->density.wcount_dh -= hydro_dimension * kernel_root;
/* Finish the calculation by inserting the missing h-factors */
p->rho *= h_inv_dim;
p->density.rho_dh *= h_inv_dim_plus_one;
p->pressure_bar *= (h_inv_dim * hydro_gamma_minus_one);
p->density.pressure_bar_dh *= (h_inv_dim_plus_one * hydro_gamma_minus_one);
p->density.wcount *= h_inv_dim;
p->density.wcount_dh *= h_inv_dim;
p->pressure_bar *= h_inv_dim;
p->pressure_bar_dh *= h_inv_dim;
/* Insert gamma factors */
p->pressure_bar *= hydro_gamma_minus_one;
p->pressure_bar_dh *= hydro_gamma_minus_one;
/* Finish calculation of the velocity curl components */
const float factor = h_inv_dim_plus_one / p->rho;
p->density.rot_v[0] *= factor;
p->density.rot_v[1] *= factor;
p->density.rot_v[2] *= factor;
/* Finish calculation of the velocity divergence */
p->density.div_v *= factor;
p->density.wcount_dh *= h_inv_dim_plus_one;
}
/**
* @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
*
* In the desperate case where a particle has no neighbours (likely because
* of the h_max ceiling), set the particle fields to something sensible to avoid
* NaNs in the next calculations.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
struct part *restrict p, struct xpart *restrict xp,
......@@ -346,68 +362,54 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
/* Re-set problematic values */
p->rho = p->mass * kernel_root * h_inv_dim;
p->pressure_bar = p->mass * p->u * hydro_gamma_minus_one * kernel_root * h_inv_dim;
p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
p->density.rho_dh = 0.f;
p->density.wcount_dh = 0.f;
p->pressure_bar =
kernel_root * p->mass * p->u * hydro_gamma_minus_one * h_inv_dim;
p->rho = kernel_root * p->mass * h_inv_dim;
p->pressure_bar_dh = 0;
p->density.div_v = 0.f;
p->density.rot_v[0] = 0.f;
p->density.rot_v[1] = 0.f;
p->density.rot_v[2] = 0.f;
p->density.pressure_bar_dh = 0.f;
}
/**
* @brief Prepare a particle for the force calculation.
*
* Computes viscosity term, conduction term and smoothing length gradient terms.
* This function is called in the ghost task to convert some quantities coming
* from the density loop over neighbours into quantities ready to be used in the
* force loop over neighbours. Quantities are typically read from the density
* sub-structure and written to the force sub-structure.
* Examples of calculations done here include the calculation of viscosity term
* constants, thermal conduction terms, hydro conversions, etc.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) {
const float fac_mu = 1.f; /* Will change with cosmological integration */
/* Compute the norm of the curl */
const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
p->density.rot_v[1] * p->density.rot_v[1] +
p->density.rot_v[2] * p->density.rot_v[2]);
/* Compute the norm of div v */
const float abs_div_v = fabsf(p->density.div_v);
/* Compute the pressure */
const float pressure =
gas_pressure_from_internal_energy(p->rho, p->u);
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
/* Compute the sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
/* Compute this `useful' property */
/* Compute the "grad h" term */
const float rho_inv = 1.f / p->rho;
const float P_over_rho2 = pressure * rho_inv * rho_inv;
/* Compute the Balsara switch */
const float balsara =
abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);
const float grad_h_term =
1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
/* Update variables */
/* Update variables. */
p->force.f = grad_h_term;
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
p->force.balsara = balsara;
p->force.P_over_rho2 = P_over_rho2;
p->force.f = 1.f; /* This cannot be computed individually in Pressure-Energy */
}
/**
* @brief Reset acceleration fields of a particle
*
* Resets all hydro acceleration and time derivative fields in preparation
* for the sums taking place in the variaous force tasks
* for the sums taking place in the various force tasks.
*
* @param p The particle to act upon
*/
......@@ -420,10 +422,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
p->a_hydro[2] = 0.0f;
/* Reset the time derivatives. */
p->force.h_dt = 0.0f;
p->u_dt = 0.0f;
/* Reset maximal signal velocity */
p->force.h_dt = 0.0f;
p->force.v_sig = 0.0f;
}
......@@ -442,15 +442,23 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
p->v[1] = xp->v_full[1];
p->v[2] = xp->v_full[2];
/* Re-set the entropy */
p->u = xp->u_full;
}
/**
* @brief Predict additional particle fields forward in time when drifting
*
* @param p The particle
* @param xp The extended data of the particle
* @param dt The drift time-step.
* Additional hydrodynamic quantites are drifted forward in time here. These
* include thermal quantities (thermal energy or total energy or entropy, ...).
*
* Note the different time-step sizes used for the different quantities as they
* include cosmological factors.
*
* @param p The particle.
* @param xp The extended data of the particle.
* @param dt_drift The drift time-step for positions.
* @param dt_therm The drift time-step for thermal quantities.
*/
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part *restrict p, const struct xpart *restrict xp, float dt_drift,
......@@ -466,99 +474,102 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->h *= expf(w1);
/* Predict density */
const float w2 = -hydro_dimension * w1;
if (fabsf(w2) < 0.2f) {
if (fabsf(w2) < 0.2f)
p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
p->pressure_bar *= approx_expf(w2);
} else {
else
p->rho *= expf(w2);
p->pressure_bar *= expf(w2);
}
/* Predict the energy */
/* Predict the internal energy */
p->u += p->u_dt * dt_therm;
/* Compute the pressure */
/* Compute the new pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
/* Compute the new sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
}
/**
* @brief Finishes the force calculation.
*
* Multiplies the forces and accelerationsby the appropiate constants
* Multiplies the force and accelerations by the appropiate constants
* and add the self-contribution term. In most cases, there is little
* to do here.
*
* Cosmological terms are also added/multiplied here.
*
* @param p The particle to act upon
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part *restrict p, const struct cosmology *cosmo) {
p->force.h_dt *= p->h * hydro_dimension_inv;
}
/**
* @brief Kick the additional variables
*
* @param p The particle to act upon
* @param xp The particle extended data to act upon
* @param dt The time-step for this kick
* @param half_dt The half time-step for this kick
* Additional hydrodynamic quantites are kicked forward in time here. These
* include thermal quantities (thermal energy or total energy or entropy, ...).
*
* @param p The particle to act upon.
* @param xp The particle extended data to act upon.
* @param dt_therm The time-step for this kick (for thermodynamic quantities).
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt_therm) {
/* Do not decrease the energy (temperature) by more than a factor of 2*/
if (dt_therm > 0. && p->u_dt * dt_therm < -0.5f * p->u) {
p->u_dt = -0.5f * p->u / dt_therm;
/* Do not decrease the energy by more than a factor of 2*/
if (dt_therm > 0. && p->u_dt * dt_therm < -0.5f * xp->u_full) {
p->u_dt = -0.5f * xp->u_full / dt_therm;
}
xp->u_full += p->u_dt * dt_therm;
/* Compute the pressure */
const float pressure =
gas_pressure_from_internal_energy(p->rho, xp->u_full);
const float pressure = gas_pressure_from_internal_energy(p->rho, xp->u_full);
/* Compute the new sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
/* Compute the sound speed */
const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
/* Update the sound speed */
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
}
/**
* @brief Converts hydro quantity of a particle at the start of a run
* @brief Converts hydro quantity of a particle at the start of a run
*
* Requires the density to be known
* This function is called once at the end of the engine_init_particle()
* routine (at the start of a calculation) after the densities of
* particles have been computed.
* This can be used to convert internal energy into entropy for instance.
*
* @param p The particle to act upon
* @param xp The extended particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p, struct xpart *restrict xp) {
p->entropy = hydro_get_comoving_entropy(p);
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
/* Compute the sound speed */
const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
const float rho_inv = 1.f / p->rho;
const float P_over_rho2 = pressure * rho_inv * rho_inv;
/* Compute the sound speed */
const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
/* Update sound speed */
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = P_over_rho2;
}
/**
* @brief Initialises the particles for the first time
*
* This function is called only once just after the ICs have been
* read in to do some conversions.
* read in to do some conversions or assignments between the particle
* and extended particle fields.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
......@@ -567,25 +578,18 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part *restrict p, struct xpart *restrict xp) {
p->time_bin = 0;
p->rho = 0.f;
p->pressure_bar = 0.f;
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
xp->a_grav[0] = 0.f;
xp->a_grav[1] = 0.f;
xp->a_grav[2] = 0.f;
xp->u_full = p->u;
hydro_reset_acceleration(p);
hydro_init_part(p, NULL);
}
/* Begin actual helper functions for the specific hydro implementation.
*
* Here we define:
*
* hydro_h_term(*pi, *pj)
*/
/**
* @brief Calculates the f_ij factors for a given particle i and j (i.e. the
* 'h-terms'.
......@@ -602,7 +606,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
* @param hi The smoothing length of particle i
*/
__attribute__((always_inline)) INLINE static float hydro_h_term(
struct part *restrict pi, struct part *restrict pj, float hi) {
const struct part * pi, const struct part * pj, float hi) {
/* First constrct h_i / (n_D n_i) */
......@@ -612,7 +616,7 @@ __attribute__((always_inline)) INLINE static float hydro_h_term(
const float bottom_term = 1 + hi_over_density * pi->density.wcount_dh;
/* Construct the 'top term' */
const float top_term_top = pi->pressure_bar_dh * hi_over_density;
const float top_term_top = pi->density.pressure_bar_dh * hi_over_density;
const float top_term_bottom =
hydro_gamma_minus_one * pj->mass * pj->u;
......@@ -622,4 +626,5 @@ __attribute__((always_inline)) INLINE static float hydro_h_term(
return f_ij;
}
#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_H */
#endif /* SWIFT_MINIMAL_HYDRO_H */
src/hydro/PressureEnergy/hydro_debug.h
+ 20
− 14
View file @ 0596dfcf
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
......@@ -18,29 +18,35 @@
******************************************************************************/
#ifndef SWIFT_PRESSURE_ENERGY_HYDRO_DEBUG_H
#define SWIFT_PRESSURE_ENERGY_HYDRO_DEBUG_H
/**
* @file PressureEnergy/hydro_debug.h
* @brief Pressure-Energy implementation of SPH (Debugging routines)
* @file Minimal/hydro_debug.h
* @brief Minimal conservative implementation of SPH (Debugging routines)
*
* The thermal variable is the energy (u) and the pressure is smoothed over
* contact discontinuities to prevent spurious surface tension.
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
*
* Follows equations (16), (17) and (18) of Hopkins, P., MNRAS, 2013,
* Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
* 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 Physics, 2012, Volume 231, Issue 3,
* pp. 759-794.
*/
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
const struct part* p, const struct xpart* xp) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"h=%.3e, wcount=%.3f, wcount_dh=%.3e, time_bin=%d\n"
"p_bar=%.3e, p_bar_dh=%3.e, u=%3.e\n",
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
"u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
"h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, \n"
"p_dh=%.3e, p_bar=%.3e \n"
"time_bin=%d\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
p->h, p->density.wcount, p->density.wcount_dh, p->time_bin,
p->pressure_bar, p->pressure_bar_dh, p->u);
p->u, p->u_dt, p->force.v_sig, hydro_get_comoving_pressure(p), p->h,
p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
p->density.pressure_bar_dh, p->pressure_bar,
p->time_bin);
}
#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_DEBUG_H */
#endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */
src/hydro/PressureEnergy/hydro_iact.h
+ 263
− 141
View file @ 0596dfcf
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* This file is part* of SWIFT.
* Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
......@@ -16,201 +16,323 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H
#define SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H
#ifndef SWIFT_MINIMAL_HYDRO_IACT_H
#define SWIFT_MINIMAL_HYDRO_IACT_H
/**
* @file PressureEnergy/hydro_iact.h
* @brief Pressure-Energy implementation of SPH (Neighbour loop equations)
* @file Minimal/hydro_iact.h
* @brief Minimal conservative implementation of SPH (Neighbour loop equations)
*
* The thermal variable is the energy (u) and the pressure is smoothed over
* contact discontinuities to prevent spurious surface tension.
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
*
* Follows equations (16), (17) and (18) of Hopkins, P., MNRAS, 2013,
* Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
* 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
* Physics, 2012, Volume 231, Issue 3, pp. 759-794.
*/
#include "adiabatic_index.h"
#include "minmax.h"
/**
* @brief Density loop (non-symmetric version)
* @brief Density interaction between two part*icles.
*
* @param r2 Comoving square distance between the two part*icles.
* @param dx Comoving vector separating both part*icles (pi - pj).
* @param hi Comoving smoothing-length of part*icle i.
* @param hj Comoving smoothing-length of part*icle j.
* @param pi First part*icle.
* @param pj Second part*icle.
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
float a, float H) {
__attribute__((always_inline)) INLINE static void runner_iact_density(
float r2, const float *dx, float hi, float hj, struct part* pi,
struct part* pj, float a, float H) {
float wi, wi_dx;
float dv[3], curlvr[3];
float wi, wj, wi_dx, wj_dx;
/* Get r and r inverse. */
const float r_inv = 1.0f / sqrtf(r2);
const float r = r2 * r_inv;
const float r = sqrtf(r2);
/* Get the masses and energies */
/* Get the masses. */
const float mi = pi->mass;
const float mj = pj->mass;
const float int_energy_j = pj->u;
/* Compute the kernel function */
const float hi_inv = 1.0f / hi;
const float hj_inv = 1.0f / hj;
/* Compute density of pi. */
const float hi_inv = 1.f / hi;
const float ui = r * hi_inv;
kernel_deval(ui, &wi, &wi_dx);
/* Compute contribution to the number of neighbours */
/* Note that wcount is \bar{n} */
pi->rho += mj * wi;
pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
pi->pressure_bar += mj * wi * pj->u;
pi->density.pressure_bar_dh -= mj * pj->u * (hydro_dimension * wi + ui * wi_dx);
pi->density.wcount += wi;
pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx) * hj_inv;
/* Compute the 'real' SPH density for other uses */
const float wimj = wi * mj;
pi->rho += wimj;
/* Compute the contribution to the local pressure */
/* The (gamma - 1) factor is added in hydro_end_density. */
pi->pressure_bar += int_energy_j * wimj;
pi->pressure_bar_dh -=
mj * int_energy_j * hj_inv * (hydro_dimension * wi + ui * wi_dx);
/* The following is lifted directly from the PressureEntropy code */
const float fac = mj * wi_dx * r_inv;
/* Compute dv dot r */
dv[0] = pi->v[0] - pj->v[0];
dv[1] = pi->v[1] - pj->v[1];
dv[2] = pi->v[2] - pj->v[2];
const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
pi->density.div_v -= fac * dvdr;
/* Compute dv cross r */
curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
pi->density.rot_v[0] += fac * curlvr[0];
pi->density.rot_v[1] += fac * curlvr[1];
pi->density.rot_v[2] += fac * curlvr[2];
pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
/* Compute density of pj. */
const float hj_inv = 1.f / hj;
const float uj = r * hj_inv;
kernel_deval(uj, &wj, &wj_dx);
pj->rho += mi * wj;
pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
pj->pressure_bar += mi * wj * pi->u;
pj->density.pressure_bar_dh -= mi * pi->u * (hydro_dimension * wj + uj * wj_dx);
pj->density.wcount += wj;
pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
}
/**
* @brief Density loop
* @brief Density interaction between two part*icles (non-symmetric).
*
* @param r2 Comoving square distance between the two part*icles.
* @param dx Comoving vector separating both part*icles (pi - pj).
* @param hi Comoving smoothing-length of part*icle i.
* @param hj Comoving smoothing-length of part*icle j.
* @param pi First part*icle.
* @param pj Second part*icle (not updated).
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_density(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
float a, float H) {
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
float r2, const float *dx, float hi, float hj, struct part* pi,
const struct part* pj, float a, float H) {
float wi, wi_dx;
/* Get the masses. */
const float mj = pj->mass;
/* For now, just do the simple case */
float negative_dx[3];
/* Get r and r inverse. */
const float r = sqrtf(r2);
negative_dx[0] = -1.f * dx[0];
negative_dx[1] = -1.f * dx[1];
negative_dx[2] = -1.f * dx[2];
const float h_inv = 1.f / hi;
const float ui = r * h_inv;
kernel_deval(ui, &wi, &wi_dx);
runner_iact_nonsym_density(r2, dx, hi, hj, pi, pj, a, H);
runner_iact_nonsym_density(r2, negative_dx, hj, hi, pj, pi, a, H);
pi->rho += mj * wi;
pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
pi->pressure_bar += mj * wi * pj->u;
pi->density.pressure_bar_dh -= mj * pj->u * (hydro_dimension * wi + ui * wi_dx);
pi->density.wcount += wi;
pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
}
/**
* @brief Force loop (non-symmetric version)
* @brief Force interaction between two part*icles.
*
* @param r2 Comoving square distance between the two part*icles.
* @param dx Comoving vector separating both part*icles (pi - pj).
* @param hi Comoving smoothing-length of part*icle i.
* @param hj Comoving smoothing-length of part*icle j.
* @param pi First part*icle.
* @param pj Second part*icle.
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
float a, float H) {
float wi, wi_dx, wj, wj_dx;
const float fac_mu = 1.f; /* Will change with cosmological integration */
float dv[3];
/* Masses */
const float mj = pj->mass;
__attribute__((always_inline)) INLINE static void runner_iact_force(
float r2, const float *dx, float hi, float hj, struct part* pi,
struct part* pj, float a, float H) {
/* Get r and r inverse */
const float r = sqrtf(r2);
const float ri = 1.0f / r;
/* Cosmological factors entering the EoMs */
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
/* Get velocity difference */
dv[0] = pi->v[0] - pj->v[0];
dv[1] = pi->v[1] - pj->v[1];
dv[2] = pi->v[2] - pj->v[2];
const float r = sqrtf(r2);
const float r_inv = 1.0f / r;
const float dvdr = dv[0] * dx[0] +
dv[1] * dx[1] +
dv[2] * dx[2];
/* Recover some data */
const float mi = pi->mass;
const float mj = pj->mass;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
/* Compute f terms */
const float f_ij = hydro_h_term(pi, pj, hi);
const float f_ji = hydro_h_term(pj, pi, hj);
/* Compute kernel function */
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
const float xi = r * hi_inv;
float wi, wi_dx;
kernel_deval(xi, &wi, &wi_dx);
const float wi_dr = hid_inv * wi_dx;
/* Get the kernel for hj. */
const float hj_inv = 1.0f / hj;
const float hid_inv = pow_dimension_plus_one(hi_inv);
const float hjd_inv = pow_dimension_plus_one(hj_inv);
const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
const float xj = r * hj_inv;
float wj, wj_dx;
kernel_deval(xj, &wj, &wj_dx);
const float wj_dr = hjd_inv * wj_dx;
const float ui = r * hi_inv;
const float uj = r * hj_inv;
kernel_deval(ui, &wi, &wi_dx);
kernel_deval(uj, &wj, &wj_dx);
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
(pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
const float wi_dr = hid_inv * wi_dx;
const float wj_dr = hjd_inv * wj_dx;
/* Are the part*icles moving towards each others ? */
const float omega_ij = min(dvdr, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute gradient terms */
/* Compute sound speeds and signal velocity */
const float ci = pi->force.soundspeed;
const float cj = pj->force.soundspeed;
const float v_sig = ci + cj - 3.f * mu_ij;
/* Construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
/* SPH acceleration term */
const float sph_acc_term = pj->u * pi->u *
hydro_gamma_minus_one * hydro_gamma_minus_one *
((f_ij/pi->pressure_bar) * wi_dr + (f_ji/pj->pressure_bar) * wj_dr) * r_inv;
/* Assemble the acceleration */
const float acc = sph_acc_term + visc_acc_term;
/* Use the force Luke ! */
pi->a_hydro[0] -= mj * acc * dx[0];
pi->a_hydro[1] -= mj * acc * dx[1];
pi->a_hydro[2] -= mj * acc * dx[2];
pj->a_hydro[0] += mi * acc * dx[0];
pj->a_hydro[1] += mi * acc * dx[1];
pj->a_hydro[2] += mi * acc * dx[2];
/* Get the time derivative for u. */
const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
pj->u * pi->u * (f_ij / pi->pressure_bar) * wi_dr * dvdr * r_inv;
const float sph_du_term_j = hydro_gamma_minus_one * hydro_gamma_minus_one *
pi->u * pj->u * (f_ji / pj->pressure_bar) * wj_dr * dvdr * r_inv;
/* Viscosity term */
const float visc_du_term = 0.5f * visc_acc_term * dvdr;
/* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term;
const float du_dt_j = sph_du_term_j + visc_du_term;
/* Internal energy time derivatibe */
pi->u_dt += du_dt_i * mj;
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;
/* Update the signal velocity. */
pi->force.v_sig = max(pi->force.v_sig, v_sig);
pj->force.v_sig = max(pj->force.v_sig, v_sig);
}
/**
* @brief Force interaction between two part*icles (non-symmetric).
*
* @param r2 Comoving square distance between the two part*icles.
* @param dx Comoving vector separating both part*icles (pi - pj).
* @param hi Comoving smoothing-length of part*icle i.
* @param hj Comoving smoothing-length of part*icle j.
* @param pi First part*icle.
* @param pj Second part*icle (not updated).
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
float r2, const float *dx, float hi, float hj, struct part* pi,
const struct part* pj, float a, float H) {
/* Cosmological factors entering the EoMs */
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
const float r = sqrtf(r2);
const float r_inv = 1.0f / r;
/* Recover some data */
// const float mi = pi->mass;
const float mj = pj->mass;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
/* Compute gradient terms */
const float f_ij = hydro_h_term(pi, pj, hi);
const float f_ji = hydro_h_term(pj, pi, hj);
/* Artificial viscosity */
const float omega_ij = fminf(dvdr, 0.f);
const float mu_ij = fac_mu * ri * omega_ij;
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
const float xi = r * hi_inv;
float wi, wi_dx;
kernel_deval(xi, &wi, &wi_dx);
const float wi_dr = hid_inv * wi_dx;
const float v_sig = pi->force.soundspeed + pj->force.soundspeed - 3.f * mu_ij;
/* Get the kernel for hj. */
const float hj_inv = 1.0f / hj;
const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
const float xj = r * hj_inv;
float wj, wj_dx;
kernel_deval(xj, &wj, &wj_dx);
const float wj_dr = hjd_inv * wj_dx;
const float rho_ij = 0.5f * (pi->rho + pj->rho);
const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
(pi->force.balsara + pj->force.balsara) / rho_ij;
const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
(pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
/* Calculate the equation of motion H13:15 */
/* Are the part*icles moving towards each others ? */
const float omega_ij = min(dvdr, 0.f);
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
const float sph_term = mj * pj->u * pi->u *
hydro_gamma_minus_one * hydro_gamma_minus_one *
((f_ij/pi->pressure_bar) * wi_dr + (f_ji/pj->pressure_bar) * wj_dr);
/* Compute sound speeds and signal velocity */
const float ci = pi->force.soundspeed;
const float cj = pj->force.soundspeed;
const float v_sig = ci + cj - 3.f * mu_ij;
const float acc = (visc_term + sph_term) * ri;
/* Construct the full viscosity term */
const float rho_ij = 0.5f * (rhoi + rhoj);
const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
/* Use the force Luke ! */
pi->a_hydro[0] -= acc * dx[0];
pi->a_hydro[1] -= acc * dx[1];
pi->a_hydro[2] -= acc * dx[2];
/* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
/* Get the time derivative for h */
pi->force.h_dt -= mj * dvdr * ri / pj->rho * wi_dr;
/* SPH acceleration term */
const float sph_acc_term = pj->u * pi->u *
hydro_gamma_minus_one * hydro_gamma_minus_one *
((f_ij/pi->pressure_bar) * wi_dr + (f_ji/pj->pressure_bar) * wj_dr) * r_inv;
/* Update the signal velocity */
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
/* Assemble the acceleration */
const float acc = sph_acc_term + visc_acc_term;
/* Calcualte change in energy from viscosity */
const float u_dt_visc = mj * visc_term * dvdr * ri;
/* Use the force Luke ! */
pi->a_hydro[0] -= mj * acc * dx[0];
pi->a_hydro[1] -= mj * acc * dx[1];
pi->a_hydro[2] -= mj * acc * dx[2];
/* Calculate the change in internal energy from hydro H13:16 */
const float u_dt_hydro =
hydro_gamma_minus_one * hydro_gamma_minus_one *
mj * pj->u * pi->u * (f_ij / pi->pressure_bar) * wi_dr * dvdr * ri;
/* Get the time derivative for u. */
const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
pj->u * pi->u * (f_ij / pi->pressure_bar) * wi_dr * dvdr * r_inv;
pi->u_dt += (u_dt_hydro + u_dt_visc);
}
/* Viscosity term */
const float visc_du_term = 0.5f * visc_acc_term * dvdr;
/**
* @brief Force loop
*/
__attribute__((always_inline)) INLINE static void runner_iact_force(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
float a, float H) {
/* Assemble the energy equation term */
const float du_dt_i = sph_du_term_i + visc_du_term;
/* For now, just do the simple case */
float negative_dx[3];
/* Internal energy time derivatibe */
pi->u_dt += du_dt_i * mj;
negative_dx[0] = -1.f * dx[0];
negative_dx[1] = -1.f * dx[1];
negative_dx[2] = -1.f * dx[2];
/* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
runner_iact_nonsym_force(r2, dx, hi, hj, pi, pj, a, H);
runner_iact_nonsym_force(r2, negative_dx, hj, hi, pj, pi, a, H);
/* Update the signal velocity. */
pi->force.v_sig = max(pi->force.v_sig, v_sig);
}
#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H */
#endif /* SWIFT_MINIMAL_HYDRO_IACT_H */
src/hydro/PressureEnergy/hydro_io.h
+ 50
− 41
View file @ 0596dfcf
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
......@@ -16,18 +16,20 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_PRESSURE_ENERGY_HYDRO_IO_H
#define SWIFT_PRESSURE_ENERGY_HYDRO_IO_H
#ifndef SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H
#define SWIFT_PRESSURE_ENERGY_HYDRO_IACT_H
/**
* @file PressureEnergy/hydro_io.h
* @brief Pressure-Energy implementation of SPH (i/o routines)
* @file Minimal/hydro_io.h
* @brief Minimal conservative implementation of SPH (i/o routines)
*
* The thermal variable is the energy (u) and the pressure is smoothed over
* contact discontinuities to prevent spurious surface tension.
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
*
* Follows equations (17) and (18) of Hopkins, P., MNRAS, 2013,
* Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
* 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 Physics, 2012, Volume 231, Issue 3,
* pp. 759-794.
*/
#include "adiabatic_index.h"
......@@ -52,21 +54,30 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_LENGTH, parts, x);
list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
UNIT_CONV_SPEED, parts, v);
list[2] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
parts, mass);
list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
UNIT_CONV_LENGTH, parts, h);
list[3] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
parts, mass);
list[4] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
UNIT_CONV_NO_UNITS, parts, id);
list[5] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
UNIT_CONV_ENERGY_PER_UNIT_MASS, parts,
u);
list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
UNIT_CONV_ACCELERATION, parts, a_hydro);
list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
UNIT_CONV_DENSITY, parts, rho);
}
void convert_S(const struct engine* e, const struct part* p, float* ret) {
ret[0] = hydro_get_comoving_entropy(p);
}
void convert_P(const struct engine* e, const struct part* p, float* ret) {
ret[0] = hydro_get_comoving_pressure(p);
}
void convert_part_pos(const struct engine* e, const struct part* p,
double* ret) {
......@@ -94,26 +105,24 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
*num_fields = 9;
/* List what we want to write */
list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
UNIT_CONV_LENGTH, parts,
convert_part_pos);
list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED,
parts, v);
list[2] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
list[0] = io_make_output_field_convert_part(
"Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
list[1] =
io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
list[2] =
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
parts, h);
list[3] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
list[4] = io_make_output_field("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
parts, u);
list[5] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts,
mass);
list[6] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
rho);
list[7] = io_make_output_field("Accelerations", FLOAT, 3,
UNIT_CONV_ACCELERATION, parts, a_hydro);
list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
parts, pressure_bar);
UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
list[6] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
list[7] = io_make_output_field_convert_part(
"Entropy", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, convert_S);
list[8] = io_make_output_field_convert_part(
"Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
}
/**
......@@ -121,12 +130,12 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
* @param h_grpsph The HDF5 group in which to write
*/
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");
io_write_attribute_s(
h_grpsph, "Viscosity Model",
"as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
io_write_attribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
io_write_attribute_f(h_grpsph, "Viscosity beta", 3.f);
io_write_attribute_s(h_grpsph, "Viscosity Model",
"Minimal treatment as in Monaghan (1992)");
/* Time integration properties */
io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
......@@ -140,4 +149,4 @@ void hydro_write_flavour(hid_t h_grpsph) {
*/
int writeEntropyFlag() { return 0; }
#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_IO_H */
#endif /* SWIFT_MINIMAL_HYDRO_IO_H */
src/hydro/PressureEnergy/hydro_part.h
+ 69
− 55
View file @ 0596dfcf
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
......@@ -18,37 +18,44 @@
******************************************************************************/
#ifndef SWIFT_PRESSURE_ENERGY_HYDRO_PART_H
#define SWIFT_PRESSURE_ENERGY_HYDRO_PART_H
/**
* @file PressureEnergy/hydro_part.h
* @brief Pressure-Energy implementation of SPH (Particle definition)
* @file Minimal/hydro_part.h
* @brief Minimal conservative implementation of SPH (Particle definition)
*
* The thermal variable is the energy (u) and the pressure is smoothed over
* contact discontinuities to prevent spurious surface tension.
* The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction
* term is implemented.
*
* Follows equations (16), (17) and (18) of Hopkins, P., MNRAS, 2013,
* Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
* 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
* Physics, 2012, Volume 231, Issue 3, pp. 759-794.
*/
#include "chemistry_struct.h"
#include "cooling_struct.h"
/* Extra particle data not needed during the SPH loops over neighbours. */
/**
* @brief Particle fields not needed during the SPH loops over neighbours.
*
* This structure contains the particle fields that are not used in the
* density or force loops. Quantities should be used in the kick, drift and
* potentially ghost tasks only.
*/
struct xpart {
/*! Offset between current position and position at last tree rebuild. */
float x_diff[3];
/* Offset between the current position and position at the last sort. */
/*! Offset between the current position and position at the last sort. */
float x_diff_sort[3];
/*! Velocity at the last full step. */
float v_full[3];
/*! Gravitational acceleration */
/*! Gravitational acceleration at the last full step. */
float a_grav[3];
/*! Full internal energy */
/*! Internal energy at the last full step. */
float u_full;
/*! Additional data used to record cooling information */
......@@ -56,10 +63,16 @@ struct xpart {
} SWIFT_STRUCT_ALIGN;
/* Data of a single particle. */
/**
* @brief Particle fields for the SPH particles
*
* The density and force substructures are used to contain variables only used
* within the density and force loops over neighbours. All more permanent
* variables should be declared in the main part of the part structure,
*/
struct part {
/*! Particle ID. */
/*! Particle unique ID. */
long long id;
/*! Pointer to corresponding gravity part. */
......@@ -77,77 +90,78 @@ struct part {
/*! Particle mass. */
float mass;
/*! SPH Density */
float rho;
/*! Smoothed particle pressure. */
float pressure_bar;
/*! Smoothed particle pressure's spatial derivative */
float pressure_bar_dh;
/*! Particle smoothing length. */
float h;
/*! Particle internal energy */
/*! Particle internal energy. */
float u;
/*! Differential of the internal energy with respect to time */
/*! Time derivative of the internal energy. */
float u_dt;
/*! Entropy (for if people want it) */
float entropy;
/*! Particle cutoff radius. */
float h;
/*! Particle density. */
float rho;
/*! Particle pressure (weighted) */
float pressure_bar;
/* Store density/force specific stuff. */
/*union {*/
/**
* @brief Structure for the variables only used in the density loop over
* neighbours.
*
* Quantities in this sub-structure should only be accessed in the density
* loop over neighbours and the ghost task.
*/
struct {
/*! Number of neighbours. */
/*! Neighbour number count. */
float wcount;
/*! Number of neighbours spatial derivative. */
/*! Derivative of the neighbour number with respect to h. */
float wcount_dh;
/*! Velocity curl */
float rot_v[3];
/*! Velocity divergence */
float div_v;
/*! d\rho/dh */
/*! Derivative of density with respect to h */
float rho_dh;
/*! dP/dh */
float pressure_dh;
/*! Derivative of the weighted pressure with respect to h */
float pressure_bar_dh;
} density;
/**
* @brief Structure for the variables only used in the force loop over
* neighbours.
*
* Quantities in this sub-structure should only be accessed in the force
* loop over neighbours and the ghost, drift and kick tasks.
*/
struct {
/*! Signal velocity. */
float v_sig;
/*! "Grad h" term */
float f;
/*! Time derivative of the smoothing length */
float h_dt;
/*! Particle pressure. */
float pressure;
/*! Sound speed */
/*! Particle soundspeed. */
float soundspeed;
/*! Balsara switch */
float balsara;
/*! F_{ij} -- not actually possible to set this. */
float f;
/*! Particle signal velocity */
float v_sig;
/*! P/\rho^2 -- not actually required for Pressure Energy but this is
needed for cross-compatibility with the unit tests */
float P_over_rho2;
/*! Time derivative of smoothing length */
float h_dt;
} force;
/*};*/
/* Chemistry information */
struct chemistry_part_data chemistry_data;
/* Time-step length */
/*! Time-step length */
timebin_t time_bin;
#ifdef SWIFT_DEBUG_CHECKS
......@@ -162,4 +176,4 @@ struct part {
} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_PRESSURE_ENERGY_HYDRO_PART_H */
#endif /* SWIFT_MINIMAL_HYDRO_PART_H */
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment