Commit 5736f328 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'master' into gravity_multi_dt

parents ef469ebb ab6d8549
......@@ -73,6 +73,9 @@ tests/testRiemannExact
tests/testRiemannTRRS
tests/testRiemannHLLC
tests/testMatrixInversion
tests/testVoronoi1D
tests/testVoronoi2D
tests/testVoronoi3D
tests/testDump
tests/testLogger
tests/benchmarkInteractions
......
......@@ -460,7 +460,6 @@ fi
# Check for HDF5. This is required.
AX_LIB_HDF5
if test "$with_hdf5" != "yes"; then
AC_MSG_ERROR([Could not find a working HDF5 library])
fi
......@@ -468,8 +467,6 @@ fi
# We want to know if this HDF5 supports MPI and whether we should use it.
# The default is to use MPI support if it is available, i.e. this is
# a parallel HDF5.
# To do this need to ask the HDF5 compiler about its configuration,
# -showconfig should have yes/no.
have_parallel_hdf5="no"
if test "$with_hdf5" = "yes"; then
AC_ARG_ENABLE([parallel-hdf5],
......@@ -482,7 +479,14 @@ if test "$with_hdf5" = "yes"; then
if test "$enable_parallel_hdf5" = "yes"; then
AC_MSG_CHECKING([for HDF5 parallel support])
parallel=`$H5CC -showconfig | grep "Parallel HDF5:" | awk '{print $3}'`
# Check if the library is capable, the header should define H5_HAVE_PARALLEL.
AC_COMPILE_IFELSE([AC_LANG_SOURCE([[
#include "hdf5.h"
#ifndef H5_HAVE_PARALLEL
# error macro not defined
#endif
]])], [parallel="yes"], [parallel="no"])
if test "$parallel" = "yes"; then
have_parallel_hdf5="yes"
AC_DEFINE([HAVE_PARALLEL_HDF5],1,[HDF5 library supports parallel access])
......@@ -583,7 +587,7 @@ fi
# Hydro scheme.
AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>],
[Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo default: gadget2@:>@]
[Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo, shadowfax default: gadget2@:>@]
)],
[with_hydro="$withval"],
[with_hydro="gadget2"]
......@@ -604,6 +608,9 @@ case "$with_hydro" in
gizmo)
AC_DEFINE([GIZMO_SPH], [1], [GIZMO SPH])
;;
shadowfax)
AC_DEFINE([SHADOWFAX_SPH], [1], [Shadowfax SPH])
;;
*)
AC_MSG_ERROR([Unknown hydrodynamics scheme: $with_hydro])
......
......@@ -22,7 +22,10 @@
# yes - do check for HDF5 library in standard locations.
# path - complete path to the HDF5 helper script h5cc or h5pcc.
#
# SWIFT modification: HDF5 is required, so only path is described.
# SWIFT modifications: HDF5 is required, so only path is described,
# when the h5cc or h5pcc commands are not available, we check if
# HDF5 can be used anyway and no macros are defined except HAVE_HDF5
# and with_hdf5.
#
# If HDF5 is successfully found, this macro calls
#
......@@ -155,11 +158,37 @@ if test "$with_hdf5" = "yes"; then
AC_MSG_CHECKING([Using provided HDF5 C wrapper])
AC_MSG_RESULT([$H5CC])
fi
AC_MSG_CHECKING([for HDF5 libraries])
if test ! -f "$H5CC" || test ! -x "$H5CC"; then
AC_MSG_RESULT([no])
AC_MSG_WARN(m4_case(m4_normalize([$1]),
[serial], [
dnl Check if we already have HDF5 for C.
AC_CHECK_HEADER([hdf5.h], [ac_cv_hhdf5_h=yes], [ac_cv_hhdf5_h=no], [AC_INCLUDES_DEFAULT])
AC_CHECK_LIB([hdf5], [H5Fcreate], [ac_cv_libhdf5=yes],
[ac_cv_libhdf5=no])
if test "$ac_cv_hhdf5_h" = "yes" && test "$ac_cv_libhdf5" = "yes" ; then
dnl Can compile and link, so we have a HDF5, just don't know which version.
AC_MSG_CHECKING([for HDF5 libraries])
AC_MSG_RESULT([yes])
with_hdf5="yes"
HDF5_VERSION="unknown"
HDF5_LIBS="-lhdf5"
AC_SUBST([HDF5_VERSION])
AC_SUBST([HDF5_CC])
AC_SUBST([HDF5_CFLAGS])
AC_SUBST([HDF5_CPPFLAGS])
AC_SUBST([HDF5_LDFLAGS])
AC_SUBST([HDF5_LIBS])
AC_SUBST([HDF5_FC])
AC_SUBST([HDF5_FFLAGS])
AC_SUBST([HDF5_FLIBS])
AC_DEFINE([HAVE_HDF5], [1], [Defined if you have HDF5 support])
else
dnl Time to give up.
AC_MSG_CHECKING([for HDF5 libraries])
AC_MSG_RESULT([no])
AC_MSG_WARN(m4_case(m4_normalize([$1]),
[serial], [
Unable to locate serial HDF5 compilation helper script 'h5cc'.
Please specify --with-hdf5=<LOCATION> as the full path to h5cc.
HDF5 support is being disabled.
......@@ -172,9 +201,12 @@ Unable to locate HDF5 compilation helper scripts 'h5cc' or 'h5pcc'.
Please specify --with-hdf5=<LOCATION> as the full path to h5cc or h5pcc.
HDF5 support is being disabled.
]))
with_hdf5="no"
with_hdf5_fortran="no"
with_hdf5="no"
with_hdf5_fortran="no"
fi
else
AC_MSG_CHECKING([for HDF5 libraries])
dnl Get the h5cc output
HDF5_SHOW=$(eval $H5CC -show)
......
......@@ -46,7 +46,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
vector_power.h
vector_power.h hydro_space.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
......@@ -56,7 +56,8 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
physical_constants.c potential.c hydro_properties.c \
runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \
statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
part_type.c xmf.c gravity_properties.c gravity.c
part_type.c xmf.c gravity_properties.c gravity.c \
hydro_space.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
......
......@@ -54,6 +54,23 @@
//#define GIZMO_FIX_PARTICLES
//#define GIZMO_TOTAL_ENERGY
/* Types of gradients to use for SHADOWFAX_SPH */
/* If no option is chosen, no gradients are used (first order scheme) */
#define SHADOWFAX_GRADIENTS
/* SHADOWFAX_SPH slope limiters */
#define SHADOWFAX_SLOPE_LIMITER_PER_FACE
#define SHADOWFAX_SLOPE_LIMITER_CELL_WIDE
/* Options to control SHADOWFAX_SPH */
/* This option disables cell movement */
//#define SHADOWFAX_FIX_CELLS
/* This option enables cell steering, i.e. trying to keep the cells regular by
adding a correction to the cell velocities.*/
#define SHADOWFAX_STEER_CELL_MOTION
/* This option evolves the total energy instead of the thermal energy */
//#define SHADOWFAX_TOTAL_ENERGY
/* Source terms */
#define SOURCETERMS_NONE
//#define SOURCETERMS_SN_FEEDBACK
......
......@@ -49,6 +49,8 @@
#include "./hydro/Default/hydro_debug.h"
#elif defined(GIZMO_SPH)
#include "./hydro/Gizmo/hydro_debug.h"
#elif defined(SHADOWFAX_SPH)
#include "./hydro/Shadowswift/hydro_debug.h"
#else
#error "Invalid choice of SPH variant"
#endif
......
......@@ -47,6 +47,11 @@
#include "./hydro/Gizmo/hydro.h"
#include "./hydro/Gizmo/hydro_iact.h"
#define SPH_IMPLEMENTATION "GIZMO (Hopkins 2015)"
#elif defined(SHADOWFAX_SPH)
#include "./hydro/Shadowswift/hydro.h"
#include "./hydro/Shadowswift/hydro_iact.h"
#define SPH_IMPLEMENTATION \
"Shadowfax moving mesh (Vandenbroucke and De Rijcke 2016)"
#else
#error "Invalid choice of SPH variant"
#endif
......
......@@ -22,6 +22,7 @@
#include "adiabatic_index.h"
#include "approx_math.h"
#include "equation_of_state.h"
#include "hydro_space.h"
#include "minmax.h"
#include <float.h>
......@@ -165,9 +166,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
* the variaous density tasks
*
* @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) {
struct part *restrict p, const struct hydro_space *hs) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
......@@ -400,7 +402,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->u_full = p->u;
hydro_reset_acceleration(p);
hydro_init_part(p);
hydro_init_part(p, NULL);
}
#endif /* SWIFT_DEFAULT_HYDRO_H */
......@@ -36,6 +36,7 @@
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "hydro_space.h"
#include "kernel_hydro.h"
#include "minmax.h"
......@@ -169,9 +170,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
* the variaous density tasks
*
* @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) {
struct part *restrict p, const struct hydro_space *hs) {
p->rho = 0.f;
p->density.wcount = 0.f;
......@@ -456,7 +458,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->entropy_full = p->entropy;
hydro_reset_acceleration(p);
hydro_init_part(p);
hydro_init_part(p, NULL);
}
#endif /* SWIFT_GADGET2_HYDRO_H */
......@@ -23,6 +23,7 @@
#include "approx_math.h"
#include "equation_of_state.h"
#include "hydro_gradients.h"
#include "hydro_space.h"
#include "minmax.h"
#include "riemann.h"
......@@ -145,9 +146,10 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
* Simply makes sure all necessary variables are initialized to zero.
*
* @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* p) {
struct part* p, const struct hydro_space* hs) {
p->density.wcount = 0.0f;
p->density.wcount_dh = 0.0f;
......
......@@ -38,6 +38,7 @@
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "hydro_space.h"
#include "kernel_hydro.h"
#include "minmax.h"
......@@ -183,9 +184,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
* 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) {
struct part *restrict p, const struct hydro_space *hs) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
......@@ -429,7 +431,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->u_full = p->u;
hydro_reset_acceleration(p);
hydro_init_part(p);
hydro_init_part(p, NULL);
}
#endif /* SWIFT_MINIMAL_HYDRO_H */
......@@ -36,6 +36,7 @@
#include "dimension.h"
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "hydro_space.h"
#include "kernel_hydro.h"
#include "minmax.h"
......@@ -169,9 +170,10 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
* the variaous density tasks
*
* @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) {
struct part *restrict p, const struct hydro_space *hs) {
p->rho = 0.f;
p->rho_bar = 0.f;
......@@ -474,7 +476,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->v_full[2] = p->v[2];
hydro_reset_acceleration(p);
hydro_init_part(p);
hydro_init_part(p, NULL);
}
#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_H */
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com).
*
* 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/>.
*
******************************************************************************/
#include <float.h>
#include "adiabatic_index.h"
#include "approx_math.h"
#include "equation_of_state.h"
#include "hydro_gradients.h"
#include "hydro_space.h"
#include "voronoi_algorithm.h"
/**
* @brief Computes the hydro time-step of a given particle
*
* @param p Pointer to the particle data.
* @param xp Pointer to the extended particle data.
* @param hydro_properties Pointer to the hydro parameters.
*/
__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 float CFL_condition = hydro_properties->CFL_condition;
float R = get_radius_dimension_sphere(p->cell.volume);
if (p->timestepvars.vmax == 0.) {
/* vmax can be zero in vacuum. We force the time step to become the maximal
time step in this case */
return FLT_MAX;
} else {
return CFL_condition * R / fabsf(p->timestepvars.vmax);
}
}
/**
* @brief Does some extra hydro operations once the actual physical time step
* for the particle is known.
*
* We use this to store the physical time step, since it is used for the flux
* exchange during the force loop.
*
* We also set the active flag of the particle to inactive. It will be set to
* active in hydro_init_part, which is called the next time the particle becomes
* active.
*
* @param p The particle to act upon.
* @param dt Physical time step of the particle during the next step.
*/
__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
struct part* p, float dt) {
p->force.dt = dt;
p->force.active = 0;
}
/**
* @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.
*
* In this case, we copy the particle velocities into the corresponding
* primitive variable field. We do this because the particle velocities in GIZMO
* can be independent of the actual fluid velocity. The latter is stored as a
* primitive variable and integrated using the linear momentum, a conserved
* variable.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part* p, struct xpart* xp) {
const float mass = p->conserved.mass;
p->primitives.v[0] = p->v[0];
p->primitives.v[1] = p->v[1];
p->primitives.v[2] = p->v[2];
p->conserved.momentum[0] = mass * p->primitives.v[0];
p->conserved.momentum[1] = mass * p->primitives.v[1];
p->conserved.momentum[2] = mass * p->primitives.v[2];
#ifdef EOS_ISOTHERMAL_GAS
p->conserved.energy = mass * const_isothermal_internal_energy;
#else
p->conserved.energy *= mass;
#endif
#ifdef SHADOWFAX_TOTAL_ENERGY
p->conserved.energy += 0.5f * (p->conserved.momentum[0] * p->primitives.v[0] +
p->conserved.momentum[1] * p->primitives.v[1] +
p->conserved.momentum[2] * p->primitives.v[2]);
#endif
#if defined(SHADOWFAX_FIX_CELLS)
p->v[0] = 0.;
p->v[1] = 0.;
p->v[2] = 0.;
#endif
/* set the initial velocity of the cells */
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
}
/**
* @brief Prepares a particle for the volume calculation.
*
* Simply makes sure all necessary variables are initialized to zero.
* Initializes the Voronoi cell.
*
* @param p The particle to act upon
* @param hs #hydro_space containing extra information about the space.
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part* p, const struct hydro_space* hs) {
p->density.wcount = 0.0f;
p->density.wcount_dh = 0.0f;
voronoi_cell_init(&p->cell, p->x, hs->anchor, hs->side);
/* Set the active flag to active. */
p->force.active = 1;
}
/**
* @brief Finishes the volume calculation.
*
* Calls the finalize method on the Voronoi cell, which calculates the volume
* and centroid of the cell. We use the return value of this function to set
* a new value for the smoothing length and possibly force another iteration
* of the volume calculation for this particle. We then use the volume to
* convert conserved variables into primitive variables.
*
* This method also initializes the gradient variables (if gradients are used).
*
* @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_end_density(
struct part* restrict p) {
float volume;
float m, momentum[3], energy;
hydro_gradients_init(p);
float hnew = voronoi_cell_finalize(&p->cell);
/* Enforce hnew as new smoothing length in the iteration
This is annoyingly difficult, as we do not have access to the variables
that govern the loop...
So here's an idea: let's force in some method somewhere that makes sure
r->e->hydro_properties->target_neighbours is 1, and
r->e->hydro_properties->delta_neighbours is 0.
This way, we can accept an iteration by setting p->density.wcount to 1.
To get the right correction for h, we set wcount to something else
(say 0), and then set p->density.wcount_dh to 1/(hnew-h). */
if (hnew < p->h) {
/* Iteration succesful: we accept, but manually set h to a smaller value
for the next time step */
p->density.wcount = 1.0f;
p->h = 1.1f * hnew;
} else {
/* Iteration not succesful: we force h to become 1.1*hnew */
p->density.wcount = 0.0f;
p->density.wcount_dh = 1.0f / (1.1f * hnew - p->h);
return;
}
volume = p->cell.volume;
#ifdef SWIFT_DEBUG_CHECKS
/* the last condition checks for NaN: a NaN value always evaluates to false,
even when checked against itself */
if (volume == 0. || volume == INFINITY || volume != volume) {
error("Invalid value for cell volume (%g)!", volume);
}
#endif
/* compute primitive variables */
/* eqns (3)-(5) */
m = p->conserved.mass;
if (m > 0.) {
momentum[0] = p->conserved.momentum[0];
momentum[1] = p->conserved.momentum[1];
momentum[2] = p->conserved.momentum[2];
p->primitives.rho = m / volume;
p->primitives.v[0] = momentum[0] / m;
p->primitives.v[1] = momentum[1] / m;
p->primitives.v[2] = momentum[2] / m;
energy = p->conserved.energy;
#ifdef SHADOWFAX_TOTAL_ENERGY
energy -= 0.5f * (momentum[0] * p->primitives.v[0] +
momentum[1] * p->primitives.v[1] +
momentum[2] * p->primitives.v[2]);
#endif
energy /= m;
p->primitives.P =
gas_pressure_from_internal_energy(p->primitives.rho, energy);
} else {
p->primitives.rho = 0.;
p->primitives.v[0] = 0.;
p->primitives.v[1] = 0.;
p->primitives.v[2] = 0.;
p->primitives.P = 0.;
}
#ifdef SWIFT_DEBUG_CHECKS
if (p->primitives.rho < 0.) {
error("Negative density!");
}
if (p->primitives.P < 0.) {
error("Negative pressure!");
}
#endif
}
/**
* @brief Prepare a particle for the gradient calculation.
*
* The name of this method is confusing, as this method is really called after
* the density loop and before the gradient loop.
*
* We use it to set the physical timestep for the particle and to copy the
* actual velocities, which we need to boost our interfaces during the flux
* calculation. We also initialize the variables used for the time step
* calculation.
*
* @param p The particle to act upon.
* @param xp The extended particle data to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* restrict p, struct xpart* restrict xp) {
/* Initialize time step criterion variables */
p->timestepvars.vmax = 0.0f;
/* Set the actual velocity of the particle */
p->force.v_full[0] = xp->v_full[0];
p->force.v_full[1] = xp->v_full[1];
p->force.v_full[2] = xp->v_full[2];
}
/**
* @brief Finishes the gradient calculation.
*
* Just a wrapper around hydro_gradients_finalize, which can be an empty method,
* in which case no gradients are used.
*
* @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_end_gradient(
struct part* p) {
hydro_gradients_finalize(p);
}
/**
* @brief Reset acceleration fields of a particle
*
* This is actually not necessary for Shadowswift, since we just set the
* accelerations after the flux calculation.
*
* @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
struct part* p) {
/* Reset the acceleration. */
p->a_hydro[0] = 0.0f;
p->a_hydro[1] = 0.0f;
p->a_hydro[2] = 0.0f;
/* Reset the time derivatives. */
p->force.h_dt = 0.0f;