Commit e9295f04 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'shadowswift_new' into 'master'

1D and 2D moving mesh algorithm

I have added a new hydro scheme, `SHADOWFAX_SPH`, that uses a finite volume method on a moving mesh. The scheme is very similar to `GIZMO_SPH`, except that during the density loop an actual Voronoi cell is constructed for every particle, which is then used to get a volume to update the primitive quantities, and faces through which fluxes are exchanged. There are also some small changes in the way gradients are computed, and in the way the velocity of the particles is set.

The Voronoi cell construction is completely different for 1D, 2D, and 3D, and currently only 1D and 2D have a stable implementation. 3D should compile, but does not produce reliable results (yet). However, I think it would be good to merge this in already, to prevent it from drifting too far apart from `master`.

The current implementation is hydro only, but it should be very straightforward to add gravity, similarly to the way it is done in `GIZMO_SPH`.

Below the results for the 2D Sod shock and Sedov blast test (the noise on the density is caused by the fact that the initial conditions contain mass and not density). The 2D Noh test crashes because of time stepping issues; I think we really need a time step limiter to make it work properly.

![SodShock](/uploads/641796ae26cf045be5ce824bb698b77a/SodShock.png)

![Sedov](/uploads/96be4bbeb7635c329b78f89e0e79239c/Sedov.png)

See merge request !321
parents 3e0294f6 cf5cb6ad
......@@ -73,6 +73,9 @@ tests/testRiemannExact
tests/testRiemannTRRS
tests/testRiemannHLLC
tests/testMatrixInversion
tests/testVoronoi1D
tests/testVoronoi2D
tests/testVoronoi3D
tests/testDump
tests/testLogger
tests/benchmarkInteractions
......
......@@ -587,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"]
......@@ -608,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])
......
......@@ -45,7 +45,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
physical_constants.h physical_constants_cgs.h potential.h version.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
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
hydro_space.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
......@@ -55,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 \
......
......@@ -57,6 +57,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;
}
/**
* @brief Sets the values to be predicted in the drifts to their values at a
* kick time
*
* @param p The particle.
* @param xp The extended data of this particle.
*/
__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
struct part* restrict p, const struct xpart* restrict xp) {}
/**
* @brief Converts the hydrodynamic variables from the initial condition file to
* conserved variables that can be used during the integration
*
* Requires the volume to be known.
*
* The initial condition file contains a mixture of primitive and conserved
* variables. Mass is a conserved variable, and we just copy the particle
* mass into the corresponding conserved quantity. We need the volume to
* also derive a density, which is then used to convert the internal energy
* to a pressure. However, we do not actually use these variables anymore.
* We do need to initialize the linear momentum, based on the mass and the
* velocity of the particle.
*
* @param p The particle to act upon.
* @param xp The extended particle data to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part* p, struct xpart* xp) {}
/**
* @brief Extra operations to be done during the drift
*
* Not used for Shadowswift.
*
* @param p Particle to act upon.
* @param xp The extended particle data to act upon.
* @param dt The drift time-step.
*/
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part* p, struct xpart* xp, float dt) {}
/**
* @brief Set the particle acceleration after the flux loop.
*
* @param p Particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part* p) {}
/**
* @brief Extra operations done during the kick
*
* Not used for Shadowswift.
*
* @param p Particle to act upon.
* @param xp Extended particle data to act upon.
* @param dt Physical time step.
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt) {
float vcell[3];
/* Update the conserved variables. We do this here and not in the kick,
since we need the updated variables below. */
p->conserved.mass += p->conserved.flux.mass;
p->conserved.momentum[0] += p->conserved.flux.momentum[0];
p->conserved.momentum[1] += p->conserved.flux.momentum[1];
p->conserved.momentum[2] += p->conserved.flux.momentum[2];
p->conserved.energy += p->conserved.flux.energy;
#ifdef EOS_ISOTHERMAL_GAS
/* reset the thermal energy */
p->conserved.energy = p->conserved.mass * const_isothermal_internal_energy;
#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
#endif
/* reset fluxes */
/* we can only do this here, since we need to keep the fluxes for inactive
particles */
p->conserved.flux.mass = 0.0f;