Commit 489556aa authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'gamma_power_speed_up' into 'master'

Gamma power speed up

"Fix" to #185. See also the discussion in #166.

Two improvements:

 - Speed-up the calculation of x^\gamma, x^{\gamma-1} and x^-(\gamma-1) . Currently support \gamma=5/3, \gamma=4/3 and \gamma=2.
 - Write the details of the hydro scheme in the snapshots.

Branch can be removed.

See merge request !195
parents 42808b91 fa2b946e
......@@ -1428,7 +1428,7 @@ FORMULA_TRANSPARENT = YES
# The default value is: NO.
# This tag requires that the tag GENERATE_HTML is set to YES.
USE_MATHJAX = NO
USE_MATHJAX = YES
# When MathJax is enabled you can set the default output format to be used for
# the MathJax output. See the MathJax site (see:
......
......@@ -54,11 +54,11 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
# Include files for distribution, not installation.
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
vector.h runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h kick.h \
timestep.h drift.h \
timestep.h drift.h adiabatic_index.h \
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
hydro.h hydro_io.h \
hydro.h hydro_io.h \
hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \
hydro/Default/hydro.h hydro/Default/hydro_iact.h hydro/Default/hydro_io.h \
......
/*******************************************************************************
* 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
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_ADIABATIC_INDEX_H
#define SWIFT_ADIABATIC_INDEX_H
/* Config parameters. */
#include "../config.h"
/* Some standard headers. */
#include <math.h>
/* Local headers. */
#include "const.h"
#include "debug.h"
#include "inline.h"
/* First define some constants */
#if defined(HYDRO_GAMMA_5_3)
#define hydro_gamma 1.66666666666666667f
#define hydro_gamma_minus_one 0.66666666666666667f
#define hydro_one_over_gamma_minus_one 1.5f
#elif defined(HYDRO_GAMMA_4_3)
#define hydro_gamma 1.33333333333333333f
#define hydro_gamma_minus_one 0.33333333333333333f
#define hydro_one_over_gamma_minus_one 3.f
#elif defined(HYDRO_GAMMA_2_1)
#define hydro_gamma 2.f
#define hydro_gamma_minus_one 1.f
#define hydro_one_over_gamma_minus_one 1.f
#endif
/**
* @brief Returns the argument to the power given by the adiabatic index
*
* Computes \f$x^\gamma\f$.
*/
__attribute__((always_inline)) INLINE static float pow_gamma(float x) {
#if defined(HYDRO_GAMMA_5_3)
const float x2 = x * x;
const float x5 = x2 * x2 * x;
return cbrtf(x5);
#elif defined(HYDRO_GAMMA_4_3)
const float x2 = x * x;
const float x4 = x2 * x2;
return cbrtf(x4);
#elif defined(HYDRO_GAMMA_2_1)
return x * x;
#else
error("The adiabatic index is not defined !");
return 0.f;
#endif
}
/**
* @brief Returns the argument to the power given by the adiabatic index minus
* one
*
* Computes \f$x^{(\gamma-1)}\f$.
*/
__attribute__((always_inline)) INLINE static float pow_gamma_minus_one(
float x) {
#if defined(HYDRO_GAMMA_5_3)
return cbrtf(x * x);
#elif defined(HYDRO_GAMMA_4_3)
return cbrtf(x);
#elif defined(HYDRO_GAMMA_2_1)
return x;
#else
error("The adiabatic index is not defined !");
return 0.f;
#endif
}
/**
* @brief Returns one over the argument to the power given by the adiabatic
* index minus one
*
* Computes \f$x^{-(\gamma-1)}\f$.
*/
__attribute__((always_inline)) INLINE static float pow_minus_gamma_minus_one(
float x) {
#if defined(HYDRO_GAMMA_5_3)
return 1.f / cbrtf(x * x);
#elif defined(HYDRO_GAMMA_4_3)
return 1.f / cbrtf(x);
#elif defined(HYDRO_GAMMA_2_1)
return 1.f / x;
#else
error("The adiabatic index is not defined !");
return 0.f;
#endif
}
#endif /* SWIFT_ADIABATIC_INDEX_H */
......@@ -20,9 +20,6 @@
#ifndef SWIFT_CONST_H
#define SWIFT_CONST_H
/* Hydrodynamical constants. */
#define const_hydro_gamma (5.0f / 3.0f)
/* SPH Viscosity constants. */
#define const_viscosity_alpha 0.8f
#define const_viscosity_alpha_min \
......@@ -44,6 +41,11 @@
#define const_G 6.672e-8f /* Gravitational constant. */
#define const_epsilon 0.0014f /* Gravity blending distance. */
/* Hydrodynamical adiabatic index. */
#define HYDRO_GAMMA_5_3
//#define HYDRO_GAMMA_4_3
//#define HYDRO_GAMMA_2_1
/* Kernel function to use */
#define CUBIC_SPLINE_KERNEL
//#define QUARTIC_SPLINE_KERNEL
......
......@@ -17,6 +17,7 @@
*
******************************************************************************/
#include "adiabatic_index.h"
#include "approx_math.h"
/**
......@@ -138,12 +139,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute this particle's sound speed. */
const float u = p->u;
const float fc = p->force.c =
sqrtf(const_hydro_gamma * (const_hydro_gamma - 1.0f) * u);
const float fc = p->force.c = sqrtf(hydro_gamma * hydro_gamma_minus_one * u);
/* Compute the P/Omega/rho2. */
xp->omega = 1.0f + 0.3333333333f * h * p->rho_dh / p->rho;
p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (p->rho * xp->omega);
p->force.POrho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega);
/* Balsara switch */
p->force.balsara = normDiv_v / (normDiv_v + normCurl_v + 0.0001f * fc * ih);
......@@ -207,7 +207,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
u = p->u *= expf(w);
/* Predict gradient term */
p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (p->rho * xp->omega);
p->force.POrho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega);
}
/**
......
......@@ -20,6 +20,8 @@
#ifndef SWIFT_RUNNER_IACT_H
#define SWIFT_RUNNER_IACT_H
#include "adiabatic_index.h"
/**
* @brief SPH interaction functions following the Gadget-2 version of SPH.
*
......@@ -408,7 +410,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
Pi_ij *= (pi->force.balsara + pj->force.balsara);
/* Thermal conductivity */
v_sig_u = sqrtf(2.f * (const_hydro_gamma - 1.f) *
v_sig_u = sqrtf(2.f * hydro_gamma_minus_one *
fabs(rhoi * pi->u - rhoj * pj->u) / (rhoi + rhoj));
tc = const_conductivity_alpha * v_sig_u / (rhoi + rhoj);
tc *= (wi_dr + wj_dr);
......@@ -608,7 +610,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
Pi_ij.v *= (wi_dr.v + wj_dr.v);
/* Thermal conductivity */
v_sig_u.v = vec_sqrt(vec_set1(2.f * (const_hydro_gamma - 1.f)) *
v_sig_u.v = vec_sqrt(vec_set1(2.f * hydro_gamma_minus_one) *
vec_fabs(pirho.v * piu.v - pjrho.v * pju.v) /
(pirho.v + pjrho.v));
tc.v = vec_set1(const_conductivity_alpha) * v_sig_u.v / (pirho.v + pjrho.v);
......@@ -721,7 +723,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
Pi_ij *= (pi->force.balsara + pj->force.balsara);
/* Thermal conductivity */
v_sig_u = sqrtf(2.f * (const_hydro_gamma - 1.f) *
v_sig_u = sqrtf(2.f * hydro_gamma_minus_one *
fabs(rhoi * pi->u - rhoj * pj->u) / (rhoi + rhoj));
tc = const_conductivity_alpha * v_sig_u / (rhoi + rhoj);
tc *= (wi_dr + wj_dr);
......@@ -912,7 +914,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
Pi_ij.v *= (wi_dr.v + wj_dr.v);
/* Thermal conductivity */
v_sig_u.v = vec_sqrt(vec_set1(2.f * (const_hydro_gamma - 1.f)) *
v_sig_u.v = vec_sqrt(vec_set1(2.f * hydro_gamma_minus_one) *
vec_fabs(pirho.v * piu.v - pjrho.v * pju.v) /
(pirho.v + pjrho.v));
tc.v = vec_set1(const_conductivity_alpha) * v_sig_u.v / (pirho.v + pjrho.v);
......
......@@ -102,9 +102,6 @@ __attribute__((always_inline)) INLINE static void hydro_write_particles(
*/
void writeSPHflavour(hid_t h_grpsph) {
/* Kernel function description */
writeAttribute_s(h_grpsph, "Kernel", kernel_name);
/* Viscosity and thermal conduction */
writeAttribute_s(h_grpsph, "Thermal Conductivity Model",
"Price (2008) without switch");
......
......@@ -17,6 +17,8 @@
*
******************************************************************************/
#include "adiabatic_index.h"
/**
* @brief Computes the hydro time-step of a given particle
*
......@@ -132,11 +134,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the pressure */
const float dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
p->force.pressure =
(p->entropy + p->entropy_dt * dt) * powf(p->rho, const_hydro_gamma);
p->force.pressure = (p->entropy + p->entropy_dt * dt) * pow_gamma(p->rho);
/* Compute the sound speed */
p->force.soundspeed = sqrtf(const_hydro_gamma * p->force.pressure / p->rho);
p->force.soundspeed = sqrtf(hydro_gamma * p->force.pressure / p->rho);
}
/**
......@@ -179,10 +180,10 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Drift the pressure */
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
p->force.pressure =
(p->entropy + p->entropy_dt * dt_entr) * powf(p->rho, const_hydro_gamma);
(p->entropy + p->entropy_dt * dt_entr) * pow_gamma(p->rho);
/* Compute the new sound speed */
p->force.soundspeed = sqrtf(const_hydro_gamma * p->force.pressure / p->rho);
p->force.soundspeed = sqrtf(hydro_gamma * p->force.pressure / p->rho);
}
/**
......@@ -195,8 +196,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part* p) {
p->entropy_dt *=
(const_hydro_gamma - 1.f) * powf(p->rho, -(const_hydro_gamma - 1.f));
p->entropy_dt *= hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
}
/**
......@@ -232,8 +232,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part* p) {
p->entropy = (const_hydro_gamma - 1.f) * p->entropy *
powf(p->rho, -(const_hydro_gamma - 1.f));
p->entropy =
hydro_gamma_minus_one * p->entropy * pow_minus_gamma_minus_one(p->rho);
}
/**
......@@ -247,6 +247,5 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const float entropy = p->entropy + p->entropy_dt * dt;
return entropy * powf(p->rho, const_hydro_gamma - 1.f) *
(1.f / (const_hydro_gamma - 1.f));
return entropy * pow_gamma_minus_one(p->rho) * hydro_one_over_gamma_minus_one;
}
......@@ -102,9 +102,6 @@ __attribute__((always_inline)) INLINE static void hydro_write_particles(
*/
void writeSPHflavour(hid_t h_grpsph) {
/* Kernel function description */
writeAttribute_s(h_grpsph, "Kernel", kernel_name);
/* Viscosity and thermal conduction */
writeAttribute_s(h_grpsph, "Thermal Conductivity Model",
"(No treatment) Legacy Gadget-2 as in Springel (2005)");
......
......@@ -17,6 +17,7 @@
*
******************************************************************************/
#include "adiabatic_index.h"
#include "approx_math.h"
/**
......@@ -132,7 +133,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* p, struct xpart* xp, int ti_current, double timeBase) {
p->force.pressure = p->rho * p->u * (const_hydro_gamma - 1.f);
p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
}
/**
......@@ -175,7 +176,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->u = xp->u_full;
/* Need to recompute the pressure as well */
p->force.pressure = p->rho * p->u * (const_hydro_gamma - 1.f);
p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
}
/**
......
......@@ -19,6 +19,8 @@
#ifndef SWIFT_RUNNER_IACT_MINIMAL_H
#define SWIFT_RUNNER_IACT_MINIMAL_H
#include "adiabatic_index.h"
/**
* @brief Minimal conservative implementation of SPH
*
......@@ -150,8 +152,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
(pi->v[2] - pj->v[2]) * dx[2];
/* Compute sound speeds */
const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
const float v_sig = ci + cj;
/* SPH acceleration term */
......@@ -233,8 +235,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
(pi->v[2] - pj->v[2]) * dx[2];
/* Compute sound speeds */
const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
const float v_sig = ci + cj;
/* SPH acceleration term */
......
......@@ -102,10 +102,6 @@ __attribute__((always_inline)) INLINE static void hydro_write_particles(
*/
void writeSPHflavour(hid_t h_grpsph) {
/* Kernel function description */
writeAttribute_s(h_grpsph, "Kernel", kernel_name);
writeAttribute_f(h_grpsph, "Hydro gamma", const_hydro_gamma);
/* Viscosity and thermal conduction */
/* Nothing in this minimal model... */
writeAttribute_s(h_grpsph, "Thermal Conductivity Model", "No model");
......
......@@ -25,6 +25,8 @@
#include <math.h>
/* Local headers. */
#include "adiabatic_index.h"
#include "common_io.h"
#include "error.h"
#include "hydro.h"
#include "kernel_hydro.h"
......@@ -54,6 +56,7 @@ void hydro_props_init(struct hydro_props *p,
void hydro_props_print(const struct hydro_props *p) {
message("Adiabatic index gamma: %f.", hydro_gamma);
message("Hydrodynamic scheme: %s.", SPH_IMPLEMENTATION);
message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours (eta=%f).",
kernel_name, p->target_neighbours, p->delta_neighbours,
......@@ -68,3 +71,21 @@ void hydro_props_print(const struct hydro_props *p) {
message("Maximal iterations in ghost task set to %d (default is %d)",
p->max_smoothing_iterations, hydro_props_default_max_iterations);
}
#if defined(HAVE_HDF5)
void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
writeAttribute_f(h_grpsph, "Adiabatic index", hydro_gamma);
writeAttribute_s(h_grpsph, "Scheme", SPH_IMPLEMENTATION);
writeAttribute_s(h_grpsph, "Kernel function", kernel_name);
writeAttribute_f(h_grpsph, "Kernel target N_ngb", p->target_neighbours);
writeAttribute_f(h_grpsph, "Kernel delta N_ngb", p->delta_neighbours);
writeAttribute_f(h_grpsph, "Kernel eta", p->eta_neighbours);
writeAttribute_f(h_grpsph, "CFL parameter", p->CFL_condition);
writeAttribute_f(h_grpsph, "Volume log(max(delta h))", p->log_max_h_change);
writeAttribute_f(h_grpsph, "Volume max change time-step",
powf(expf(p->log_max_h_change), 3.f));
writeAttribute_f(h_grpsph, "Max ghost iterations",
p->max_smoothing_iterations);
}
#endif
......@@ -23,6 +23,10 @@
/* Config parameters. */
#include "../config.h"
#if defined(HAVE_HDF5)
#include <hdf5.h>
#endif
/* Local includes. */
#include "const.h"
#include "parser.h"
......@@ -53,4 +57,8 @@ struct hydro_props {
void hydro_props_print(const struct hydro_props *p);
void hydro_props_init(struct hydro_props *p, const struct swift_params *params);
#if defined(HAVE_HDF5)
void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p);
#endif
#endif /* SWIFT_HYDRO_PROPERTIES */
......@@ -39,7 +39,7 @@
#include "common_io.h"
#include "engine.h"
#include "error.h"
#include "kernel_hydro.h"
#include "hydro_properties.h"
#include "part.h"
#include "units.h"
......@@ -639,8 +639,10 @@ void write_output_parallel(struct engine* e, const char* baseName,
writeCodeDescription(h_file);
/* Print the SPH parameters */
h_grp = H5Gcreate(h_file, "/SPH", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
h_grp =
H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (h_grp < 0) error("Error while creating SPH group");
hydro_props_print_snapshot(h_grp, e->hydro_properties);
writeSPHflavour(h_grp);
H5Gclose(h_grp);
......
......@@ -39,7 +39,7 @@
#include "common_io.h"
#include "engine.h"
#include "error.h"
#include "kernel_hydro.h"
#include "hydro_properties.h"
#include "part.h"
#include "units.h"
......@@ -714,8 +714,10 @@ void write_output_serial(struct engine* e, const char* baseName,
writeCodeDescription(h_file);
/* Print the SPH parameters */
h_grp = H5Gcreate(h_file, "/SPH", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
h_grp = H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT,
H5P_DEFAULT);
if (h_grp < 0) error("Error while creating SPH group");
hydro_props_print_snapshot(h_grp, e->hydro_properties);
writeSPHflavour(h_grp);
H5Gclose(h_grp);
......
......@@ -38,7 +38,7 @@
#include "common_io.h"
#include "engine.h"
#include "error.h"
#include "kernel_hydro.h"
#include "hydro_properties.h"
#include "part.h"
#include "units.h"
......@@ -564,8 +564,10 @@ void write_output_single(struct engine* e, const char* baseName,
writeCodeDescription(h_file);
/* Print the SPH parameters */
h_grp = H5Gcreate(h_file, "/SPH", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
h_grp =
H5Gcreate(h_file, "/HydroScheme", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (h_grp < 0) error("Error while creating SPH group");
hydro_props_print_snapshot(h_grp, e->hydro_properties);
writeSPHflavour(h_grp);
H5Gclose(h_grp);
......
......@@ -37,6 +37,7 @@
#include "units.h"
/* Includes. */
#include "adiabatic_index.h"
#include "const.h"
#include "error.h"
......@@ -188,14 +189,14 @@ void units_get_base_unit_exponants_array(float baseUnitsExp[5],
break;
case UNIT_CONV_ENTROPY:
baseUnitsExp[UNIT_MASS] = 1.f - const_hydro_gamma;
baseUnitsExp[UNIT_LENGTH] = 3.f * const_hydro_gamma - 1.f;
baseUnitsExp[UNIT_MASS] = 1.f - hydro_gamma;
baseUnitsExp[UNIT_LENGTH] = 3.f * hydro_gamma - 1.f;
baseUnitsExp[UNIT_TIME] = -2.f;
break;
case UNIT_CONV_ENTROPY_PER_UNIT_MASS:
baseUnitsExp[UNIT_MASS] = -const_hydro_gamma;
baseUnitsExp[UNIT_LENGTH] = 3.f * const_hydro_gamma - 1.f;
baseUnitsExp[UNIT_MASS] = -hydro_gamma;
baseUnitsExp[UNIT_LENGTH] = 3.f * hydro_gamma - 1.f;
baseUnitsExp[UNIT_TIME] = -2.f;
break;
......
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