From 00128697949fe59e68233931b36335beefb04b59 Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <bert.vandenbroucke@ugent.be>
Date: Wed, 17 Aug 2016 21:08:37 +0100
Subject: [PATCH] Fully documented 1D moving mesh implementation.
---
src/engine.c | 2 +-
src/hydro/Shadowswift/hydro.h | 30 ++++----
src/hydro/Shadowswift/hydro_debug.h | 2 +-
.../Shadowswift/hydro_gradients_shadowfax.h | 14 ++++
src/hydro/Shadowswift/hydro_iact.h | 33 +++-----
src/hydro/Shadowswift/hydro_io.h | 23 +++++-
src/hydro/Shadowswift/hydro_part.h | 4 +-
.../Shadowswift/hydro_slope_limiters_cell.h | 11 +++
src/hydro/Shadowswift/voronoi1d_algorithm.h | 76 +++++++++++++++++--
src/hydro/Shadowswift/voronoi1d_cell.h | 1 +
10 files changed, 145 insertions(+), 51 deletions(-)
diff --git a/src/engine.c b/src/engine.c
index 84e5655dd5..29d35e087b 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2724,7 +2724,7 @@ void engine_step(struct engine *e) {
mask |= 1 << task_type_sub_pair;
mask |= 1 << task_type_ghost;
mask |= 1 << task_type_extra_ghost; /* Adding unnecessary things to the mask
- does not harm */
+ does not harm */
submask |= 1 << task_subtype_density;
submask |= 1 << task_subtype_gradient;
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 258115f48a..b70ed690d3 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -72,6 +72,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
* @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
*/
@@ -88,19 +89,13 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
/**
* @brief Finishes the volume calculation.
*
- * Multiplies the density and number of neighbours by the appropiate constants
- * and adds the self-contribution term. Calculates the volume and uses it to
- * update the primitive variables (based on the conserved variables). The latter
- * should only be done for active particles. This is okay, since this method is
- * only called for active particles.
+ * 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.
*
- * Multiplies the components of the matrix E with the appropriate constants and
- * inverts it. Initializes the variables used during the gradient loop. This
- * cannot be done in hydro_prepare_force, since that method is called for all
- * particles, and not just the active ones. If we would initialize the
- * variables there, gradients for passive particles would be zero, while we
- * actually use the old gradients in the flux calculation between active and
- * passive particles.
+ * This method also initializes the gradient variables (if gradients are used).
*
* @param p The particle to act upon.
* @param The current physical time.
@@ -200,8 +195,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
/**
* @brief Reset acceleration fields of a particle
*
- * This is actually not necessary for GIZMO, since we just set the accelerations
- * after the flux calculation.
+ * This is actually not necessary for Shadowswift, since we just set the
+ * accelerations after the flux calculation.
*
* @param p The particle to act upon.
*/
@@ -256,7 +251,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
/**
* @brief Extra operations to be done during the drift
*
- * Not used for GIZMO.
+ * Not used for Shadowswift.
*
* @param p Particle to act upon.
* @param xp The extended particle data to act upon.
@@ -272,7 +267,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
*
* We use the new conserved variables to calculate the new velocity of the
* particle, and use that to derive the change of the velocity over the particle
- * time step.
+ * time step. We also add a correction to the velocity which steers the
+ * generator towards the centroid of the cell.
*
* If the particle time step is zero, we set the accelerations to zero. This
* should only happen at the start of the simulation.
@@ -356,7 +352,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
/**
* @brief Extra operations done during the kick
*
- * Not used for GIZMO.
+ * Not used for Shadowswift.
*
* @param p Particle to act upon.
* @param xp Extended particle data to act upon.
diff --git a/src/hydro/Shadowswift/hydro_debug.h b/src/hydro/Shadowswift/hydro_debug.h
index fa6c25dd30..d867152a99 100644
--- a/src/hydro/Shadowswift/hydro_debug.h
+++ b/src/hydro/Shadowswift/hydro_debug.h
@@ -1,6 +1,6 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Coypright (c) 2016 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
diff --git a/src/hydro/Shadowswift/hydro_gradients_shadowfax.h b/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
index 4a7df6560f..9977c36d94 100644
--- a/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
+++ b/src/hydro/Shadowswift/hydro_gradients_shadowfax.h
@@ -50,6 +50,20 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_init(
hydro_slope_limit_cell_init(p);
}
+/**
+ * @brief Add the gradient estimate for a single quantity due to a particle pair
+ * to the total gradient for that quantity
+ *
+ * This corresponds to one term of equation (21) in Springel (2010).
+ *
+ * @param qL Value of the quantity on the left.
+ * @param qR Value of the quantity on the right.
+ * @param cLR Vector pointing from the midpoint of the particle pair to the
+ * geometrical centroid of the face in between the particles.
+ * @param xLR Vector pointing from the right particle to the left particle.
+ * @param A Surface area of the face in between the particles.
+ * @param grad Current value of the gradient for the quantity (is updated).
+ */
__attribute__((always_inline)) INLINE void hydro_gradients_single_quantity(
float qL, float qR, float *cLR, float *xLR, float rLR, float A,
float *grad) {
diff --git a/src/hydro/Shadowswift/hydro_iact.h b/src/hydro/Shadowswift/hydro_iact.h
index 40d5da9595..1c1ec9e1ad 100644
--- a/src/hydro/Shadowswift/hydro_iact.h
+++ b/src/hydro/Shadowswift/hydro_iact.h
@@ -1,8 +1,6 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
- * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- * 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ * 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
@@ -25,14 +23,9 @@
#include "voronoi_algorithm.h"
/**
- * @brief Calculate the volume interaction between particle i and particle j
+ * @brief Calculate the Voronoi cell by interacting particle pi and pj
*
- * The volume is in essence the same as the weighted number of neighbours in a
- * classical SPH density calculation.
- *
- * We also calculate the components of the matrix E, which is used for second
- * order accurate gradient calculations and for the calculation of the interface
- * surface areas.
+ * This method wraps around voronoi_cell_interact().
*
* @param r2 Squared distance between particle i and particle j.
* @param dx Distance vector between the particles (dx = pi->x - pj->x).
@@ -54,15 +47,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
}
/**
- * @brief Calculate the volume interaction between particle i and particle j:
- * non-symmetric version
- *
- * The volume is in essence the same as the weighted number of neighbours in a
- * classical SPH density calculation.
+ * @brief Calculate the Voronoi cell by interacting particle pi with pj
*
- * We also calculate the components of the matrix E, which is used for second
- * order accurate gradient calculations and for the calculation of the interface
- * surface areas.
+ * This method wraps around voronoi_cell_interact().
*
* @param r2 Squared distance between particle i and particle j.
* @param dx Distance vector between the particles (dx = pi->x - pj->x).
@@ -125,11 +112,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
* active), both runner_iact_force and runner_iact_nonsym_force call this
* method, with an appropriate mode.
*
- * This method calculates the surface area of the interface between particle i
- * and particle j, as well as the interface position and velocity. These are
- * then used to reconstruct and predict the primitive variables, which are then
- * fed to a Riemann solver that calculates a flux. This flux is used to update
- * the conserved variables of particle i or both particles.
+ * This method retrieves the oriented surface area and face midpoint for the
+ * Voronoi face between pi and pj (if it exists). It uses the midpoint position
+ * to reconstruct the primitive quantities (if gradients are used) at the face
+ * and then uses the face quantities to estimate a flux through the face using
+ * a Riemann solver.
*
* This method also calculates the maximal velocity used to calculate the time
* step.
diff --git a/src/hydro/Shadowswift/hydro_io.h b/src/hydro/Shadowswift/hydro_io.h
index b3e01bf120..ffde309562 100644
--- a/src/hydro/Shadowswift/hydro_io.h
+++ b/src/hydro/Shadowswift/hydro_io.h
@@ -1,6 +1,6 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Coypright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ * 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
@@ -55,14 +55,35 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, primitives.rho);
}
+/**
+ * @brief Get the internal energy of a particle
+ *
+ * @param e #engine.
+ * @param p Particle.
+ * @return Internal energy of the particle
+ */
float convert_u(struct engine* e, struct part* p) {
return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
}
+/**
+ * @brief Get the entropic function of a particle
+ *
+ * @param e #engine.
+ * @param p Particle.
+ * @return Entropic function of the particle
+ */
float convert_A(struct engine* e, struct part* p) {
return p->primitives.P / pow_gamma(p->primitives.rho);
}
+/**
+ * @brief Get the total energy of a particle
+ *
+ * @param e #engine.
+ * @param p Particle.
+ * @return Total energy of the particle
+ */
float convert_Etot(struct engine* e, struct part* p) {
float momentum2;
diff --git a/src/hydro/Shadowswift/hydro_part.h b/src/hydro/Shadowswift/hydro_part.h
index f0a360c3e4..46f9ec4509 100644
--- a/src/hydro/Shadowswift/hydro_part.h
+++ b/src/hydro/Shadowswift/hydro_part.h
@@ -1,8 +1,8 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
* Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- * Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ * 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
diff --git a/src/hydro/Shadowswift/hydro_slope_limiters_cell.h b/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
index 2e62bce4fc..a7b3f7511f 100644
--- a/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
+++ b/src/hydro/Shadowswift/hydro_slope_limiters_cell.h
@@ -80,6 +80,17 @@ hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
pi->primitives.limiter.maxr = fmax(r, pi->primitives.limiter.maxr);
}
+/**
+ * @brief Apply the cell wide slope limiter to the gradient of a single quantity
+ *
+ * This corresponds to equation (B2) in Hopkins (2015).
+ *
+ * @param grad Gradient to slope limit
+ * @param qval Value of the quantity at the cell generator
+ * @param qmin Minimal value of the quantity among all cell neighbours
+ * @param qmax Maximal value of the quantity among all cell neighbours
+ * @param maxr Maximal distance between the generator and all of its neighbours
+ */
__attribute__((always_inline)) INLINE static void
hydro_slope_limit_cell_quantity(float* grad, float qval, float qmin, float qmax,
float maxr) {
diff --git a/src/hydro/Shadowswift/voronoi1d_algorithm.h b/src/hydro/Shadowswift/voronoi1d_algorithm.h
index c58b2474fa..222c97b617 100644
--- a/src/hydro/Shadowswift/voronoi1d_algorithm.h
+++ b/src/hydro/Shadowswift/voronoi1d_algorithm.h
@@ -27,6 +27,16 @@
#include "inline.h"
#include "voronoi1d_cell.h"
+/**
+ * @brief Initialize a 1D Voronoi cell
+ *
+ * Sets the positions of left and right neighbours to very large values, the
+ * generator position to the given particle position, and all other quantities
+ * to zero.
+ *
+ * @param cell 1D Voronoi cell to initialize.
+ * @param x Position of the generator of the cell.
+ */
__attribute__((always_inline)) INLINE void voronoi_cell_init(
struct voronoi_cell *cell, double *x) {
cell->x = x[0];
@@ -38,6 +48,19 @@ __attribute__((always_inline)) INLINE void voronoi_cell_init(
cell->centroid = 0.0f;
}
+/**
+ * @brief Interact a 1D Voronoi cell with a particle with given relative
+ * position and ID
+ *
+ * This method checks if the given relative position is closer to the cell
+ * generator than the current left or right neighbour and updates neighbours
+ * accordingly.
+ *
+ * @param cell 1D Voronoi cell.
+ * @param dx Relative position of the interacting generator w.r.t. the cell
+ * generator (in fact: dx = generator - neighbour).
+ * @param id ID of the interacting neighbour.
+ */
__attribute__((always_inline)) INLINE void voronoi_cell_interact(
struct voronoi_cell *cell, float *dx, unsigned long long id) {
@@ -61,6 +84,22 @@ __attribute__((always_inline)) INLINE void voronoi_cell_interact(
}
}
+/**
+ * @brief Finalize a 1D Voronoi cell
+ *
+ * Calculates the relative positions of the midpoints of the faces (which in
+ * this case are just the midpoints of the segments connecting the generator
+ * with the two neighbours) w.r.t. the generator, and the cell volume (length)
+ * and centroid (midpoint of the segment connecting the midpoints of the faces).
+ * This function returns the maximal radius at which a particle could still
+ * change the structure of the cell, i.e. twice the largest distance between
+ * the cell generator and one of its faces. If the cell has been interacted with
+ * all neighbours within this radius, we know for sure that the cell is
+ * complete.
+ *
+ * @param cell 1D Voronoi cell.
+ * @return Maximal radius that could still change the structure of the cell.
+ */
__attribute__((always_inline)) INLINE float voronoi_cell_finalize(
struct voronoi_cell *cell) {
@@ -77,6 +116,28 @@ __attribute__((always_inline)) INLINE float voronoi_cell_finalize(
return max_radius;
}
+/**
+ * @brief Get the oriented surface area and midpoint of the face between a
+ * 1D Voronoi cell and the given neighbour
+ *
+ * This function also checks if the given neighbour is in fact a neighbour of
+ * this cell. Since we perform gradient and flux calculations for all neighbour
+ * pairs within the smoothing length, which assumes the cell to be spherical,
+ * it can happen that this is not the case. It is the responsibility of the
+ * routine that calls this function to check for a zero return value and
+ * deal with it appropriately.
+ *
+ * For this specific case, we simply check if the neighbour is the left or
+ * right neighbour and set the surface area to the unit vector pointing either
+ * to the right or to the left. The midpoint is set to the relative position
+ * vector of the appropriate face.
+ *
+ * @param cell 1D Voronoi cell.
+ * @param ngb ID of a particle that is possibly a neighbour of this cell.
+ * @param A Array to store the oriented surface area in.
+ * @param midpoint Array to store the relative position of the face in.
+ * @return 0 if the given neighbour is not a neighbour, 1 otherwise.
+ */
__attribute__((always_inline)) INLINE int voronoi_get_face(
struct voronoi_cell *cell, unsigned long long ngb, float *A,
float *midpoint) {
@@ -106,6 +167,15 @@ __attribute__((always_inline)) INLINE int voronoi_get_face(
return 1;
}
+/**
+ * @brief Get the centroid of a 1D Voronoi cell
+ *
+ * We store only the relevant coordinate of the centroid, but need to return
+ * a 3D vector.
+ *
+ * @param cell 1D Voronoi cell.
+ * @param centroid Array to store the centroid in.
+ */
__attribute__((always_inline)) INLINE void voronoi_get_centroid(
struct voronoi_cell *cell, float *centroid) {
@@ -114,10 +184,4 @@ __attribute__((always_inline)) INLINE void voronoi_get_centroid(
centroid[2] = 0.0f;
}
-__attribute__((always_inline)) INLINE int voronoi_is_neighbour(
- struct voronoi_cell *cell, unsigned long long ngb) {
-
- return (ngb == cell->idL || ngb == cell->idR);
-}
-
#endif // SWIFT_VORONOI1D_ALGORITHM_H
diff --git a/src/hydro/Shadowswift/voronoi1d_cell.h b/src/hydro/Shadowswift/voronoi1d_cell.h
index 2dc9bcf4e9..b286609593 100644
--- a/src/hydro/Shadowswift/voronoi1d_cell.h
+++ b/src/hydro/Shadowswift/voronoi1d_cell.h
@@ -20,6 +20,7 @@
#ifndef SWIFT_VORONOI1D_CELL_H
#define SWIFT_VORONOI1D_CELL_H
+/* 1D Voronoi cell */
struct voronoi_cell {
/* The position of the generator of the cell. */
--
GitLab