diff --git a/.gitignore b/.gitignore
index 672a71c3ed03ec73df5dbc917572c8f87a0967d1..8c3ede8f3125f5024c1fe01a405024d1cf5a7f19 100644
--- a/.gitignore
+++ b/.gitignore
@@ -119,6 +119,7 @@ tests/testPotentialPair
tests/testEOS
tests/testEOS*.txt
tests/testEOS*.png
+tests/testUtilities
theory/latex/swift.pdf
theory/SPH/Kernels/kernels.pdf
diff --git a/configure.ac b/configure.ac
index 0640d98aaaa1494da52def5a4118a9204436e608..56018272debbb9e9ccb2f425eff45b97c358aec9 100644
--- a/configure.ac
+++ b/configure.ac
@@ -855,7 +855,7 @@ fi
# Hydro scheme.
AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>],
- [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, default, gizmo, shadowfax debug default: gadget2@:>@]
+ [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, default, gizmo-mfv, gizmo-mfm, shadowfax, minimal-multi-mat, debug default: gadget2@:>@]
)],
[with_hydro="$withval"],
[with_hydro="gadget2"]
@@ -885,8 +885,11 @@ case "$with_hydro" in
default)
AC_DEFINE([DEFAULT_SPH], [1], [Default SPH])
;;
- gizmo)
- AC_DEFINE([GIZMO_SPH], [1], [GIZMO SPH])
+ gizmo-mfv)
+ AC_DEFINE([GIZMO_MFV_SPH], [1], [GIZMO MFV SPH])
+ ;;
+ gizmo-mfm)
+ AC_DEFINE([GIZMO_MFM_SPH], [1], [GIZMO MFM SPH])
;;
shadowfax)
AC_DEFINE([SHADOWFAX_SPH], [1], [Shadowfax SPH])
diff --git a/examples/MoonFormingImpact/moon_forming_impact.yml b/examples/MoonFormingImpact/moon_forming_impact.yml
index 7d726bc02ba4fb6c8dd02f74907ffc48c3ed9431..323adf7f3ac73f41b45b50eaa76a95033dca35d7 100644
--- a/examples/MoonFormingImpact/moon_forming_impact.yml
+++ b/examples/MoonFormingImpact/moon_forming_impact.yml
@@ -42,3 +42,7 @@ Gravity:
InitialConditions:
# The initial conditions file to read
file_name: moon_forming_impact.hdf5
+
+# Parameters related to the equation of state
+EoS:
+ planetary_use_Til: 1 # Whether to prepare the Tillotson EOS
diff --git a/examples/UranusImpact/uranus_impact.yml b/examples/UranusImpact/uranus_impact.yml
index aae9f66847a9f9b55984ea5d2ea1c79099c01e95..fabddca00f80fcdd79ff6114ff0544cd251046f4 100644
--- a/examples/UranusImpact/uranus_impact.yml
+++ b/examples/UranusImpact/uranus_impact.yml
@@ -41,3 +41,11 @@ Gravity:
# Parameters related to the initial conditions
InitialConditions:
file_name: uranus_impact.hdf5 # The initial conditions file to read
+
+# Parameters related to the equation of state
+EoS:
+ planetary_use_HM80: 1 # Whether to prepare the Hubbard & MacFarlane (1980) EOS
+ # Table file paths
+ planetary_HM80_HHe_table_file: /gpfs/data/dc-kege1/gihr_data/P_rho_u_HHe.txt
+ planetary_HM80_ice_table_file: /gpfs/data/dc-kege1/gihr_data/P_rho_u_ice.txt
+ planetary_HM80_rock_table_file: /gpfs/data/dc-kege1/gihr_data/P_rho_u_roc.txt
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index c9bc9811803051b274cd09cdcf6b58d1a65a0d7d..791db2758290e400d4ed9ffe5b8e0d4303057874 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -131,7 +131,16 @@ DomainDecomposition:
# Parameters related to the equation of state ------------------------------------------
EoS:
- isothermal_internal_energy: 20.26784 # Thermal energy per unit mass for the case of isothermal equation of state (in internal units).
+ isothermal_internal_energy: 20.26784 # Thermal energy per unit mass for the case of isothermal equation of state (in internal units).
+
+ planetary_use_Til: 1 # (Optional) Whether to prepare the Tillotson EOS
+ planetary_use_HM80: 0 # (Optional) Whether to prepare the Hubbard & MacFarlane (1980) EOS
+ planetary_use_ANEOS: 0 # (Optional) Whether to prepare the ANEOS EOS
+ planetary_use_SESAME: 0 # (Optional) Whether to prepare the SESAME EOS
+ # (Optional) Table file paths
+ planetary_HM80_HHe_table_file: HM80_HHe.txt
+ planetary_HM80_ice_table_file: HM80_ice.txt
+ planetary_HM80_rock_table_file: HM80_rock.txt
# Parameters related to external potentials --------------------------------------------
diff --git a/src/Makefile.am b/src/Makefile.am
index 677a4e9524a733de92d176992a67f7a18653c82d..374d2c4569efe0e07ca1b2558b679b33564da432 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -47,7 +47,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.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 \
gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h \
- chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h
+ chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
@@ -81,17 +81,28 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \
hydro/PressureEntropy/hydro.h hydro/PressureEntropy/hydro_iact.h hydro/PressureEntropy/hydro_io.h \
hydro/PressureEntropy/hydro_debug.h hydro/PressureEntropy/hydro_part.h \
- hydro/Gizmo/hydro.h hydro/Gizmo/hydro_iact.h \
- hydro/Gizmo/hydro_io.h hydro/Gizmo/hydro_debug.h \
- hydro/Gizmo/hydro_part.h \
- hydro/Gizmo/hydro_gradients_gizmo.h \
- hydro/Gizmo/hydro_gradients.h \
- hydro/Gizmo/hydro_gradients_sph.h \
- hydro/Gizmo/hydro_slope_limiters_cell.h \
- hydro/Gizmo/hydro_slope_limiters_face.h \
- hydro/Gizmo/hydro_slope_limiters.h \
- hydro/Gizmo/hydro_unphysical.h \
- hydro/Gizmo/hydro_velocities.h \
+ hydro/GizmoMFV/hydro.h hydro/GizmoMFV/hydro_iact.h \
+ hydro/GizmoMFV/hydro_io.h hydro/GizmoMFV/hydro_debug.h \
+ hydro/GizmoMFV/hydro_part.h \
+ hydro/GizmoMFV/hydro_gradients_gizmo.h \
+ hydro/GizmoMFV/hydro_gradients.h \
+ hydro/GizmoMFV/hydro_gradients_sph.h \
+ hydro/GizmoMFV/hydro_slope_limiters_cell.h \
+ hydro/GizmoMFV/hydro_slope_limiters_face.h \
+ hydro/GizmoMFV/hydro_slope_limiters.h \
+ hydro/GizmoMFV/hydro_unphysical.h \
+ hydro/GizmoMFV/hydro_velocities.h \
+ hydro/GizmoMFM/hydro.h hydro/GizmoMFM/hydro_iact.h \
+ hydro/GizmoMFM/hydro_io.h hydro/GizmoMFM/hydro_debug.h \
+ hydro/GizmoMFM/hydro_part.h \
+ hydro/GizmoMFM/hydro_gradients_gizmo.h \
+ hydro/GizmoMFM/hydro_gradients.h \
+ hydro/GizmoMFM/hydro_gradients_sph.h \
+ hydro/GizmoMFM/hydro_slope_limiters_cell.h \
+ hydro/GizmoMFM/hydro_slope_limiters_face.h \
+ hydro/GizmoMFM/hydro_slope_limiters.h \
+ hydro/GizmoMFM/hydro_unphysical.h \
+ hydro/GizmoMFM/hydro_velocities.h \
hydro/Shadowswift/hydro_debug.h \
hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
hydro/Shadowswift/hydro_iact.h \
diff --git a/src/const.h b/src/const.h
index 928835c5ae0ac70d641243b0d183c2c0652cdc8a..6c5b5299c08efb7935b046ecfd0b3d67b7dc4c7a 100644
--- a/src/const.h
+++ b/src/const.h
@@ -53,7 +53,8 @@
/* This option disables particle movement */
//#define GIZMO_FIX_PARTICLES
/* Try to keep cells regular by adding a correction velocity. */
-#define GIZMO_STEER_MOTION
+//#define GIZMO_STEER_MOTION
+/* Use the total energy instead of the thermal energy as conserved variable. */
//#define GIZMO_TOTAL_ENERGY
/* Options to control handling of unphysical values (GIZMO_SPH only). */
diff --git a/src/debug.c b/src/debug.c
index 2ca6a9bf832e2da1860d4f89f032c8a4f7bef86b..93d14952f523be5f1d1fa90484e9e7951f8e3f6e 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -50,8 +50,10 @@
#include "./hydro/PressureEnergy/hydro_debug.h"
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_debug.h"
-#elif defined(GIZMO_SPH)
-#include "./hydro/Gizmo/hydro_debug.h"
+#elif defined(GIZMO_MFV_SPH)
+#include "./hydro/GizmoMFV/hydro_debug.h"
+#elif defined(GIZMO_MFM_SPH)
+#include "./hydro/GizmoMFM/hydro_debug.h"
#elif defined(SHADOWFAX_SPH)
#include "./hydro/Shadowswift/hydro_debug.h"
#elif defined(MINIMAL_MULTI_MAT_SPH)
diff --git a/src/equation_of_state/planetary/equation_of_state.h b/src/equation_of_state/planetary/equation_of_state.h
index fcc1c959152801253d465915d8b28b1552c913b8..d7a7f64f87bb126094163ada16dbd05e48d8cf8d 100644
--- a/src/equation_of_state/planetary/equation_of_state.h
+++ b/src/equation_of_state/planetary/equation_of_state.h
@@ -25,7 +25,7 @@
*
* For any/all of the planetary EOS. Each EOS type's functions are set in its
* own header file: `equation_of_state/planetary/<eos_type>.h`.
- * See `material_id` for the available choices.
+ * See `eos_planetary_material_id` for the available choices.
*
* Not all functions are implemented for all EOS types, so not all can be used
* with all hydro formulations yet.
@@ -95,7 +95,7 @@ enum eos_planetary_material_id {
eos_planetary_id_ANEOS_iron =
eos_planetary_type_ANEOS * eos_planetary_type_factor,
- /*! ANEOS forsterite */
+ /*! MANEOS forsterite */
eos_planetary_id_MANEOS_forsterite =
eos_planetary_type_ANEOS * eos_planetary_type_factor + 1,
@@ -1077,46 +1077,62 @@ __attribute__((always_inline)) INLINE static void eos_init(
struct eos_parameters *e, const struct phys_const *phys_const,
const struct unit_system *us, const struct swift_params *params) {
+ // Table file names
+ char HM80_HHe_table_file[PARSER_MAX_LINE_SIZE];
+ char HM80_ice_table_file[PARSER_MAX_LINE_SIZE];
+ char HM80_rock_table_file[PARSER_MAX_LINE_SIZE];
+
// Set the parameters and material IDs, load tables, etc. for each material
+ // and convert to internal units
// Tillotson
- set_Til_iron(&e->Til_iron, eos_planetary_id_Til_iron);
- set_Til_granite(&e->Til_granite, eos_planetary_id_Til_granite);
- set_Til_water(&e->Til_water, eos_planetary_id_Til_water);
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_Til", 0)) {
+ set_Til_iron(&e->Til_iron, eos_planetary_id_Til_iron);
+ set_Til_granite(&e->Til_granite, eos_planetary_id_Til_granite);
+ set_Til_water(&e->Til_water, eos_planetary_id_Til_water);
+
+ convert_units_Til(&e->Til_iron, us);
+ convert_units_Til(&e->Til_granite, us);
+ convert_units_Til(&e->Til_water, us);
+ }
// Hubbard & MacFarlane (1980)
- set_HM80_HHe(&e->HM80_HHe, eos_planetary_id_HM80_HHe);
- set_HM80_ice(&e->HM80_ice, eos_planetary_id_HM80_ice);
- set_HM80_rock(&e->HM80_rock, eos_planetary_id_HM80_rock);
-
- load_HM80_table(&e->HM80_HHe, HM80_HHe_table_file);
- load_HM80_table(&e->HM80_ice, HM80_ice_table_file);
- load_HM80_table(&e->HM80_rock, HM80_rock_table_file);
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_HM80", 0)) {
+ set_HM80_HHe(&e->HM80_HHe, eos_planetary_id_HM80_HHe);
+ set_HM80_ice(&e->HM80_ice, eos_planetary_id_HM80_ice);
+ set_HM80_rock(&e->HM80_rock, eos_planetary_id_HM80_rock);
+
+ parser_get_param_string(params, "EoS:planetary_HM80_HHe_table_file",
+ HM80_HHe_table_file);
+ parser_get_param_string(params, "EoS:planetary_HM80_ice_table_file",
+ HM80_ice_table_file);
+ parser_get_param_string(params, "EoS:planetary_HM80_rock_table_file",
+ HM80_rock_table_file);
+
+ load_HM80_table(&e->HM80_HHe, HM80_HHe_table_file);
+ load_HM80_table(&e->HM80_ice, HM80_ice_table_file);
+ load_HM80_table(&e->HM80_rock, HM80_rock_table_file);
+
+ convert_units_HM80(&e->HM80_HHe, us);
+ convert_units_HM80(&e->HM80_ice, us);
+ convert_units_HM80(&e->HM80_rock, us);
+ }
// ANEOS
- set_ANEOS_iron(&e->ANEOS_iron, eos_planetary_id_ANEOS_iron);
- set_MANEOS_forsterite(&e->MANEOS_forsterite,
- eos_planetary_id_MANEOS_forsterite);
-
- // SESAME
- set_SESAME_iron(&e->SESAME_iron, eos_planetary_id_SESAME_iron);
-
- // Convert to internal units
- // Tillotson
- convert_units_Til(&e->Til_iron, us);
- convert_units_Til(&e->Til_granite, us);
- convert_units_Til(&e->Til_water, us);
-
- // Hubbard & MacFarlane (1980)
- convert_units_HM80(&e->HM80_HHe, us);
- convert_units_HM80(&e->HM80_ice, us);
- convert_units_HM80(&e->HM80_rock, us);
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_ANEOS", 0)) {
+ set_ANEOS_iron(&e->ANEOS_iron, eos_planetary_id_ANEOS_iron);
+ set_MANEOS_forsterite(&e->MANEOS_forsterite,
+ eos_planetary_id_MANEOS_forsterite);
- // ANEOS
- convert_units_ANEOS(&e->ANEOS_iron, us);
- convert_units_ANEOS(&e->MANEOS_forsterite, us);
+ convert_units_ANEOS(&e->ANEOS_iron, us);
+ convert_units_ANEOS(&e->MANEOS_forsterite, us);
+ }
// SESAME
- convert_units_SESAME(&e->SESAME_iron, us);
+ if (parser_get_opt_param_int(params, "EoS:planetary_use_SESAME", 0)) {
+ set_SESAME_iron(&e->SESAME_iron, eos_planetary_id_SESAME_iron);
+
+ convert_units_SESAME(&e->SESAME_iron, us);
+ }
}
/**
diff --git a/src/equation_of_state/planetary/hm80.h b/src/equation_of_state/planetary/hm80.h
index b0669c129b2b4d642ae27e58f2f2bec221e3df9b..0131bab6c447e5a8898e29e13dc3f8f6e1c897c6 100644
--- a/src/equation_of_state/planetary/hm80.h
+++ b/src/equation_of_state/planetary/hm80.h
@@ -41,19 +41,13 @@
// Hubbard & MacFarlane (1980) parameters
struct HM80_params {
- float **table_P_rho_u;
+ float *table_P_rho_u;
int num_rho, num_u;
float log_rho_min, log_rho_max, log_rho_step, inv_log_rho_step, log_u_min,
log_u_max, log_u_step, inv_log_u_step, bulk_mod;
enum eos_planetary_material_id mat_id;
};
-// Table file names
-/// to be read in from the parameter file instead once finished testing...
-#define HM80_HHe_table_file "/gpfs/data/dc-kege1/gihr_data/P_rho_u_HHe.txt"
-#define HM80_ice_table_file "/gpfs/data/dc-kege1/gihr_data/P_rho_u_ice.txt"
-#define HM80_rock_table_file "/gpfs/data/dc-kege1/gihr_data/P_rho_u_roc.txt"
-
// Parameter values for each material (cgs units)
INLINE static void set_HM80_HHe(struct HM80_params *mat,
enum eos_planetary_material_id mat_id) {
@@ -107,17 +101,15 @@ INLINE static void set_HM80_rock(struct HM80_params *mat,
// Read the table from file
INLINE static void load_HM80_table(struct HM80_params *mat, char *table_file) {
// Allocate table memory
- mat->table_P_rho_u = (float **)malloc(mat->num_rho * sizeof(float *));
- for (int i = 0; i < mat->num_rho; i++) {
- mat->table_P_rho_u[i] = (float *)malloc(mat->num_u * sizeof(float));
- }
+ mat->table_P_rho_u =
+ (float *)malloc(mat->num_rho * mat->num_u * sizeof(float *));
// Load table contents from file
FILE *f = fopen(table_file, "r");
int c;
for (int i = 0; i < mat->num_rho; i++) {
for (int j = 0; j < mat->num_u; j++) {
- c = fscanf(f, "%f", &mat->table_P_rho_u[i][j]);
+ c = fscanf(f, "%f", &mat->table_P_rho_u[i * mat->num_rho + j]);
if (c != 1) {
error("Failed to read EOS table");
}
@@ -147,7 +139,7 @@ INLINE static void convert_units_HM80(struct HM80_params *mat,
// Table Pressures in Mbar
for (int i = 0; i < mat->num_rho; i++) {
for (int j = 0; j < mat->num_u; j++) {
- mat->table_P_rho_u[i][j] *=
+ mat->table_P_rho_u[i * mat->num_rho + j] *=
Mbar_to_Ba / units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
}
}
@@ -194,7 +186,6 @@ INLINE static float HM80_soundspeed_from_entropy(
// gas_entropy_from_internal_energy
INLINE static float HM80_entropy_from_internal_energy(
float density, float u, const struct HM80_params *mat) {
- error("This EOS function is not yet implemented!");
return 0;
}
@@ -205,7 +196,7 @@ INLINE static float HM80_pressure_from_internal_energy(
float P;
- if (u <= 0) {
+ if (u <= 0.f) {
return 0.f;
}
@@ -226,32 +217,35 @@ INLINE static float HM80_pressure_from_internal_energy(
// Return zero pressure if below the table minimum/a
// Extrapolate the pressure for low densities
if (rho_idx < 0) { // Too-low rho
- P = expf(logf((1 - intp_u) * mat->table_P_rho_u[0][u_idx] +
- intp_u * mat->table_P_rho_u[0][u_idx + 1]) +
+ P = expf(logf((1 - intp_u) * mat->table_P_rho_u[u_idx] +
+ intp_u * mat->table_P_rho_u[u_idx + 1]) +
log_rho - mat->log_rho_min);
if (u_idx < 0) { // and too-low u
- P = 0;
+ P = 0.f;
}
} else if (u_idx < 0) { // Too-low u
- P = 0;
+ P = 0.f;
}
// Return an edge value if above the table maximum/a
else if (rho_idx >= mat->num_rho - 1) { // Too-high rho
if (u_idx >= mat->num_u - 1) { // and too-high u
- P = mat->table_P_rho_u[mat->num_rho - 1][mat->num_u - 1];
+ P = mat->table_P_rho_u[(mat->num_rho - 1) * mat->num_u + mat->num_u - 1];
} else {
- P = mat->table_P_rho_u[mat->num_rho - 1][u_idx];
+ P = mat->table_P_rho_u[(mat->num_rho - 1) * mat->num_u + u_idx];
}
} else if (u_idx >= mat->num_u - 1) { // Too-high u
- P = mat->table_P_rho_u[rho_idx][mat->num_u - 1];
+ P = mat->table_P_rho_u[rho_idx * mat->num_u + mat->num_u - 1];
}
// Normal interpolation within the table
else {
P = (1.f - intp_rho) *
- ((1.f - intp_u) * mat->table_P_rho_u[rho_idx][u_idx] +
- intp_u * mat->table_P_rho_u[rho_idx][u_idx + 1]) +
- intp_rho * ((1 - intp_u) * mat->table_P_rho_u[rho_idx + 1][u_idx] +
- intp_u * mat->table_P_rho_u[rho_idx + 1][u_idx + 1]);
+ ((1.f - intp_u) * mat->table_P_rho_u[rho_idx * mat->num_u + u_idx] +
+ intp_u * mat->table_P_rho_u[rho_idx * mat->num_u + u_idx + 1]) +
+ intp_rho *
+ ((1 - intp_u) *
+ mat->table_P_rho_u[(rho_idx + 1) * mat->num_u + u_idx] +
+ intp_u *
+ mat->table_P_rho_u[(rho_idx + 1) * mat->num_u + u_idx + 1]);
}
return P;
diff --git a/src/equation_of_state/planetary/tillotson.h b/src/equation_of_state/planetary/tillotson.h
index ed280079b040776dfc20480c50fbc9a3b6ff4470..d5b6d5c35d5edf9e114fe7f010c4f5b1e2327a83 100644
--- a/src/equation_of_state/planetary/tillotson.h
+++ b/src/equation_of_state/planetary/tillotson.h
@@ -147,7 +147,6 @@ INLINE static float Til_soundspeed_from_entropy(float density, float entropy,
// gas_entropy_from_internal_energy
INLINE static float Til_entropy_from_internal_energy(
float density, float u, const struct Til_params *mat) {
- error("This EOS function is not yet implemented!");
return 0;
}
diff --git a/src/hydro.h b/src/hydro.h
index cc6fa46ad374b39f5c9777874d2d85608b5eb2d6..950f63526a1590fa0fdcf2bfb5e650a2dfe14431 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -49,10 +49,14 @@
#include "./hydro/Default/hydro.h"
#include "./hydro/Default/hydro_iact.h"
#define SPH_IMPLEMENTATION "Default version of SPH"
-#elif defined(GIZMO_SPH)
-#include "./hydro/Gizmo/hydro.h"
-#include "./hydro/Gizmo/hydro_iact.h"
-#define SPH_IMPLEMENTATION "GIZMO (Hopkins 2015)"
+#elif defined(GIZMO_MFV_SPH)
+#include "./hydro/GizmoMFV/hydro.h"
+#include "./hydro/GizmoMFV/hydro_iact.h"
+#define SPH_IMPLEMENTATION "GIZMO MFV (Hopkins 2015)"
+#elif defined(GIZMO_MFM_SPH)
+#include "./hydro/GizmoMFM/hydro.h"
+#include "./hydro/GizmoMFM/hydro_iact.h"
+#define SPH_IMPLEMENTATION "GIZMO MFM (Hopkins 2015)"
#elif defined(SHADOWFAX_SPH)
#include "./hydro/Shadowswift/hydro.h"
#include "./hydro/Shadowswift/hydro_iact.h"
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/GizmoMFM/hydro.h
similarity index 99%
rename from src/hydro/Gizmo/hydro.h
rename to src/hydro/GizmoMFM/hydro.h
index 079e66669eb4e3b574742b09ca6b7e80d5386fb7..9c4be6af359d7236e483b712065b357c6ed35402 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/GizmoMFM/hydro.h
@@ -17,8 +17,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_GIZMO_HYDRO_H
-#define SWIFT_GIZMO_HYDRO_H
+#ifndef SWIFT_GIZMO_MFM_HYDRO_H
+#define SWIFT_GIZMO_MFM_HYDRO_H
#include "adiabatic_index.h"
#include "approx_math.h"
@@ -972,4 +972,4 @@ hydro_set_init_internal_energy(struct part* p, float u_init) {
p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u_init;
}
-#endif /* SWIFT_GIZMO_HYDRO_H */
+#endif /* SWIFT_GIZMO_MFM_HYDRO_H */
diff --git a/src/hydro/Gizmo/hydro_debug.h b/src/hydro/GizmoMFM/hydro_debug.h
similarity index 96%
rename from src/hydro/Gizmo/hydro_debug.h
rename to src/hydro/GizmoMFM/hydro_debug.h
index 0516068d3452e179c076f40a7c2db859be1c73c6..6603bc216b986b40513383120587d3caec1adc87 100644
--- a/src/hydro/Gizmo/hydro_debug.h
+++ b/src/hydro/GizmoMFM/hydro_debug.h
@@ -16,8 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_GIZMO_HYDRO_DEBUG_H
-#define SWIFT_GIZMO_HYDRO_DEBUG_H
+#ifndef SWIFT_GIZMO_MFM_HYDRO_DEBUG_H
+#define SWIFT_GIZMO_MFM_HYDRO_DEBUG_H
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
const struct part* p, const struct xpart* xp) {
@@ -78,4 +78,4 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
p->timestepvars.vmax, p->density.wcount_dh, p->density.wcount);
}
-#endif /* SWIFT_GIZMO_HYDRO_DEBUG_H */
+#endif /* SWIFT_GIZMO_MFM_HYDRO_DEBUG_H */
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/GizmoMFM/hydro_gradients.h
similarity index 97%
rename from src/hydro/Gizmo/hydro_gradients.h
rename to src/hydro/GizmoMFM/hydro_gradients.h
index c5d95e4dab4f574dd20e8355aacdd176530ec63d..964a2adcfe09b95c2a221af540e5e3ff0830dd67 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/GizmoMFM/hydro_gradients.h
@@ -17,8 +17,8 @@
*
******************************************************************************/
-#ifndef SWIFT_HYDRO_GIZMO_GRADIENTS_H
-#define SWIFT_HYDRO_GIZMO_GRADIENTS_H
+#ifndef SWIFT_HYDRO_GIZMO_MFM_GRADIENTS_H
+#define SWIFT_HYDRO_GIZMO_MFM_GRADIENTS_H
#include "hydro_slope_limiters.h"
#include "hydro_unphysical.h"
@@ -156,4 +156,4 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
Wj[3], Wj[4]);
}
-#endif /* SWIFT_HYDRO_GIZMO_GRADIENTS_H */
+#endif /* SWIFT_HYDRO_GIZMO_MFM_GRADIENTS_H */
diff --git a/src/hydro/Gizmo/hydro_gradients_gizmo.h b/src/hydro/GizmoMFM/hydro_gradients_gizmo.h
similarity index 99%
rename from src/hydro/Gizmo/hydro_gradients_gizmo.h
rename to src/hydro/GizmoMFM/hydro_gradients_gizmo.h
index 9f23f82e870cffb82a73964707bd28fa68ed00d7..1c3b68bb28375259628e09f16730710fbbd80149 100644
--- a/src/hydro/Gizmo/hydro_gradients_gizmo.h
+++ b/src/hydro/GizmoMFM/hydro_gradients_gizmo.h
@@ -22,8 +22,8 @@
*
* @param p Particle.
*/
-#ifndef SWIFT_GIZMO_HYDRO_GRADIENTS_H
-#define SWIFT_GIZMO_HYDRO_GRADIENTS_H
+#ifndef SWIFT_GIZMO_MFM_HYDRO_GRADIENTS_H
+#define SWIFT_GIZMO_MFM_HYDRO_GRADIENTS_H
__attribute__((always_inline)) INLINE static void hydro_gradients_init(
struct part *p) {
@@ -485,4 +485,4 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
hydro_slope_limit_cell(p);
}
-#endif /* SWIFT_GIZMO_HYDRO_GRADIENTS_H */
+#endif /* SWIFT_GIZMO_MFM_HYDRO_GRADIENTS_H */
diff --git a/src/hydro/Gizmo/hydro_gradients_sph.h b/src/hydro/GizmoMFM/hydro_gradients_sph.h
similarity index 98%
rename from src/hydro/Gizmo/hydro_gradients_sph.h
rename to src/hydro/GizmoMFM/hydro_gradients_sph.h
index fbaa056443df56691f1414f2d661ba9edc459734..169bed74f0b1b7e966f9880248f811d100bec13b 100644
--- a/src/hydro/Gizmo/hydro_gradients_sph.h
+++ b/src/hydro/GizmoMFM/hydro_gradients_sph.h
@@ -22,8 +22,8 @@
*
* @param p Particle.
*/
-#ifndef SWIFT_GIZMO_HYDRO_SPH_GRADIENTS_H
-#define SWIFT_GIZMO_HYDRO_SPH_GRADIENTS_H
+#ifndef SWIFT_GIZMO_MFM_HYDRO_SPH_GRADIENTS_H
+#define SWIFT_GIZMO_MFM_HYDRO_SPH_GRADIENTS_H
__attribute__((always_inline)) INLINE static void hydro_gradients_init(
struct part *p) {
@@ -253,4 +253,4 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
hydro_slope_limit_cell(p);
}
-#endif /* SWIFT_GIZMO_HYDRO_SPH_GRADIENTS_H */
+#endif /* SWIFT_GIZMO_MFM_HYDRO_SPH_GRADIENTS_H */
diff --git a/src/hydro/GizmoMFM/hydro_iact.h b/src/hydro/GizmoMFM/hydro_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..a1a82e514baa60d5895343ea84ef1c3aedc8a6b7
--- /dev/null
+++ b/src/hydro/GizmoMFM/hydro_iact.h
@@ -0,0 +1,517 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ *
+ * 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_GIZMO_MFM_HYDRO_IACT_H
+#define SWIFT_GIZMO_MFM_HYDRO_IACT_H
+
+#include "adiabatic_index.h"
+#include "hydro_gradients.h"
+#include "riemann.h"
+
+#define GIZMO_VOLUME_CORRECTION
+
+/**
+ * @brief Calculate the volume interaction between particle i and particle j
+ *
+ * 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.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_density(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ float wi, wj, wi_dx, wj_dx;
+
+ /* Get r and h inverse. */
+ const float r = sqrtf(r2);
+
+ /* Compute density of pi. */
+ const float hi_inv = 1.f / hi;
+ const float xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx);
+
+ /* these are eqns. (1) and (2) in the summary */
+ pi->geometry.volume += wi;
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
+
+ pi->geometry.centroid[0] -= dx[0] * wi;
+ pi->geometry.centroid[1] -= dx[1] * wi;
+ pi->geometry.centroid[2] -= dx[2] * wi;
+
+ /* Compute density of pj. */
+ const float hj_inv = 1.f / hj;
+ const float xj = r * hj_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ pj->density.wcount += wj;
+ pj->density.wcount_dh -= (hydro_dimension * wj + xj * wj_dx);
+
+ /* these are eqns. (1) and (2) in the summary */
+ pj->geometry.volume += wj;
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ pj->geometry.matrix_E[k][l] += dx[k] * dx[l] * wj;
+
+ pj->geometry.centroid[0] += dx[0] * wj;
+ pj->geometry.centroid[1] += dx[1] * wj;
+ pj->geometry.centroid[2] += dx[2] * wj;
+}
+
+/**
+ * @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.
+ *
+ * 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.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ const struct part *restrict pj, float a, float H) {
+
+ float wi, wi_dx;
+
+ /* Get r and h inverse. */
+ const float r = sqrtf(r2);
+
+ const float hi_inv = 1.f / hi;
+ const float xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx);
+
+ /* these are eqns. (1) and (2) in the summary */
+ pi->geometry.volume += wi;
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ pi->geometry.matrix_E[k][l] += dx[k] * dx[l] * wi;
+
+ pi->geometry.centroid[0] -= dx[0] * wi;
+ pi->geometry.centroid[1] -= dx[1] * wi;
+ pi->geometry.centroid[2] -= dx[2] * wi;
+}
+
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j
+ *
+ * This method wraps around hydro_gradients_collect, which can be an empty
+ * method, in which case no gradients are used.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_gradient(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ hydro_gradients_collect(r2, dx, hi, hj, pi, pj);
+}
+
+/**
+ * @brief Calculate the gradient interaction between particle i and particle j:
+ * non-symmetric version
+ *
+ * This method wraps around hydro_gradients_nonsym_collect, which can be an
+ * empty method, in which case no gradients are used.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_gradient(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ hydro_gradients_nonsym_collect(r2, dx, hi, hj, pi, pj);
+}
+
+/**
+ * @brief Common part of the flux calculation between particle i and j
+ *
+ * Since the only difference between the symmetric and non-symmetric version
+ * of the flux calculation is in the update of the conserved variables at the
+ * very end (which is not done for particle j if mode is 0), 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 also calculates the maximal velocity used to calculate the time
+ * step.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, int mode, float a, float H) {
+
+ const float r_inv = 1.f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ /* Initialize local variables */
+ float Bi[3][3];
+ float Bj[3][3];
+ float vi[3], vj[3];
+ for (int k = 0; k < 3; k++) {
+ for (int l = 0; l < 3; l++) {
+ Bi[k][l] = pi->geometry.matrix_E[k][l];
+ Bj[k][l] = pj->geometry.matrix_E[k][l];
+ }
+ vi[k] = pi->v[k]; /* particle velocities */
+ vj[k] = pj->v[k];
+ }
+ const float Vi = pi->geometry.volume;
+ const float Vj = pj->geometry.volume;
+ float Wi[5], Wj[5];
+ Wi[0] = pi->primitives.rho;
+ Wi[1] = pi->primitives.v[0];
+ Wi[2] = pi->primitives.v[1];
+ Wi[3] = pi->primitives.v[2];
+ Wi[4] = pi->primitives.P;
+ Wj[0] = pj->primitives.rho;
+ Wj[1] = pj->primitives.v[0];
+ Wj[2] = pj->primitives.v[1];
+ Wj[3] = pj->primitives.v[2];
+ Wj[4] = pj->primitives.P;
+
+ /* calculate the maximal signal velocity */
+ float vmax;
+ if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
+ const float ci = gas_soundspeed_from_pressure(Wi[0], Wi[4]);
+ const float cj = gas_soundspeed_from_pressure(Wj[0], Wj[4]);
+ vmax = ci + cj;
+ } else
+ vmax = 0.0f;
+
+ 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];
+
+ /* Velocity on the axis linking the particles */
+ float dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
+ (Wi[3] - Wj[3]) * dx[2];
+ dvdotdx = min(dvdotdx, dvdr);
+
+ /* We only care about this velocity for particles moving towards each others
+ */
+ dvdotdx = min(dvdotdx, 0.f);
+
+ /* Get the signal velocity */
+ /* the magical factor 3 also appears in Gadget2 */
+ vmax -= 3.f * dvdotdx * r_inv;
+
+ /* Store the signal velocity */
+ pi->timestepvars.vmax = max(pi->timestepvars.vmax, vmax);
+ if (mode == 1) pj->timestepvars.vmax = max(pj->timestepvars.vmax, vmax);
+
+ /* Compute kernel of pi. */
+ float wi, wi_dx;
+ const float hi_inv = 1.f / hi;
+ const float hi_inv_dim = pow_dimension(hi_inv);
+ const float xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ /* Compute kernel of pj. */
+ float wj, wj_dx;
+ const float hj_inv = 1.f / hj;
+ const float hj_inv_dim = pow_dimension(hj_inv);
+ const float xj = r * hj_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ /* Compute h_dt. We are going to use an SPH-like estimate of div_v for that */
+ const float hidp1 = pow_dimension_plus_one(hi_inv);
+ const float hjdp1 = pow_dimension_plus_one(hj_inv);
+ const float wi_dr = hidp1 * wi_dx;
+ const float wj_dr = hjdp1 * wj_dx;
+ dvdr *= r_inv;
+ if (pj->primitives.rho > 0.)
+ pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
+ if (mode == 1 && pi->primitives.rho > 0.)
+ pj->force.h_dt -= pi->conserved.mass * dvdr / pi->primitives.rho * wj_dr;
+
+ /* Compute (square of) area */
+ /* eqn. (7) */
+ float Anorm2 = 0.0f;
+ float A[3];
+ if (pi->density.wcorr > const_gizmo_min_wcorr &&
+ pj->density.wcorr > const_gizmo_min_wcorr) {
+ /* in principle, we use Vi and Vj as weights for the left and right
+ contributions to the generalized surface vector.
+ However, if Vi and Vj are very different (because they have very
+ different
+ smoothing lengths), then the expressions below are more stable. */
+ float Xi = Vi;
+ float Xj = Vj;
+#ifdef GIZMO_VOLUME_CORRECTION
+ if (fabsf(Vi - Vj) / min(Vi, Vj) > 1.5f * hydro_dimension) {
+ Xi = (Vi * hj + Vj * hi) / (hi + hj);
+ Xj = Xi;
+ }
+#endif
+ for (int k = 0; k < 3; k++) {
+ /* we add a minus sign since dx is pi->x - pj->x */
+ A[k] = -Xi * (Bi[k][0] * dx[0] + Bi[k][1] * dx[1] + Bi[k][2] * dx[2]) *
+ wi * hi_inv_dim -
+ Xj * (Bj[k][0] * dx[0] + Bj[k][1] * dx[1] + Bj[k][2] * dx[2]) *
+ wj * hj_inv_dim;
+ Anorm2 += A[k] * A[k];
+ }
+ } else {
+ /* ill condition gradient matrix: revert to SPH face area */
+ const float Anorm =
+ -(hidp1 * Vi * Vi * wi_dx + hjdp1 * Vj * Vj * wj_dx) * r_inv;
+ A[0] = -Anorm * dx[0];
+ A[1] = -Anorm * dx[1];
+ A[2] = -Anorm * dx[2];
+ Anorm2 = Anorm * Anorm * r2;
+ }
+
+ /* if the interface has no area, nothing happens and we return */
+ /* continuing results in dividing by zero and NaN's... */
+ if (Anorm2 == 0.f) return;
+
+ /* Compute the area */
+ const float Anorm_inv = 1. / sqrtf(Anorm2);
+ const float Anorm = Anorm2 * Anorm_inv;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* For stability reasons, we do require A and dx to have opposite
+ directions (basically meaning that the surface normal for the surface
+ always points from particle i to particle j, as it would in a real
+ moving-mesh code). If not, our scheme is no longer upwind and hence can
+ become unstable. */
+ const float dA_dot_dx = A[0] * dx[0] + A[1] * dx[1] + A[2] * dx[2];
+ /* In GIZMO, Phil Hopkins reverts to an SPH integration scheme if this
+ happens. We curently just ignore this case and display a message. */
+ const float rdim = pow_dimension(r);
+ if (dA_dot_dx > 1.e-6f * rdim) {
+ message("Ill conditioned gradient matrix (%g %g %g %g %g)!", dA_dot_dx,
+ Anorm, Vi, Vj, r);
+ }
+#endif
+
+ /* compute the normal vector of the interface */
+ const float n_unit[3] = {A[0] * Anorm_inv, A[1] * Anorm_inv,
+ A[2] * Anorm_inv};
+
+ /* Compute interface position (relative to pi, since we don't need the actual
+ * position) eqn. (8) */
+ const float xfac = -hi / (hi + hj);
+ const float xij_i[3] = {xfac * dx[0], xfac * dx[1], xfac * dx[2]};
+
+ /* Compute interface velocity */
+ /* eqn. (9) */
+ float xijdotdx = xij_i[0] * dx[0] + xij_i[1] * dx[1] + xij_i[2] * dx[2];
+ xijdotdx *= r_inv * r_inv;
+ const float vij[3] = {vi[0] + (vi[0] - vj[0]) * xijdotdx,
+ vi[1] + (vi[1] - vj[1]) * xijdotdx,
+ vi[2] + (vi[2] - vj[2]) * xijdotdx};
+
+ /* complete calculation of position of interface */
+ /* NOTE: dx is not necessarily just pi->x - pj->x but can also contain
+ correction terms for periodicity. If we do the interpolation,
+ we have to use xij w.r.t. the actual particle.
+ => we need a separate xij for pi and pj... */
+ /* tldr: we do not need the code below, but we do need the same code as above
+ but then with i and j swapped */
+ // for ( k = 0 ; k < 3 ; k++ )
+ // xij[k] += pi->x[k];
+
+ /* Boost the primitive variables to the frame of reference of the interface */
+ /* Note that velocities are indices 1-3 in W */
+ Wi[1] -= vij[0];
+ Wi[2] -= vij[1];
+ Wi[3] -= vij[2];
+ Wj[1] -= vij[0];
+ Wj[2] -= vij[1];
+ Wj[3] -= vij[2];
+
+ hydro_gradients_predict(pi, pj, hi, hj, dx, r, xij_i, Wi, Wj);
+
+ /* we don't need to rotate, we can use the unit vector in the Riemann problem
+ * itself (see GIZMO) */
+
+ float totflux[5];
+ riemann_solve_for_middle_state_flux(Wi, Wj, n_unit, vij, totflux);
+
+ /* Multiply with the interface surface area */
+ totflux[0] *= Anorm;
+ totflux[1] *= Anorm;
+ totflux[2] *= Anorm;
+ totflux[3] *= Anorm;
+ totflux[4] *= Anorm;
+
+ /* Store mass flux */
+ const float mflux_i = totflux[0];
+ pi->gravity.mflux[0] += mflux_i * dx[0];
+ pi->gravity.mflux[1] += mflux_i * dx[1];
+ pi->gravity.mflux[2] += mflux_i * dx[2];
+
+ /* Update conserved variables */
+ /* eqn. (16) */
+ pi->conserved.flux.mass -= totflux[0];
+ pi->conserved.flux.momentum[0] -= totflux[1];
+ pi->conserved.flux.momentum[1] -= totflux[2];
+ pi->conserved.flux.momentum[2] -= totflux[3];
+ pi->conserved.flux.energy -= totflux[4];
+
+#ifndef GIZMO_TOTAL_ENERGY
+ const float ekin_i = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
+ pi->primitives.v[1] * pi->primitives.v[1] +
+ pi->primitives.v[2] * pi->primitives.v[2]);
+ pi->conserved.flux.energy += totflux[1] * pi->primitives.v[0];
+ pi->conserved.flux.energy += totflux[2] * pi->primitives.v[1];
+ pi->conserved.flux.energy += totflux[3] * pi->primitives.v[2];
+ pi->conserved.flux.energy -= totflux[0] * ekin_i;
+#endif
+
+ /* Note that this used to be much more complicated in early implementations of
+ * the GIZMO scheme, as we wanted manifest conservation of conserved variables
+ * and had to do symmetric flux exchanges. Now we don't care about manifest
+ * conservation anymore and just assume the current fluxes are representative
+ * for the flux over the entire time step. */
+ if (mode == 1) {
+ /* Store mass flux */
+ const float mflux_j = totflux[0];
+ pj->gravity.mflux[0] -= mflux_j * dx[0];
+ pj->gravity.mflux[1] -= mflux_j * dx[1];
+ pj->gravity.mflux[2] -= mflux_j * dx[2];
+
+ pj->conserved.flux.mass += totflux[0];
+ pj->conserved.flux.momentum[0] += totflux[1];
+ pj->conserved.flux.momentum[1] += totflux[2];
+ pj->conserved.flux.momentum[2] += totflux[3];
+ pj->conserved.flux.energy += totflux[4];
+
+#ifndef GIZMO_TOTAL_ENERGY
+ const float ekin_j = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
+ pj->primitives.v[1] * pj->primitives.v[1] +
+ pj->primitives.v[2] * pj->primitives.v[2]);
+ pj->conserved.flux.energy -= totflux[1] * pj->primitives.v[0];
+ pj->conserved.flux.energy -= totflux[2] * pj->primitives.v[1];
+ pj->conserved.flux.energy -= totflux[3] * pj->primitives.v[2];
+ pj->conserved.flux.energy += totflux[0] * ekin_j;
+#endif
+ }
+}
+
+/**
+ * @brief Flux calculation between particle i and particle j
+ *
+ * This method calls runner_iact_fluxes_common with mode 1.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_force(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1, a, H);
+}
+
+/**
+ * @brief Flux calculation between particle i and particle j: non-symmetric
+ * version
+ *
+ * This method calls runner_iact_fluxes_common with mode 0.
+ *
+ * @param r2 Comoving squared distance between particle i and particle j.
+ * @param dx Comoving distance vector between the particles (dx = pi->x -
+ * pj->x).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @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 *restrict pi,
+ struct part *restrict pj, float a, float H) {
+
+ runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0, a, H);
+}
+
+#endif /* SWIFT_GIZMO_MFM_HYDRO_IACT_H */
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/GizmoMFM/hydro_io.h
similarity index 98%
rename from src/hydro/Gizmo/hydro_io.h
rename to src/hydro/GizmoMFM/hydro_io.h
index 7aa56b3bde5247f4839b8f4f60cdb3de67253392..171132eacfcd43feeec57d3c16b5a458171b1d79 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/GizmoMFM/hydro_io.h
@@ -16,8 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_GIZMO_HYDRO_IO_H
-#define SWIFT_GIZMO_HYDRO_IO_H
+#ifndef SWIFT_GIZMO_MFM_HYDRO_IO_H
+#define SWIFT_GIZMO_MFM_HYDRO_IO_H
#include "adiabatic_index.h"
#include "hydro.h"
@@ -241,4 +241,4 @@ void hydro_write_flavour(hid_t h_grpsph) {
*/
int writeEntropyFlag() { return 0; }
-#endif /* SWIFT_GIZMO_HYDRO_IO_H */
+#endif /* SWIFT_GIZMO_MFM_HYDRO_IO_H */
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/GizmoMFM/hydro_part.h
similarity index 97%
rename from src/hydro/Gizmo/hydro_part.h
rename to src/hydro/GizmoMFM/hydro_part.h
index d5bc3f84238fecbf3b2f120f031996f55bbd65fe..857c429ec933ad3eea730e8f0b6f830782cdf77b 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/GizmoMFM/hydro_part.h
@@ -16,8 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_GIZMO_HYDRO_PART_H
-#define SWIFT_GIZMO_HYDRO_PART_H
+#ifndef SWIFT_GIZMO_MFM_HYDRO_PART_H
+#define SWIFT_GIZMO_MFM_HYDRO_PART_H
#include "chemistry_struct.h"
#include "cooling_struct.h"
@@ -210,4 +210,4 @@ struct part {
} SWIFT_STRUCT_ALIGN;
-#endif /* SWIFT_GIZMO_HYDRO_PART_H */
+#endif /* SWIFT_GIZMO_MFM_HYDRO_PART_H */
diff --git a/src/hydro/Gizmo/hydro_slope_limiters.h b/src/hydro/GizmoMFM/hydro_slope_limiters.h
similarity index 100%
rename from src/hydro/Gizmo/hydro_slope_limiters.h
rename to src/hydro/GizmoMFM/hydro_slope_limiters.h
diff --git a/src/hydro/Gizmo/hydro_slope_limiters_cell.h b/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h
similarity index 98%
rename from src/hydro/Gizmo/hydro_slope_limiters_cell.h
rename to src/hydro/GizmoMFM/hydro_slope_limiters_cell.h
index 599be632b64e3b06ab1e0339a6ad4b560b2abece..7dec6f499da31de1f10652a31781a788166957cc 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters_cell.h
+++ b/src/hydro/GizmoMFM/hydro_slope_limiters_cell.h
@@ -16,8 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_GIZMO_SLOPE_LIMITER_CELL_H
-#define SWIFT_GIZMO_SLOPE_LIMITER_CELL_H
+#ifndef SWIFT_GIZMO_MFM_SLOPE_LIMITER_CELL_H
+#define SWIFT_GIZMO_MFM_SLOPE_LIMITER_CELL_H
#include <float.h>
@@ -183,4 +183,4 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
}
}
-#endif /* SWIFT_GIZMO_SLOPE_LIMITER_CELL_H */
+#endif /* SWIFT_GIZMO_MFM_SLOPE_LIMITER_CELL_H */
diff --git a/src/hydro/Gizmo/hydro_slope_limiters_face.h b/src/hydro/GizmoMFM/hydro_slope_limiters_face.h
similarity index 97%
rename from src/hydro/Gizmo/hydro_slope_limiters_face.h
rename to src/hydro/GizmoMFM/hydro_slope_limiters_face.h
index 39a453b3d21383fbf075a3c1985dfb182d9a2dcc..8f0636c992ccbd615cc62ca09c1b0ca9b08ea56b 100644
--- a/src/hydro/Gizmo/hydro_slope_limiters_face.h
+++ b/src/hydro/GizmoMFM/hydro_slope_limiters_face.h
@@ -29,8 +29,8 @@
* @return The slope limited difference between the quantity at the particle
* position and the quantity at the interface position.
*/
-#ifndef SWIFT_GIZMO_SLOPE_LIMITER_FACE_H
-#define SWIFT_GIZMO_SLOPE_LIMITER_FACE_H
+#ifndef SWIFT_GIZMO_MFM_SLOPE_LIMITER_FACE_H
+#define SWIFT_GIZMO_MFM_SLOPE_LIMITER_FACE_H
/* Some standard headers. */
#include <float.h>
@@ -125,4 +125,4 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
xij_j_norm, r_inv);
}
-#endif /* SWIFT_GIZMO_SLOPE_LIMITER_FACE_H */
+#endif /* SWIFT_GIZMO_MFM_SLOPE_LIMITER_FACE_H */
diff --git a/src/hydro/Gizmo/hydro_unphysical.h b/src/hydro/GizmoMFM/hydro_unphysical.h
similarity index 100%
rename from src/hydro/Gizmo/hydro_unphysical.h
rename to src/hydro/GizmoMFM/hydro_unphysical.h
diff --git a/src/hydro/Gizmo/hydro_velocities.h b/src/hydro/GizmoMFM/hydro_velocities.h
similarity index 100%
rename from src/hydro/Gizmo/hydro_velocities.h
rename to src/hydro/GizmoMFM/hydro_velocities.h
diff --git a/src/hydro/GizmoMFV/hydro.h b/src/hydro/GizmoMFV/hydro.h
new file mode 100644
index 0000000000000000000000000000000000000000..1d5abeaaf63d88b02817a691f160c537d0b1915b
--- /dev/null
+++ b/src/hydro/GizmoMFV/hydro.h
@@ -0,0 +1,975 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * 2016, 2017 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_MFV_HYDRO_H
+#define SWIFT_GIZMO_MFV_HYDRO_H
+
+#include "adiabatic_index.h"
+#include "approx_math.h"
+#include "cosmology.h"
+#include "equation_of_state.h"
+#include "hydro_gradients.h"
+#include "hydro_properties.h"
+#include "hydro_space.h"
+#include "hydro_unphysical.h"
+#include "hydro_velocities.h"
+#include "minmax.h"
+#include "riemann.h"
+
+#include <float.h>
+
+//#define GIZMO_LLOYD_ITERATION
+
+/**
+ * @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.
+ * @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;
+
+#ifdef GIZMO_LLOYD_ITERATION
+ return CFL_condition;
+#endif
+
+ /* v_full is the actual velocity of the particle, primitives.v is its
+ hydrodynamical velocity. The time step depends on the relative difference
+ of the two. */
+ float vrel[3];
+ vrel[0] = p->primitives.v[0] - xp->v_full[0];
+ vrel[1] = p->primitives.v[1] - xp->v_full[1];
+ vrel[2] = p->primitives.v[2] - xp->v_full[2];
+ float vmax =
+ sqrtf(vrel[0] * vrel[0] + vrel[1] * vrel[1] + vrel[2] * vrel[2]) +
+ sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+ vmax = max(vmax, p->timestepvars.vmax);
+
+ // MATTHIEU: Bert is this correct? Do we need more cosmology terms here?
+ const float psize =
+ cosmo->a * powf(p->geometry.volume / hydro_dimension_unit_sphere,
+ hydro_dimension_inv);
+ float dt = FLT_MAX;
+ if (vmax > 0.) {
+ dt = psize / vmax;
+ }
+ return CFL_condition * dt;
+}
+
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * This method is no longer used, as Gizmo is now unaware of the actual particle
+ * time step.
+ *
+ * @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) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (dt == 0.) {
+ error("Zero time step assigned to particle!");
+ }
+
+ if (dt != dt) {
+ error("NaN time step assigned to particle!");
+ }
+#endif
+}
+
+/**
+ * @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];
+
+ /* we can already initialize the momentum */
+ 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];
+
+/* and the thermal energy */
+/* remember that we store the total thermal energy, not the specific thermal
+ energy (as in Gadget) */
+#if defined(EOS_ISOTHERMAL_GAS)
+ /* this overwrites the internal energy from the initial condition file
+ * Note that we call the EoS function just to get the constant u here. */
+ p->conserved.energy = mass * gas_internal_energy_from_entropy(0.f, 0.f);
+#else
+ p->conserved.energy *= mass;
+#endif
+
+#ifdef GIZMO_TOTAL_ENERGY
+ /* add the total kinetic 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
+
+#ifdef GIZMO_LLOYD_ITERATION
+ /* overwrite all variables to make sure they have safe values */
+ p->primitives.rho = 1.;
+ p->primitives.v[0] = 0.;
+ p->primitives.v[1] = 0.;
+ p->primitives.v[2] = 0.;
+ p->primitives.P = 1.;
+
+ p->conserved.mass = 1.;
+ p->conserved.momentum[0] = 0.;
+ p->conserved.momentum[1] = 0.;
+ p->conserved.momentum[2] = 0.;
+ p->conserved.energy = 1.;
+
+ p->v[0] = 0.;
+ p->v[1] = 0.;
+ p->v[2] = 0.;
+#endif
+
+ /* initialize the particle velocity based on the primitive fluid velocity */
+ hydro_velocities_init(p, xp);
+
+ /* ignore accelerations present in the initial condition */
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+
+ /* we cannot initialize wcorr in init_part, as init_part gets called every
+ time the density loop is repeated, and the whole point of storing wcorr
+ is to have a way of remembering that we need more neighbours for this
+ particle */
+ p->density.wcorr = 1.0f;
+}
+
+/**
+ * @brief Prepares a particle for the volume calculation.
+ *
+ * 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, const struct hydro_space* hs) {
+
+ p->density.wcount = 0.0f;
+ p->density.wcount_dh = 0.0f;
+ p->geometry.volume = 0.0f;
+ p->geometry.matrix_E[0][0] = 0.0f;
+ p->geometry.matrix_E[0][1] = 0.0f;
+ p->geometry.matrix_E[0][2] = 0.0f;
+ p->geometry.matrix_E[1][0] = 0.0f;
+ p->geometry.matrix_E[1][1] = 0.0f;
+ p->geometry.matrix_E[1][2] = 0.0f;
+ p->geometry.matrix_E[2][0] = 0.0f;
+ p->geometry.matrix_E[2][1] = 0.0f;
+ p->geometry.matrix_E[2][2] = 0.0f;
+ p->geometry.centroid[0] = 0.0f;
+ p->geometry.centroid[1] = 0.0f;
+ p->geometry.centroid[2] = 0.0f;
+}
+
+/**
+ * @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.
+ *
+ * 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.
+ *
+ * @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) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float ih = 1.0f / h;
+ const float ihdim = pow_dimension(ih);
+ const float ihdim_plus_one = ihdim * ih;
+
+ /* Final operation on the density. */
+ p->density.wcount += kernel_root;
+ p->density.wcount *= ihdim;
+
+ p->density.wcount_dh -= hydro_dimension * kernel_root;
+ p->density.wcount_dh *= ihdim_plus_one;
+
+ /* Final operation on the geometry. */
+ /* we multiply with the smoothing kernel normalization ih3 and calculate the
+ * volume */
+ const float volume = 1.f / (ihdim * (p->geometry.volume + kernel_root));
+ p->geometry.volume = volume;
+
+ /* we multiply with the smoothing kernel normalization */
+ p->geometry.matrix_E[0][0] = ihdim * p->geometry.matrix_E[0][0];
+ p->geometry.matrix_E[0][1] = ihdim * p->geometry.matrix_E[0][1];
+ p->geometry.matrix_E[0][2] = ihdim * p->geometry.matrix_E[0][2];
+ p->geometry.matrix_E[1][0] = ihdim * p->geometry.matrix_E[1][0];
+ p->geometry.matrix_E[1][1] = ihdim * p->geometry.matrix_E[1][1];
+ p->geometry.matrix_E[1][2] = ihdim * p->geometry.matrix_E[1][2];
+ p->geometry.matrix_E[2][0] = ihdim * p->geometry.matrix_E[2][0];
+ p->geometry.matrix_E[2][1] = ihdim * p->geometry.matrix_E[2][1];
+ p->geometry.matrix_E[2][2] = ihdim * p->geometry.matrix_E[2][2];
+
+ p->geometry.centroid[0] *= kernel_norm;
+ p->geometry.centroid[1] *= kernel_norm;
+ p->geometry.centroid[2] *= kernel_norm;
+
+ p->geometry.centroid[0] /= p->density.wcount;
+ p->geometry.centroid[1] /= p->density.wcount;
+ p->geometry.centroid[2] /= p->density.wcount;
+
+ /* Check the condition number to see if we have a stable geometry. */
+ float condition_number_E = 0.0f;
+ int i, j;
+ for (i = 0; i < 3; ++i) {
+ for (j = 0; j < 3; ++j) {
+ condition_number_E +=
+ p->geometry.matrix_E[i][j] * p->geometry.matrix_E[i][j];
+ }
+ }
+
+ invert_dimension_by_dimension_matrix(p->geometry.matrix_E);
+
+ float condition_number_Einv = 0.0f;
+ for (i = 0; i < 3; ++i) {
+ for (j = 0; j < 3; ++j) {
+ condition_number_Einv +=
+ p->geometry.matrix_E[i][j] * p->geometry.matrix_E[i][j];
+ }
+ }
+
+ float condition_number =
+ hydro_dimension_inv * sqrtf(condition_number_E * condition_number_Einv);
+
+ if (condition_number > const_gizmo_max_condition_number &&
+ p->density.wcorr > const_gizmo_min_wcorr) {
+#ifdef GIZMO_PATHOLOGICAL_ERROR
+ error("Condition number larger than %g (%g)!",
+ const_gizmo_max_condition_number, condition_number);
+#endif
+#ifdef GIZMO_PATHOLOGICAL_WARNING
+ message("Condition number too large: %g (> %g, p->id: %llu)!",
+ condition_number, const_gizmo_max_condition_number, p->id);
+#endif
+ /* add a correction to the number of neighbours for this particle */
+ p->density.wcorr *= const_gizmo_w_correction_factor;
+ }
+
+ hydro_gradients_init(p);
+
+ /* compute primitive variables */
+ /* eqns (3)-(5) */
+ const float m = p->conserved.mass;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (m < 0.) {
+ error("Mass is negative!");
+ }
+
+ if (volume == 0.) {
+ error("Volume is 0!");
+ }
+#endif
+
+ // MATTHIEU: Bert is this correct? Do we need cosmology terms here?
+ float momentum[3];
+ momentum[0] = p->conserved.momentum[0];
+ momentum[1] = p->conserved.momentum[1];
+ momentum[2] = p->conserved.momentum[2];
+ p->primitives.rho = m / volume;
+ if (m == 0.) {
+ p->primitives.v[0] = 0.;
+ p->primitives.v[1] = 0.;
+ p->primitives.v[2] = 0.;
+ } else {
+ p->primitives.v[0] = momentum[0] / m;
+ p->primitives.v[1] = momentum[1] / m;
+ p->primitives.v[2] = momentum[2] / m;
+ }
+
+#ifdef EOS_ISOTHERMAL_GAS
+ /* although the pressure is not formally used anywhere if an isothermal eos
+ has been selected, we still make sure it is set to the correct value */
+ p->primitives.P = gas_pressure_from_internal_energy(p->primitives.rho, 0.);
+#else
+
+ float energy = p->conserved.energy;
+
+#ifdef GIZMO_TOTAL_ENERGY
+ /* subtract the kinetic energy; we want the thermal energy */
+ energy -= 0.5f * (momentum[0] * p->primitives.v[0] +
+ momentum[1] * p->primitives.v[1] +
+ momentum[2] * p->primitives.v[2]);
+#endif
+
+ /* energy contains the total thermal energy, we want the specific energy.
+ this is why we divide by the volume, and not by the density */
+ p->primitives.P = hydro_gamma_minus_one * energy / volume;
+#endif
+
+ /* sanity checks */
+ gizmo_check_physical_quantities("density", "pressure", p->primitives.rho,
+ p->primitives.v[0], p->primitives.v[1],
+ p->primitives.v[2], p->primitives.P);
+
+#ifdef GIZMO_LLOYD_ITERATION
+ /* overwrite primitive variables to make sure they still have safe values */
+ p->primitives.rho = 1.;
+ p->primitives.v[0] = 0.;
+ p->primitives.v[1] = 0.;
+ p->primitives.v[2] = 0.;
+ p->primitives.P = 1.;
+#endif
+
+ /* Add a correction factor to wcount (to force a neighbour number increase if
+ the geometry matrix is close to singular) */
+ p->density.wcount *= p->density.wcorr;
+ p->density.wcount_dh *= p->density.wcorr;
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * @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,
+ const struct cosmology* cosmo) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+
+ /* Re-set problematic values */
+ p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+ p->density.wcount_dh = 0.f;
+ p->geometry.volume = 1.0f;
+ p->geometry.matrix_E[0][0] = 1.0f;
+ p->geometry.matrix_E[0][1] = 0.0f;
+ p->geometry.matrix_E[0][2] = 0.0f;
+ p->geometry.matrix_E[1][0] = 0.0f;
+ p->geometry.matrix_E[1][1] = 1.0f;
+ p->geometry.matrix_E[1][2] = 0.0f;
+ p->geometry.matrix_E[2][0] = 0.0f;
+ p->geometry.matrix_E[2][1] = 0.0f;
+ p->geometry.matrix_E[2][2] = 1.0f;
+ /* centroid is relative w.r.t. particle position */
+ /* by setting the centroid to 0.0f, we make sure no velocity correction is
+ applied */
+ p->geometry.centroid[0] = 0.0f;
+ p->geometry.centroid[1] = 0.0f;
+ p->geometry.centroid[2] = 0.0f;
+}
+
+/**
+ * @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.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+ struct part* restrict p, struct xpart* restrict xp,
+ const struct cosmology* cosmo) {
+
+ /* Initialize time step criterion variables */
+ p->timestepvars.vmax = 0.;
+
+ // MATTHIEU: Bert is this correct? Do we need cosmology terms here?
+
+ /* Set the actual velocity of the particle */
+ hydro_velocities_prepare_force(p, xp);
+}
+
+/**
+ * @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.
+ *
+ * This method also initializes the force loop variables.
+ *
+ * @param p The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_gradient(
+ struct part* p) {
+
+ hydro_gradients_finalize(p);
+
+ p->gravity.mflux[0] = 0.0f;
+ p->gravity.mflux[1] = 0.0f;
+ p->gravity.mflux[2] = 0.0f;
+
+ p->conserved.flux.mass = 0.0f;
+ p->conserved.flux.momentum[0] = 0.0f;
+ p->conserved.flux.momentum[1] = 0.0f;
+ p->conserved.flux.momentum[2] = 0.0f;
+ p->conserved.flux.energy = 0.0f;
+
+#ifdef GIZMO_LLOYD_ITERATION
+ /* reset the gradients to zero, as we don't want them */
+ hydro_gradients_init(p);
+#endif
+}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * This is actually not necessary for GIZMO, 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
+ *
+ * We no longer do this, as the mass needs to be provided in the initial
+ * condition file, and the mass alone is enough to initialize all conserved
+ * variables. This is now done in hydro_first_init_part.
+ *
+ * @param p The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+ struct part* p, struct xpart* xp, const struct cosmology* cosmo) {}
+
+/**
+ * @brief Extra operations to be done during the drift
+ *
+ * @param p Particle to act upon.
+ * @param xp The extended particle data to act upon.
+ * @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* p, struct xpart* xp, float dt_drift, float dt_therm) {
+
+#ifdef GIZMO_LLOYD_ITERATION
+ return;
+#endif
+
+ const float h_inv = 1.0f / p->h;
+
+ /* Predict smoothing length */
+ const float w1 = p->force.h_dt * h_inv * dt_drift;
+ float h_corr;
+ if (fabsf(w1) < 0.2f)
+ h_corr = approx_expf(w1); /* 4th order expansion of exp(w) */
+ else
+ h_corr = expf(w1);
+
+ /* Limit the smoothing length correction (and make sure it is always
+ positive). */
+ if (h_corr < 2.0f && h_corr > 0.) {
+ p->h *= h_corr;
+ }
+
+ /* drift the primitive variables based on the old fluxes */
+ if (p->geometry.volume > 0.) {
+ p->primitives.rho += p->conserved.flux.mass * dt_drift / p->geometry.volume;
+ }
+
+ if (p->conserved.mass > 0.) {
+ p->primitives.v[0] +=
+ p->conserved.flux.momentum[0] * dt_drift / p->conserved.mass;
+ p->primitives.v[1] +=
+ p->conserved.flux.momentum[1] * dt_drift / p->conserved.mass;
+ p->primitives.v[2] +=
+ p->conserved.flux.momentum[2] * dt_drift / p->conserved.mass;
+
+#if !defined(EOS_ISOTHERMAL_GAS)
+ const float u = p->conserved.energy + p->conserved.flux.energy * dt_therm;
+ p->primitives.P =
+ hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass;
+#endif
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (p->h <= 0.) {
+ error("Zero or negative smoothing length (%g)!", p->h);
+ }
+#endif
+
+ gizmo_check_physical_quantities("density", "pressure", p->primitives.rho,
+ p->primitives.v[0], p->primitives.v[1],
+ p->primitives.v[2], p->primitives.P);
+}
+
+/**
+ * @brief Set the particle acceleration after the flux loop
+ *
+ * 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.
+ *
+ * If the particle time step is zero, we set the accelerations to zero. This
+ * should only happen at the start of the simulation.
+ *
+ * @param p Particle to act upon.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_force(
+ struct part* p, const struct cosmology* cosmo) {
+
+ /* set the variables that are used to drift the primitive variables */
+
+ // MATTHIEU: Bert is this correct? Do we need cosmology terms here?
+ hydro_velocities_end_force(p);
+}
+
+/**
+ * @brief Extra operations done during the kick
+ *
+ * Not used for GIZMO.
+ *
+ * @param p Particle to act upon.
+ * @param xp Extended particle data to act upon.
+ * @param dt Physical time step.
+ * @param half_dt Half the physical time step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+ struct part* p, struct xpart* xp, float dt, const struct cosmology* cosmo,
+ const struct hydro_props* hydro_props) {
+
+ float a_grav[3];
+
+ /* Update conserved variables. */
+ p->conserved.mass += p->conserved.flux.mass * dt;
+ p->conserved.momentum[0] += p->conserved.flux.momentum[0] * dt;
+ p->conserved.momentum[1] += p->conserved.flux.momentum[1] * dt;
+ p->conserved.momentum[2] += p->conserved.flux.momentum[2] * dt;
+#if defined(EOS_ISOTHERMAL_GAS)
+ /* We use the EoS equation in a sneaky way here just to get the constant u */
+ p->conserved.energy =
+ p->conserved.mass * gas_internal_energy_from_entropy(0.f, 0.f);
+#else
+ p->conserved.energy += p->conserved.flux.energy * dt;
+#endif
+
+ /* Apply the minimal energy limit */
+ const float min_energy =
+ hydro_props->minimal_internal_energy * cosmo->a_factor_internal_energy;
+ if (p->conserved.energy < min_energy * p->conserved.mass) {
+ p->conserved.energy = min_energy * p->conserved.mass;
+ p->conserved.flux.energy = 0.f;
+ }
+
+ gizmo_check_physical_quantities(
+ "mass", "energy", p->conserved.mass, p->conserved.momentum[0],
+ p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.energy);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Note that this check will only have effect if no GIZMO_UNPHYSICAL option
+ was selected. */
+ if (p->conserved.mass < 0.) {
+ error(
+ "Negative mass after conserved variables update (mass: %g, dmass: %g)!",
+ p->conserved.mass, p->conserved.flux.mass);
+ }
+
+ if (p->conserved.energy < 0.) {
+ error(
+ "Negative energy after conserved variables update (energy: %g, "
+ "denergy: %g)!",
+ p->conserved.energy, p->conserved.flux.energy);
+ }
+#endif
+
+ /* Add gravity. We only do this if we have gravity activated. */
+ if (p->gpart) {
+ /* Retrieve the current value of the gravitational acceleration from the
+ gpart. We are only allowed to do this because this is the kick. We still
+ need to check whether gpart exists though.*/
+ a_grav[0] = p->gpart->a_grav[0];
+ a_grav[1] = p->gpart->a_grav[1];
+ a_grav[2] = p->gpart->a_grav[2];
+
+ /* Make sure the gpart knows the mass has changed. */
+ p->gpart->mass = p->conserved.mass;
+
+ /* Kick the momentum for half a time step */
+ /* Note that this also affects the particle movement, as the velocity for
+ the particles is set after this. */
+ p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0];
+ p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
+ p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
+
+ p->conserved.energy += dt * (p->gravity.mflux[0] * a_grav[0] +
+ p->gravity.mflux[1] * a_grav[1] +
+ p->gravity.mflux[2] * a_grav[2]);
+ }
+
+ hydro_velocities_set(p, xp);
+
+#ifdef GIZMO_LLOYD_ITERATION
+ /* reset conserved variables to safe values */
+ p->conserved.mass = 1.;
+ p->conserved.momentum[0] = 0.;
+ p->conserved.momentum[1] = 0.;
+ p->conserved.momentum[2] = 0.;
+ p->conserved.energy = 1.;
+
+ /* set the particle velocities to the Lloyd velocities */
+ /* note that centroid is the relative position of the centroid w.r.t. the
+ particle position (position - centroid) */
+ xp->v_full[0] = -p->geometry.centroid[0] / p->force.dt;
+ xp->v_full[1] = -p->geometry.centroid[1] / p->force.dt;
+ xp->v_full[2] = -p->geometry.centroid[2] / p->force.dt;
+ p->v[0] = xp->v_full[0];
+ p->v[1] = xp->v_full[1];
+ p->v[2] = xp->v_full[2];
+#endif
+
+ /* reset wcorr */
+ p->density.wcorr = 1.0f;
+}
+
+/**
+ * @brief Returns the comoving internal energy of a particle
+ *
+ * @param p The particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_internal_energy(const struct part* restrict p) {
+
+ if (p->primitives.rho > 0.)
+ return gas_internal_energy_from_pressure(p->primitives.rho,
+ p->primitives.P);
+ else
+ return 0.;
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle
+ *
+ * @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) {
+
+ return cosmo->a_factor_internal_energy *
+ hydro_get_comoving_internal_energy(p);
+}
+
+/**
+ * @brief Returns the comoving entropy of a particle
+ *
+ * @param p The particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy(
+ const struct part* restrict p) {
+
+ if (p->primitives.rho > 0.) {
+ return gas_entropy_from_pressure(p->primitives.rho, p->primitives.P);
+ } else {
+ return 0.;
+ }
+}
+
+/**
+ * @brief Returns the physical internal energy of a particle
+ *
+ * @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: no cosmological conversion required here with our choice of
+ * coordinates. */
+ return hydro_get_comoving_entropy(p);
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest.
+ */
+__attribute__((always_inline)) INLINE static float
+hydro_get_comoving_soundspeed(const struct part* restrict p) {
+
+ if (p->primitives.rho > 0.)
+ return gas_soundspeed_from_pressure(p->primitives.rho, p->primitives.P);
+ else
+ return 0.;
+}
+
+/**
+ * @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) {
+
+ return cosmo->a_factor_sound_speed * hydro_get_comoving_soundspeed(p);
+}
+
+/**
+ * @brief Returns the comoving pressure of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
+ const struct part* restrict p) {
+
+ return p->primitives.P;
+}
+
+/**
+ * @brief Returns the comoving pressure of a particle
+ *
+ * @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) {
+
+ return cosmo->a_factor_pressure * p->primitives.P;
+}
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_mass(
+ const struct part* restrict p) {
+
+ return p->conserved.mass;
+}
+
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part* restrict p, float m) {
+
+ p->conserved.mass = m;
+}
+
+/**
+ * @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_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(
+ const struct part* restrict p, const struct xpart* xp, float dt_kick_hydro,
+ float dt_kick_grav, float v[3]) {
+
+ if (p->conserved.mass > 0.) {
+ v[0] = p->primitives.v[0] +
+ p->conserved.flux.momentum[0] * dt_kick_hydro / p->conserved.mass;
+ v[1] = p->primitives.v[1] +
+ p->conserved.flux.momentum[1] * dt_kick_hydro / p->conserved.mass;
+ v[2] = p->primitives.v[2] +
+ p->conserved.flux.momentum[2] * dt_kick_hydro / p->conserved.mass;
+ } else {
+ v[0] = p->primitives.v[0];
+ v[1] = p->primitives.v[1];
+ v[2] = p->primitives.v[2];
+ }
+
+ // MATTHIEU: Bert is this correct?
+ v[0] += xp->a_grav[0] * dt_kick_grav;
+ v[1] += xp->a_grav[1] * dt_kick_grav;
+ v[2] += xp->a_grav[2] * dt_kick_grav;
+}
+
+/**
+ * @brief Returns the comoving density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_comoving_density(
+ const struct part* restrict p) {
+
+ return p->primitives.rho;
+}
+
+/**
+ * @brief Returns the physical 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* cosmo) {
+
+ return cosmo->a3_inv * p->primitives.rho;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed internal
+ * energy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param u The new internal energy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
+ struct part* restrict p, float u) {
+
+ /* conserved.energy is NOT the specific energy (u), but the total thermal
+ energy (u*m) */
+ p->conserved.energy = u * p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+ /* add the kinetic energy */
+ p->conserved.energy += 0.5f * p->conserved.mass *
+ (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
+ p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed entropy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param S The new entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_entropy(
+ struct part* restrict p, float S) {
+
+ p->conserved.energy = S * pow_gamma_minus_one(p->primitives.rho) *
+ hydro_one_over_gamma_minus_one * p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+ /* add the kinetic energy */
+ p->conserved.energy += 0.5f * p->conserved.mass *
+ (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
+ p->primitives.P = S * pow_gamma(p->primitives.rho);
+}
+
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part* p, float u_init) {
+
+ p->conserved.energy = u_init * p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+ /* add the kinetic energy */
+ p->conserved.energy += 0.5f * p->conserved.mass *
+ (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
+ p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u_init;
+}
+
+#endif /* SWIFT_GIZMO_MFV_HYDRO_H */
diff --git a/src/hydro/GizmoMFV/hydro_debug.h b/src/hydro/GizmoMFV/hydro_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..8af3f824666529efad833c3bd520ace779718449
--- /dev/null
+++ b/src/hydro/GizmoMFV/hydro_debug.h
@@ -0,0 +1,81 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * 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
+ * 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_GIZMO_MFV_HYDRO_DEBUG_H
+#define SWIFT_GIZMO_MFV_HYDRO_DEBUG_H
+
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+ const struct part* p, const struct xpart* xp) {
+ printf(
+ "x=[%.16e,%.16e,%.16e], "
+ "v=[%.3e,%.3e,%.3e], "
+ "a=[%.3e,%.3e,%.3e], "
+ "h=%.3e, "
+ "time_bin=%d, "
+ "primitives={"
+ "v=[%.3e,%.3e,%.3e], "
+ "rho=%.3e, "
+ "P=%.3e, "
+ "gradients={"
+ "rho=[%.3e,%.3e,%.3e], "
+ "v=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]], "
+ "P=[%.3e,%.3e,%.3e]}, "
+ "limiter={"
+ "rho=[%.3e,%.3e], "
+ "v=[[%.3e,%.3e],[%.3e,%.3e],[%.3e,%.3e]], "
+ "P=[%.3e,%.3e], "
+ "maxr=%.3e}}, "
+ "conserved={"
+ "momentum=[%.3e,%.3e,%.3e], "
+ "mass=%.3e, "
+ "energy=%.3e}, "
+ "geometry={"
+ "volume=%.3e, "
+ "matrix_E=[[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e],[%.3e,%.3e,%.3e]]}, "
+ "timestepvars={"
+ "vmax=%.3e},"
+ "density={"
+ "wcount_dh=%.3e, "
+ "wcount=%.3e}\n",
+ p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], p->a_hydro[0],
+ p->a_hydro[1], p->a_hydro[2], p->h, p->time_bin, p->primitives.v[0],
+ p->primitives.v[1], p->primitives.v[2], p->primitives.rho,
+ p->primitives.P, p->primitives.gradients.rho[0],
+ p->primitives.gradients.rho[1], p->primitives.gradients.rho[2],
+ p->primitives.gradients.v[0][0], p->primitives.gradients.v[0][1],
+ p->primitives.gradients.v[0][2], p->primitives.gradients.v[1][0],
+ p->primitives.gradients.v[1][1], p->primitives.gradients.v[1][2],
+ p->primitives.gradients.v[2][0], p->primitives.gradients.v[2][1],
+ p->primitives.gradients.v[2][2], p->primitives.gradients.P[0],
+ p->primitives.gradients.P[1], p->primitives.gradients.P[2],
+ p->primitives.limiter.rho[0], p->primitives.limiter.rho[1],
+ p->primitives.limiter.v[0][0], p->primitives.limiter.v[0][1],
+ p->primitives.limiter.v[1][0], p->primitives.limiter.v[1][1],
+ p->primitives.limiter.v[2][0], p->primitives.limiter.v[2][1],
+ p->primitives.limiter.P[0], p->primitives.limiter.P[1],
+ p->primitives.limiter.maxr, p->conserved.momentum[0],
+ p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.mass,
+ p->conserved.energy, p->geometry.volume, p->geometry.matrix_E[0][0],
+ p->geometry.matrix_E[0][1], p->geometry.matrix_E[0][2],
+ p->geometry.matrix_E[1][0], p->geometry.matrix_E[1][1],
+ p->geometry.matrix_E[1][2], p->geometry.matrix_E[2][0],
+ p->geometry.matrix_E[2][1], p->geometry.matrix_E[2][2],
+ p->timestepvars.vmax, p->density.wcount_dh, p->density.wcount);
+}
+
+#endif /* SWIFT_GIZMO_MFV_HYDRO_DEBUG_H */
diff --git a/src/hydro/GizmoMFV/hydro_gradients.h b/src/hydro/GizmoMFV/hydro_gradients.h
new file mode 100644
index 0000000000000000000000000000000000000000..387c263775ddfb37e4c7cb31a624ba5dc673beb2
--- /dev/null
+++ b/src/hydro/GizmoMFV/hydro_gradients.h
@@ -0,0 +1,159 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_HYDRO_GIZMO_MFV_GRADIENTS_H
+#define SWIFT_HYDRO_GIZMO_MFV_GRADIENTS_H
+
+#include "hydro_slope_limiters.h"
+#include "hydro_unphysical.h"
+#include "riemann.h"
+
+#if defined(GRADIENTS_SPH)
+
+#define HYDRO_GRADIENT_IMPLEMENTATION "SPH gradients (Price 2012)"
+#include "hydro_gradients_sph.h"
+
+#elif defined(GRADIENTS_GIZMO)
+
+#define HYDRO_GRADIENT_IMPLEMENTATION "GIZMO gradients (Hopkins 2015)"
+#include "hydro_gradients_gizmo.h"
+
+#else
+
+/* No gradients. Perfectly acceptable, but we have to provide empty functions */
+#define HYDRO_GRADIENT_IMPLEMENTATION "No gradients (first order scheme)"
+
+/**
+ * @brief Initialize gradient variables
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_init(
+ struct part *p) {}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj) {}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop: non-symmetric
+ * version
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
+ struct part *restrict pi,
+ struct part *restrict pj) {}
+
+/**
+ * @brief Finalize the gradient variables after all data have been collected
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
+ struct part *p) {}
+
+#endif
+
+/**
+ * @brief Gradients reconstruction. Is the same for all gradient types (although
+ * gradients_none does nothing, since all gradients are zero -- are they?).
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_predict(
+ struct part* restrict pi, struct part* restrict pj, float hi, float hj,
+ const float* dx, float r, const float* xij_i, float* Wi, float* Wj) {
+
+ /* perform gradient reconstruction in space and time */
+ /* Compute interface position (relative to pj, since we don't need the actual
+ * position) eqn. (8) */
+ const float xfac = hj / (hi + hj);
+ const float xij_j[3] = {xfac * dx[0], xfac * dx[1], xfac * dx[2]};
+
+ float dWi[5];
+ dWi[0] = pi->primitives.gradients.rho[0] * xij_i[0] +
+ pi->primitives.gradients.rho[1] * xij_i[1] +
+ pi->primitives.gradients.rho[2] * xij_i[2];
+ dWi[1] = pi->primitives.gradients.v[0][0] * xij_i[0] +
+ pi->primitives.gradients.v[0][1] * xij_i[1] +
+ pi->primitives.gradients.v[0][2] * xij_i[2];
+ dWi[2] = pi->primitives.gradients.v[1][0] * xij_i[0] +
+ pi->primitives.gradients.v[1][1] * xij_i[1] +
+ pi->primitives.gradients.v[1][2] * xij_i[2];
+ dWi[3] = pi->primitives.gradients.v[2][0] * xij_i[0] +
+ pi->primitives.gradients.v[2][1] * xij_i[1] +
+ pi->primitives.gradients.v[2][2] * xij_i[2];
+ dWi[4] = pi->primitives.gradients.P[0] * xij_i[0] +
+ pi->primitives.gradients.P[1] * xij_i[1] +
+ pi->primitives.gradients.P[2] * xij_i[2];
+
+ float dWj[5];
+ dWj[0] = pj->primitives.gradients.rho[0] * xij_j[0] +
+ pj->primitives.gradients.rho[1] * xij_j[1] +
+ pj->primitives.gradients.rho[2] * xij_j[2];
+ dWj[1] = pj->primitives.gradients.v[0][0] * xij_j[0] +
+ pj->primitives.gradients.v[0][1] * xij_j[1] +
+ pj->primitives.gradients.v[0][2] * xij_j[2];
+ dWj[2] = pj->primitives.gradients.v[1][0] * xij_j[0] +
+ pj->primitives.gradients.v[1][1] * xij_j[1] +
+ pj->primitives.gradients.v[1][2] * xij_j[2];
+ dWj[3] = pj->primitives.gradients.v[2][0] * xij_j[0] +
+ pj->primitives.gradients.v[2][1] * xij_j[1] +
+ pj->primitives.gradients.v[2][2] * xij_j[2];
+ dWj[4] = pj->primitives.gradients.P[0] * xij_j[0] +
+ pj->primitives.gradients.P[1] * xij_j[1] +
+ pj->primitives.gradients.P[2] * xij_j[2];
+
+ /* Apply the slope limiter at this interface */
+ hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
+
+ Wi[0] += dWi[0];
+ Wi[1] += dWi[1];
+ Wi[2] += dWi[2];
+ Wi[3] += dWi[3];
+ Wi[4] += dWi[4];
+
+ Wj[0] += dWj[0];
+ Wj[1] += dWj[1];
+ Wj[2] += dWj[2];
+ Wj[3] += dWj[3];
+ Wj[4] += dWj[4];
+
+ gizmo_check_physical_quantities("density", "pressure", Wi[0], Wi[1], Wi[2],
+ Wi[3], Wi[4]);
+ gizmo_check_physical_quantities("density", "pressure", Wj[0], Wj[1], Wj[2],
+ Wj[3], Wj[4]);
+}
+
+#endif /* SWIFT_HYDRO_GIZMO_MFV_GRADIENTS_H */
diff --git a/src/hydro/GizmoMFV/hydro_gradients_gizmo.h b/src/hydro/GizmoMFV/hydro_gradients_gizmo.h
new file mode 100644
index 0000000000000000000000000000000000000000..2592f46da9a8c118d18beba933f98f477ec5a3b2
--- /dev/null
+++ b/src/hydro/GizmoMFV/hydro_gradients_gizmo.h
@@ -0,0 +1,488 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/**
+ * @brief Initialize gradient variables
+ *
+ * @param p Particle.
+ */
+#ifndef SWIFT_GIZMO_MFV_HYDRO_GRADIENTS_H
+#define SWIFT_GIZMO_MFV_HYDRO_GRADIENTS_H
+
+__attribute__((always_inline)) INLINE static void hydro_gradients_init(
+ struct part *p) {
+
+ p->primitives.gradients.rho[0] = 0.0f;
+ p->primitives.gradients.rho[1] = 0.0f;
+ p->primitives.gradients.rho[2] = 0.0f;
+
+ p->primitives.gradients.v[0][0] = 0.0f;
+ p->primitives.gradients.v[0][1] = 0.0f;
+ p->primitives.gradients.v[0][2] = 0.0f;
+
+ p->primitives.gradients.v[1][0] = 0.0f;
+ p->primitives.gradients.v[1][1] = 0.0f;
+ p->primitives.gradients.v[1][2] = 0.0f;
+
+ p->primitives.gradients.v[2][0] = 0.0f;
+ p->primitives.gradients.v[2][1] = 0.0f;
+ p->primitives.gradients.v[2][2] = 0.0f;
+
+ p->primitives.gradients.P[0] = 0.0f;
+ p->primitives.gradients.P[1] = 0.0f;
+ p->primitives.gradients.P[2] = 0.0f;
+
+ hydro_slope_limit_cell_init(p);
+}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj) {
+
+ const float r_inv = 1.f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ float wi, wj, wi_dx, wj_dx;
+ float Bi[3][3];
+ float Bj[3][3];
+ float Wi[5], Wj[5];
+
+ /* Initialize local variables */
+ for (int k = 0; k < 3; k++) {
+ for (int l = 0; l < 3; l++) {
+ Bi[k][l] = pi->geometry.matrix_E[k][l];
+ Bj[k][l] = pj->geometry.matrix_E[k][l];
+ }
+ }
+ Wi[0] = pi->primitives.rho;
+ Wi[1] = pi->primitives.v[0];
+ Wi[2] = pi->primitives.v[1];
+ Wi[3] = pi->primitives.v[2];
+ Wi[4] = pi->primitives.P;
+ Wj[0] = pj->primitives.rho;
+ Wj[1] = pj->primitives.v[0];
+ Wj[2] = pj->primitives.v[1];
+ Wj[3] = pj->primitives.v[2];
+ Wj[4] = pj->primitives.P;
+
+ /* Compute kernel of pi. */
+ const float hi_inv = 1.f / hi;
+ const float xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ if (pi->density.wcorr > const_gizmo_min_wcorr) {
+ /* Compute gradients for pi */
+ /* there is a sign difference w.r.t. eqn. (6) because of the inverse
+ * definition of dx */
+ pi->primitives.gradients.rho[0] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.rho[1] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.rho[2] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->primitives.gradients.v[0][0] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[0][1] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[0][2] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ pi->primitives.gradients.v[1][0] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[1][1] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[1][2] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ pi->primitives.gradients.v[2][0] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[2][1] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[2][2] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->primitives.gradients.P[0] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.P[1] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.P[2] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ } else {
+ /* The gradient matrix was not well-behaved, switch to SPH gradients */
+
+ pi->primitives.gradients.rho[0] -=
+ wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+ pi->primitives.gradients.rho[1] -=
+ wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+ pi->primitives.gradients.rho[2] -=
+ wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+
+ pi->primitives.gradients.v[0][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+ pi->primitives.gradients.v[0][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+ pi->primitives.gradients.v[0][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+
+ pi->primitives.gradients.v[1][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pi->primitives.gradients.v[1][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pi->primitives.gradients.v[1][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+
+ pi->primitives.gradients.v[2][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+ pi->primitives.gradients.v[2][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+ pi->primitives.gradients.v[2][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+
+ pi->primitives.gradients.P[0] -=
+ wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ pi->primitives.gradients.P[1] -=
+ wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ pi->primitives.gradients.P[2] -=
+ wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ }
+
+ hydro_slope_limit_cell_collect(pi, pj, r);
+
+ /* Compute kernel of pj. */
+ const float hj_inv = 1.f / hj;
+ const float xj = r * hj_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ if (pj->density.wcorr > const_gizmo_min_wcorr) {
+ /* Compute gradients for pj */
+ /* there is no sign difference w.r.t. eqn. (6) because dx is now what we
+ * want
+ * it to be */
+ pj->primitives.gradients.rho[0] +=
+ (Wi[0] - Wj[0]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.rho[1] +=
+ (Wi[0] - Wj[0]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.rho[2] +=
+ (Wi[0] - Wj[0]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+ pj->primitives.gradients.v[0][0] +=
+ (Wi[1] - Wj[1]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.v[0][1] +=
+ (Wi[1] - Wj[1]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.v[0][2] +=
+ (Wi[1] - Wj[1]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+ pj->primitives.gradients.v[1][0] +=
+ (Wi[2] - Wj[2]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.v[1][1] +=
+ (Wi[2] - Wj[2]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.v[1][2] +=
+ (Wi[2] - Wj[2]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+ pj->primitives.gradients.v[2][0] +=
+ (Wi[3] - Wj[3]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.v[2][1] +=
+ (Wi[3] - Wj[3]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.v[2][2] +=
+ (Wi[3] - Wj[3]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+ pj->primitives.gradients.P[0] +=
+ (Wi[4] - Wj[4]) * wj *
+ (Bj[0][0] * dx[0] + Bj[0][1] * dx[1] + Bj[0][2] * dx[2]);
+ pj->primitives.gradients.P[1] +=
+ (Wi[4] - Wj[4]) * wj *
+ (Bj[1][0] * dx[0] + Bj[1][1] * dx[1] + Bj[1][2] * dx[2]);
+ pj->primitives.gradients.P[2] +=
+ (Wi[4] - Wj[4]) * wj *
+ (Bj[2][0] * dx[0] + Bj[2][1] * dx[1] + Bj[2][2] * dx[2]);
+
+ } else {
+ /* SPH gradients */
+
+ pj->primitives.gradients.rho[0] -=
+ wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+ pj->primitives.gradients.rho[1] -=
+ wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+ pj->primitives.gradients.rho[2] -=
+ wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+
+ pj->primitives.gradients.v[0][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+ pj->primitives.gradients.v[0][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+ pj->primitives.gradients.v[0][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+
+ pj->primitives.gradients.v[1][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pj->primitives.gradients.v[1][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pj->primitives.gradients.v[1][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pj->primitives.gradients.v[2][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+ pj->primitives.gradients.v[2][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+ pj->primitives.gradients.v[2][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+
+ pj->primitives.gradients.P[0] -=
+ wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ pj->primitives.gradients.P[1] -=
+ wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ pj->primitives.gradients.P[2] -=
+ wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ }
+
+ hydro_slope_limit_cell_collect(pj, pi, r);
+}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
+ struct part *restrict pi,
+ struct part *restrict pj) {
+
+ const float r_inv = 1.f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ float Bi[3][3];
+ float Wi[5], Wj[5];
+
+ /* Initialize local variables */
+ for (int k = 0; k < 3; k++) {
+ for (int l = 0; l < 3; l++) {
+ Bi[k][l] = pi->geometry.matrix_E[k][l];
+ }
+ }
+ Wi[0] = pi->primitives.rho;
+ Wi[1] = pi->primitives.v[0];
+ Wi[2] = pi->primitives.v[1];
+ Wi[3] = pi->primitives.v[2];
+ Wi[4] = pi->primitives.P;
+ Wj[0] = pj->primitives.rho;
+ Wj[1] = pj->primitives.v[0];
+ Wj[2] = pj->primitives.v[1];
+ Wj[3] = pj->primitives.v[2];
+ Wj[4] = pj->primitives.P;
+
+ /* Compute kernel of pi. */
+ float wi, wi_dx;
+ const float hi_inv = 1.f / hi;
+ const float xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ if (pi->density.wcorr > const_gizmo_min_wcorr) {
+ /* Compute gradients for pi */
+ /* there is a sign difference w.r.t. eqn. (6) because of the inverse
+ * definition of dx */
+ pi->primitives.gradients.rho[0] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.rho[1] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.rho[2] +=
+ (Wi[0] - Wj[0]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->primitives.gradients.v[0][0] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[0][1] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[0][2] +=
+ (Wi[1] - Wj[1]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ pi->primitives.gradients.v[1][0] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[1][1] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[1][2] +=
+ (Wi[2] - Wj[2]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+ pi->primitives.gradients.v[2][0] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.v[2][1] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.v[2][2] +=
+ (Wi[3] - Wj[3]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ pi->primitives.gradients.P[0] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[0][0] * dx[0] + Bi[0][1] * dx[1] + Bi[0][2] * dx[2]);
+ pi->primitives.gradients.P[1] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[1][0] * dx[0] + Bi[1][1] * dx[1] + Bi[1][2] * dx[2]);
+ pi->primitives.gradients.P[2] +=
+ (Wi[4] - Wj[4]) * wi *
+ (Bi[2][0] * dx[0] + Bi[2][1] * dx[1] + Bi[2][2] * dx[2]);
+
+ } else {
+ /* Gradient matrix is not well-behaved, switch to SPH gradients */
+
+ pi->primitives.gradients.rho[0] -=
+ wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+ pi->primitives.gradients.rho[1] -=
+ wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+ pi->primitives.gradients.rho[2] -=
+ wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+
+ pi->primitives.gradients.v[0][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+ pi->primitives.gradients.v[0][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+ pi->primitives.gradients.v[0][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+ pi->primitives.gradients.v[1][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pi->primitives.gradients.v[1][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pi->primitives.gradients.v[1][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+
+ pi->primitives.gradients.v[2][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+ pi->primitives.gradients.v[2][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+ pi->primitives.gradients.v[2][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+
+ pi->primitives.gradients.P[0] -=
+ wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ pi->primitives.gradients.P[1] -=
+ wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ pi->primitives.gradients.P[2] -=
+ wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ }
+
+ hydro_slope_limit_cell_collect(pi, pj, r);
+}
+
+/**
+ * @brief Finalize the gradient variables after all data have been collected
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
+ struct part *p) {
+
+ /* add kernel normalization to gradients */
+ const float volume = p->geometry.volume;
+ const float h = p->h;
+ const float h_inv = 1.0f / h;
+ const float ihdim = pow_dimension(h_inv);
+ const float ihdimp1 = pow_dimension_plus_one(h_inv);
+
+ if (p->density.wcorr > const_gizmo_min_wcorr) {
+ p->primitives.gradients.rho[0] *= ihdim;
+ p->primitives.gradients.rho[1] *= ihdim;
+ p->primitives.gradients.rho[2] *= ihdim;
+
+ p->primitives.gradients.v[0][0] *= ihdim;
+ p->primitives.gradients.v[0][1] *= ihdim;
+ p->primitives.gradients.v[0][2] *= ihdim;
+ p->primitives.gradients.v[1][0] *= ihdim;
+ p->primitives.gradients.v[1][1] *= ihdim;
+ p->primitives.gradients.v[1][2] *= ihdim;
+ p->primitives.gradients.v[2][0] *= ihdim;
+ p->primitives.gradients.v[2][1] *= ihdim;
+ p->primitives.gradients.v[2][2] *= ihdim;
+
+ p->primitives.gradients.P[0] *= ihdim;
+ p->primitives.gradients.P[1] *= ihdim;
+ p->primitives.gradients.P[2] *= ihdim;
+
+ } else {
+
+ /* finalize gradients by multiplying with volume */
+ p->primitives.gradients.rho[0] *= ihdimp1 * volume;
+ p->primitives.gradients.rho[1] *= ihdimp1 * volume;
+ p->primitives.gradients.rho[2] *= ihdimp1 * volume;
+
+ p->primitives.gradients.v[0][0] *= ihdimp1 * volume;
+ p->primitives.gradients.v[0][1] *= ihdimp1 * volume;
+ p->primitives.gradients.v[0][2] *= ihdimp1 * volume;
+
+ p->primitives.gradients.v[1][0] *= ihdimp1 * volume;
+ p->primitives.gradients.v[1][1] *= ihdimp1 * volume;
+ p->primitives.gradients.v[1][2] *= ihdimp1 * volume;
+ p->primitives.gradients.v[2][0] *= ihdimp1 * volume;
+ p->primitives.gradients.v[2][1] *= ihdimp1 * volume;
+ p->primitives.gradients.v[2][2] *= ihdimp1 * volume;
+
+ p->primitives.gradients.P[0] *= ihdimp1 * volume;
+ p->primitives.gradients.P[1] *= ihdimp1 * volume;
+ p->primitives.gradients.P[2] *= ihdimp1 * volume;
+ }
+
+ hydro_slope_limit_cell(p);
+}
+
+#endif /* SWIFT_GIZMO_MFV_HYDRO_GRADIENTS_H */
diff --git a/src/hydro/GizmoMFV/hydro_gradients_sph.h b/src/hydro/GizmoMFV/hydro_gradients_sph.h
new file mode 100644
index 0000000000000000000000000000000000000000..7b2b89ae688622ff017e4188f9e8fd2c7cee9d7a
--- /dev/null
+++ b/src/hydro/GizmoMFV/hydro_gradients_sph.h
@@ -0,0 +1,256 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/**
+ * @brief Initialize gradient variables
+ *
+ * @param p Particle.
+ */
+#ifndef SWIFT_GIZMO_MFV_HYDRO_SPH_GRADIENTS_H
+#define SWIFT_GIZMO_MFV_HYDRO_SPH_GRADIENTS_H
+
+__attribute__((always_inline)) INLINE static void hydro_gradients_init(
+ struct part *p) {
+
+ p->primitives.gradients.rho[0] = 0.0f;
+ p->primitives.gradients.rho[1] = 0.0f;
+ p->primitives.gradients.rho[2] = 0.0f;
+
+ p->primitives.gradients.v[0][0] = 0.0f;
+ p->primitives.gradients.v[0][1] = 0.0f;
+ p->primitives.gradients.v[0][2] = 0.0f;
+
+ p->primitives.gradients.v[1][0] = 0.0f;
+ p->primitives.gradients.v[1][1] = 0.0f;
+ p->primitives.gradients.v[1][2] = 0.0f;
+ p->primitives.gradients.v[2][0] = 0.0f;
+ p->primitives.gradients.v[2][1] = 0.0f;
+ p->primitives.gradients.v[2][2] = 0.0f;
+
+ p->primitives.gradients.P[0] = 0.0f;
+ p->primitives.gradients.P[1] = 0.0f;
+ p->primitives.gradients.P[2] = 0.0f;
+
+ hydro_slope_limit_cell_init(p);
+}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
+ float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+ struct part *restrict pj) {
+
+ const float r_inv = 1.f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ float wi, wi_dx;
+ const float hi_inv = 1.0f / hi;
+ const float xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ /* very basic gradient estimate */
+ pi->primitives.gradients.rho[0] -=
+ wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+ pi->primitives.gradients.rho[1] -=
+ wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+ pi->primitives.gradients.rho[2] -=
+ wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+
+ pi->primitives.gradients.v[0][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+ pi->primitives.gradients.v[0][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+ pi->primitives.gradients.v[0][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+
+ pi->primitives.gradients.v[1][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pi->primitives.gradients.v[1][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pi->primitives.gradients.v[1][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+
+ pi->primitives.gradients.v[2][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+ pi->primitives.gradients.v[2][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+ pi->primitives.gradients.v[2][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+
+ pi->primitives.gradients.P[0] -=
+ wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ pi->primitives.gradients.P[1] -=
+ wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ pi->primitives.gradients.P[2] -=
+ wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+
+ hydro_slope_limit_cell_collect(pi, pj, r);
+
+ float wj, wj_dx;
+ const float hj_inv = 1.0f / hj;
+ const float xj = r * hj_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+
+ /* signs are the same as before, since we swap i and j twice */
+ pj->primitives.gradients.rho[0] -=
+ wj_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+ pj->primitives.gradients.rho[1] -=
+ wj_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+ pj->primitives.gradients.rho[2] -=
+ wj_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+
+ pj->primitives.gradients.v[0][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+ pj->primitives.gradients.v[0][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+ pj->primitives.gradients.v[0][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+
+ pj->primitives.gradients.v[1][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pj->primitives.gradients.v[1][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pj->primitives.gradients.v[1][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pj->primitives.gradients.v[2][0] -=
+ wj_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+ pj->primitives.gradients.v[2][1] -=
+ wj_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+ pj->primitives.gradients.v[2][2] -=
+ wj_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+
+ pj->primitives.gradients.P[0] -=
+ wj_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ pj->primitives.gradients.P[1] -=
+ wj_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ pj->primitives.gradients.P[2] -=
+ wj_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+
+ hydro_slope_limit_cell_collect(pj, pi, r);
+}
+
+/**
+ * @brief Gradient calculations done during the neighbour loop: non-symmetric
+ * version
+ *
+ * @param r2 Squared distance between the two particles.
+ * @param dx Distance vector (pi->x - pj->x).
+ * @param hi Smoothing length of particle i.
+ * @param hj Smoothing length of particle j.
+ * @param pi Particle i.
+ * @param pj Particle j.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_gradients_nonsym_collect(float r2, const float *dx, float hi, float hj,
+ struct part *restrict pi,
+ struct part *restrict pj) {
+
+ const float r_inv = 1.f / sqrtf(r2);
+ const float r = r2 * r_inv;
+
+ float wi, wi_dx;
+ const float hi_inv = 1.0f / hi;
+ const float xi = r * hi_inv;
+ kernel_deval(xi, &wi, &wi_dx);
+
+ /* very basic gradient estimate */
+ pi->primitives.gradients.rho[0] -=
+ wi_dx * dx[0] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+ pi->primitives.gradients.rho[1] -=
+ wi_dx * dx[1] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+ pi->primitives.gradients.rho[2] -=
+ wi_dx * dx[2] * (pi->primitives.rho - pj->primitives.rho) * r_inv;
+
+ pi->primitives.gradients.v[0][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+ pi->primitives.gradients.v[0][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+ pi->primitives.gradients.v[0][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[0] - pj->primitives.v[0]) * r_inv;
+
+ pi->primitives.gradients.v[1][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pi->primitives.gradients.v[1][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+ pi->primitives.gradients.v[1][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[1] - pj->primitives.v[1]) * r_inv;
+
+ pi->primitives.gradients.v[2][0] -=
+ wi_dx * dx[0] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+ pi->primitives.gradients.v[2][1] -=
+ wi_dx * dx[1] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+ pi->primitives.gradients.v[2][2] -=
+ wi_dx * dx[2] * (pi->primitives.v[2] - pj->primitives.v[2]) * r_inv;
+
+ pi->primitives.gradients.P[0] -=
+ wi_dx * dx[0] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ pi->primitives.gradients.P[1] -=
+ wi_dx * dx[1] * (pi->primitives.P - pj->primitives.P) * r_inv;
+ pi->primitives.gradients.P[2] -=
+ wi_dx * dx[2] * (pi->primitives.P - pj->primitives.P) * r_inv;
+
+ hydro_slope_limit_cell_collect(pi, pj, r);
+}
+
+/**
+ * @brief Finalize the gradient variables after all data have been collected
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_gradients_finalize(
+ struct part *p) {
+
+ const float h = p->h;
+ const float ih = 1.0f / h;
+ const float ihdimp1 = pow_dimension_plus_one(ih);
+ const float volume = p->geometry.volume;
+
+ /* finalize gradients by multiplying with volume */
+ p->primitives.gradients.rho[0] *= ihdimp1 * volume;
+ p->primitives.gradients.rho[1] *= ihdimp1 * volume;
+ p->primitives.gradients.rho[2] *= ihdimp1 * volume;
+
+ p->primitives.gradients.v[0][0] *= ihdimp1 * volume;
+ p->primitives.gradients.v[0][1] *= ihdimp1 * volume;
+ p->primitives.gradients.v[0][2] *= ihdimp1 * volume;
+
+ p->primitives.gradients.v[1][0] *= ihdimp1 * volume;
+ p->primitives.gradients.v[1][1] *= ihdimp1 * volume;
+ p->primitives.gradients.v[1][2] *= ihdimp1 * volume;
+
+ p->primitives.gradients.v[2][0] *= ihdimp1 * volume;
+ p->primitives.gradients.v[2][1] *= ihdimp1 * volume;
+ p->primitives.gradients.v[2][2] *= ihdimp1 * volume;
+
+ p->primitives.gradients.P[0] *= ihdimp1 * volume;
+ p->primitives.gradients.P[1] *= ihdimp1 * volume;
+ p->primitives.gradients.P[2] *= ihdimp1 * volume;
+
+ hydro_slope_limit_cell(p);
+}
+
+#endif /* SWIFT_GIZMO_MFV_HYDRO_SPH_GRADIENTS_H */
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/GizmoMFV/hydro_iact.h
similarity index 99%
rename from src/hydro/Gizmo/hydro_iact.h
rename to src/hydro/GizmoMFV/hydro_iact.h
index 0c888e4159c6da6eb07eb0afbf758efce1c697ee..bb835094acd285b109383c1f5d04a6f5e2d936df 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/GizmoMFV/hydro_iact.h
@@ -18,8 +18,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_GIZMO_HYDRO_IACT_H
-#define SWIFT_GIZMO_HYDRO_IACT_H
+#ifndef SWIFT_GIZMO_MFV_HYDRO_IACT_H
+#define SWIFT_GIZMO_MFV_HYDRO_IACT_H
#include "adiabatic_index.h"
#include "hydro_gradients.h"
@@ -514,4 +514,4 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0, a, H);
}
-#endif /* SWIFT_GIZMO_HYDRO_IACT_H */
+#endif /* SWIFT_GIZMO_MFV_HYDRO_IACT_H */
diff --git a/src/hydro/GizmoMFV/hydro_io.h b/src/hydro/GizmoMFV/hydro_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..c1b151230f3198a30d6696e36a2704156804fdce
--- /dev/null
+++ b/src/hydro/GizmoMFV/hydro_io.h
@@ -0,0 +1,244 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_MFV_HYDRO_IO_H
+#define SWIFT_GIZMO_MFV_HYDRO_IO_H
+
+#include "adiabatic_index.h"
+#include "hydro.h"
+#include "hydro_gradients.h"
+#include "hydro_slope_limiters.h"
+#include "io_properties.h"
+#include "riemann.h"
+
+/* Set the description of the particle movement. */
+#if defined(GIZMO_FIX_PARTICLES)
+#define GIZMO_PARTICLE_MOVEMENT "Fixed particles."
+#else
+#define GIZMO_PARTICLE_MOVEMENT "Particles move with flow velocity."
+#endif
+
+/**
+ * @brief Specifies which particle fields to read from a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+void hydro_read_particles(struct part* parts, struct io_props* list,
+ int* num_fields) {
+
+ *num_fields = 8;
+
+ /* List what we want to read */
+ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+ 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("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+ parts, conserved.mass);
+ list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_LENGTH, parts, h);
+ list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS, parts,
+ conserved.energy);
+ list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+ UNIT_CONV_NO_UNITS, parts, id);
+ 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, primitives.rho);
+}
+
+/**
+ * @brief Get the internal energy of a particle
+ *
+ * @param e #engine.
+ * @param p Particle.
+ * @param ret (return) Internal energy of the particle
+ */
+void convert_u(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ ret[0] = hydro_get_comoving_internal_energy(p);
+}
+
+/**
+ * @brief Get the entropic function of a particle
+ *
+ * @param e #engine.
+ * @param p Particle.
+ * @param ret (return) Entropic function of the particle
+ */
+void convert_A(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+ ret[0] = hydro_get_comoving_entropy(p);
+}
+
+/**
+ * @brief Get the total energy of a particle
+ *
+ * @param e #engine.
+ * @param p Particle.
+ * @return Total energy of the particle
+ */
+void convert_Etot(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+#ifdef GIZMO_TOTAL_ENERGY
+ ret[0] = p->conserved.energy;
+#else
+ float momentum2;
+
+ momentum2 = p->conserved.momentum[0] * p->conserved.momentum[0] +
+ p->conserved.momentum[1] * p->conserved.momentum[1] +
+ p->conserved.momentum[2] * p->conserved.momentum[2];
+
+ ret[0] = p->conserved.energy + 0.5f * momentum2 / p->conserved.mass;
+#endif
+}
+
+void convert_part_pos(const struct engine* e, const struct part* p,
+ const struct xpart* xp, double* ret) {
+
+ if (e->s->periodic) {
+ ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
+ ret[1] = box_wrap(p->x[1], 0.0, e->s->dim[1]);
+ ret[2] = box_wrap(p->x[2], 0.0, e->s->dim[2]);
+ } else {
+ ret[0] = p->x[0];
+ ret[1] = p->x[1];
+ ret[2] = p->x[2];
+ }
+}
+
+void convert_part_vel(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a2_inv;
+ ret[1] *= cosmo->a2_inv;
+ ret[2] *= cosmo->a2_inv;
+}
+
+void convert_part_potential(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ if (p->gpart != NULL)
+ ret[0] = gravity_get_comoving_potential(p->gpart);
+ else
+ ret[0] = 0.f;
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
+ struct io_props* list, int* num_fields) {
+
+ *num_fields = 11;
+
+ /* List what we want to write */
+ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+ UNIT_CONV_LENGTH, parts, xparts,
+ convert_part_pos);
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
+
+ list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts,
+ conserved.mass);
+ list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+ parts, h);
+ list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS,
+ parts, xparts, convert_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,
+ primitives.rho);
+ list[7] = io_make_output_field_convert_part(
+ "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, xparts, convert_A);
+ list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
+ parts, primitives.P);
+ list[9] = io_make_output_field_convert_part(
+ "TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, xparts, convert_Etot);
+
+ list[10] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+ UNIT_CONV_POTENTIAL, parts,
+ xparts, convert_part_potential);
+}
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grpsph The HDF5 group in which to write
+ */
+void hydro_write_flavour(hid_t h_grpsph) {
+ /* Gradient information */
+ io_write_attribute_s(h_grpsph, "Gradient reconstruction model",
+ HYDRO_GRADIENT_IMPLEMENTATION);
+
+ /* Slope limiter information */
+ io_write_attribute_s(h_grpsph, "Cell wide slope limiter model",
+ HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION);
+ io_write_attribute_s(h_grpsph, "Piecewise slope limiter model",
+ HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION);
+
+ /* Riemann solver information */
+ io_write_attribute_s(h_grpsph, "Riemann solver type",
+ RIEMANN_SOLVER_IMPLEMENTATION);
+
+ /* Particle movement information */
+ io_write_attribute_s(h_grpsph, "Particle movement", GIZMO_PARTICLE_MOVEMENT);
+}
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+int writeEntropyFlag() { return 0; }
+
+#endif /* SWIFT_GIZMO_MFV_HYDRO_IO_H */
diff --git a/src/hydro/GizmoMFV/hydro_part.h b/src/hydro/GizmoMFV/hydro_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..6248ddb11daf39a65be9a57fe51e40386ecda50b
--- /dev/null
+++ b/src/hydro/GizmoMFV/hydro_part.h
@@ -0,0 +1,213 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2014 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ *
+ * 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_GIZMO_MFV_HYDRO_PART_H
+#define SWIFT_GIZMO_MFV_HYDRO_PART_H
+
+#include "chemistry_struct.h"
+#include "cooling_struct.h"
+
+/* Extra particle data not needed during the computation. */
+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. */
+ float x_diff_sort[3];
+
+ /* Velocity at the last full step. */
+ float v_full[3];
+
+ /* Gravitational acceleration at the last full step. */
+ float a_grav[3];
+
+ /* Additional data used to record cooling information */
+ struct cooling_xpart_data cooling_data;
+
+} SWIFT_STRUCT_ALIGN;
+
+/* Data of a single particle. */
+struct part {
+
+ /* Particle ID. */
+ long long id;
+
+ /* Associated gravitas. */
+ struct gpart *gpart;
+
+ /* Particle position. */
+ double x[3];
+
+ /* Particle predicted velocity. */
+ float v[3];
+
+ /* Particle acceleration. */
+ float a_hydro[3];
+
+ /* Particle smoothing length. */
+ float h;
+
+ /* The primitive hydrodynamical variables. */
+ struct {
+
+ /* Density. */
+ float rho;
+
+ /* Fluid velocity. */
+ float v[3];
+
+ /* Pressure. */
+ float P;
+
+ /* Gradients of the primitive variables. */
+ struct {
+
+ /* Density gradients. */
+ float rho[3];
+
+ /* Fluid velocity gradients. */
+ float v[3][3];
+
+ /* Pressure gradients. */
+ float P[3];
+
+ } gradients;
+
+ /* Quantities needed by the slope limiter. */
+ struct {
+
+ /* Extreme values of the density among the neighbours. */
+ float rho[2];
+
+ /* Extreme values of the fluid velocity among the neighbours. */
+ float v[3][2];
+
+ /* Extreme values of the pressure among the neighbours. */
+ float P[2];
+
+ /* Maximal distance to all neighbouring faces. */
+ float maxr;
+
+ } limiter;
+
+ } primitives;
+
+ /* The conserved hydrodynamical variables. */
+ struct {
+
+ /* Fluid mass */
+ float mass;
+
+ /* Fluid momentum. */
+ float momentum[3];
+
+ /* Fluid thermal energy (not per unit mass!). */
+ float energy;
+
+ /* Fluxes. */
+ struct {
+
+ /* Mass flux. */
+ float mass;
+
+ /* Momentum flux. */
+ float momentum[3];
+
+ /* Energy flux. */
+ float energy;
+
+ } flux;
+
+ } conserved;
+
+ /* Geometrical quantities used for hydro. */
+ struct {
+
+ /* Volume of the particle. */
+ float volume;
+
+ /* Geometrical shear matrix used to calculate second order accurate
+ gradients */
+ float matrix_E[3][3];
+
+ /* Centroid of the "cell". */
+ float centroid[3];
+
+ } geometry;
+
+ /* Variables used for timestep calculation. */
+ struct {
+
+ /* Maximum signal velocity among all the neighbours of the particle. The
+ * signal velocity encodes information about the relative fluid velocities
+ * AND particle velocities of the neighbour and this particle, as well as
+ * the sound speed of both particles. */
+ float vmax;
+
+ } timestepvars;
+
+ /* Quantities used during the volume (=density) loop. */
+ struct {
+
+ /* Derivative of particle number density. */
+ float wcount_dh;
+
+ /* Particle number density. */
+ float wcount;
+
+ /* Correction factor for wcount. */
+ float wcorr;
+
+ } density;
+
+ /* Quantities used during the force loop. */
+ struct {
+
+ /* Needed to drift the primitive variables. */
+ float h_dt;
+
+ } force;
+
+ /* Specific stuff for the gravity-hydro coupling. */
+ struct {
+
+ /* Current value of the mass flux vector. */
+ float mflux[3];
+
+ } gravity;
+
+ /* Chemistry information */
+ struct chemistry_part_data chemistry_data;
+
+ /* Time-step length */
+ timebin_t time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Time of the last drift */
+ integertime_t ti_drift;
+
+ /* Time of the last kick */
+ integertime_t ti_kick;
+
+#endif
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_GIZMO_MFV_HYDRO_PART_H */
diff --git a/src/hydro/GizmoMFV/hydro_slope_limiters.h b/src/hydro/GizmoMFV/hydro_slope_limiters.h
new file mode 100644
index 0000000000000000000000000000000000000000..78f2785cdae5dc2334d37e3924dd5b259cca8c05
--- /dev/null
+++ b/src/hydro/GizmoMFV/hydro_slope_limiters.h
@@ -0,0 +1,94 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_HYDRO_SLOPE_LIMITERS_H
+#define SWIFT_HYDRO_SLOPE_LIMITERS_H
+
+#include "dimension.h"
+#include "kernel_hydro.h"
+
+#ifdef SLOPE_LIMITER_PER_FACE
+
+#define HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION \
+ "GIZMO piecewise slope limiter (Hopkins 2015)"
+#include "hydro_slope_limiters_face.h"
+
+#else
+
+#define HYDRO_SLOPE_LIMITER_FACE_IMPLEMENTATION "No piecewise slope limiter"
+
+/**
+ * @brief Slope limit the slopes at the interface between two particles
+ *
+ * @param Wi Hydrodynamic variables of particle i.
+ * @param Wj Hydrodynamic variables of particle j.
+ * @param dWi Difference between the hydrodynamic variables of particle i at the
+ * position of particle i and at the interface position.
+ * @param dWj Difference between the hydrodynamic variables of particle j at the
+ * position of particle j and at the interface position.
+ * @param xij_i Relative position vector of the interface w.r.t. particle i.
+ * @param xij_j Relative position vector of the interface w.r.t. partilce j.
+ * @param r Distance between particle i and particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
+ float *Wi, float *Wj, float *dWi, float *dWj, float *xij_i, float *xij_j,
+ float r) {}
+
+#endif
+
+#ifdef SLOPE_LIMITER_CELL_WIDE
+
+#define HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION \
+ "Cell wide slope limiter (Springel 2010)"
+#include "hydro_slope_limiters_cell.h"
+
+#else
+
+#define HYDRO_SLOPE_LIMITER_CELL_IMPLEMENTATION "No cell wide slope limiter"
+
+/**
+ * @brief Initialize variables for the cell wide slope limiter
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
+ struct part *p) {}
+
+/**
+ * @brief Collect information for the cell wide slope limiter during the
+ * neighbour loop
+ *
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param r Distance between particle i and particle j.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_slope_limit_cell_collect(struct part *pi, struct part *pj, float r) {}
+
+/**
+ * @brief Slope limit cell gradients
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
+ struct part *p) {}
+
+#endif
+
+#endif /* SWIFT_HYDRO_SLOPE_LIMITERS_H */
diff --git a/src/hydro/GizmoMFV/hydro_slope_limiters_cell.h b/src/hydro/GizmoMFV/hydro_slope_limiters_cell.h
new file mode 100644
index 0000000000000000000000000000000000000000..130a87def95a3198f0871ce485f0a9d377a96cd5
--- /dev/null
+++ b/src/hydro/GizmoMFV/hydro_slope_limiters_cell.h
@@ -0,0 +1,186 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GIZMO_MFV_SLOPE_LIMITER_CELL_H
+#define SWIFT_GIZMO_MFV_SLOPE_LIMITER_CELL_H
+
+#include <float.h>
+
+/**
+ * @brief Initialize variables for the cell wide slope limiter
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell_init(
+ struct part* p) {
+
+ p->primitives.limiter.rho[0] = FLT_MAX;
+ p->primitives.limiter.rho[1] = -FLT_MAX;
+ p->primitives.limiter.v[0][0] = FLT_MAX;
+ p->primitives.limiter.v[0][1] = -FLT_MAX;
+ p->primitives.limiter.v[1][0] = FLT_MAX;
+ p->primitives.limiter.v[1][1] = -FLT_MAX;
+ p->primitives.limiter.v[2][0] = FLT_MAX;
+ p->primitives.limiter.v[2][1] = -FLT_MAX;
+ p->primitives.limiter.P[0] = FLT_MAX;
+ p->primitives.limiter.P[1] = -FLT_MAX;
+
+ p->primitives.limiter.maxr = -FLT_MAX;
+}
+
+/**
+ * @brief Collect information for the cell wide slope limiter during the
+ * neighbour loop
+ *
+ * @param pi Particle i.
+ * @param pj Particle j.
+ * @param r Distance between particle i and particle j.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_slope_limit_cell_collect(struct part* pi, struct part* pj, float r) {
+
+ /* basic slope limiter: collect the maximal and the minimal value for the
+ * primitive variables among the ngbs */
+ pi->primitives.limiter.rho[0] =
+ min(pj->primitives.rho, pi->primitives.limiter.rho[0]);
+ pi->primitives.limiter.rho[1] =
+ max(pj->primitives.rho, pi->primitives.limiter.rho[1]);
+
+ pi->primitives.limiter.v[0][0] =
+ min(pj->primitives.v[0], pi->primitives.limiter.v[0][0]);
+ pi->primitives.limiter.v[0][1] =
+ max(pj->primitives.v[0], pi->primitives.limiter.v[0][1]);
+ pi->primitives.limiter.v[1][0] =
+ min(pj->primitives.v[1], pi->primitives.limiter.v[1][0]);
+ pi->primitives.limiter.v[1][1] =
+ max(pj->primitives.v[1], pi->primitives.limiter.v[1][1]);
+ pi->primitives.limiter.v[2][0] =
+ min(pj->primitives.v[2], pi->primitives.limiter.v[2][0]);
+ pi->primitives.limiter.v[2][1] =
+ max(pj->primitives.v[2], pi->primitives.limiter.v[2][1]);
+
+ pi->primitives.limiter.P[0] =
+ min(pj->primitives.P, pi->primitives.limiter.P[0]);
+ pi->primitives.limiter.P[1] =
+ max(pj->primitives.P, pi->primitives.limiter.P[1]);
+
+ pi->primitives.limiter.maxr = max(r, pi->primitives.limiter.maxr);
+}
+
+/**
+ * @brief Slope limit cell gradients
+ *
+ * @param p Particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_cell(
+ struct part* p) {
+
+ float gradtrue, gradrho[3], gradv[3][3], gradP[3];
+
+ gradrho[0] = p->primitives.gradients.rho[0];
+ gradrho[1] = p->primitives.gradients.rho[1];
+ gradrho[2] = p->primitives.gradients.rho[2];
+
+ gradv[0][0] = p->primitives.gradients.v[0][0];
+ gradv[0][1] = p->primitives.gradients.v[0][1];
+ gradv[0][2] = p->primitives.gradients.v[0][2];
+
+ gradv[1][0] = p->primitives.gradients.v[1][0];
+ gradv[1][1] = p->primitives.gradients.v[1][1];
+ gradv[1][2] = p->primitives.gradients.v[1][2];
+
+ gradv[2][0] = p->primitives.gradients.v[2][0];
+ gradv[2][1] = p->primitives.gradients.v[2][1];
+ gradv[2][2] = p->primitives.gradients.v[2][2];
+
+ gradP[0] = p->primitives.gradients.P[0];
+ gradP[1] = p->primitives.gradients.P[1];
+ gradP[2] = p->primitives.gradients.P[2];
+
+ gradtrue = sqrtf(gradrho[0] * gradrho[0] + gradrho[1] * gradrho[1] +
+ gradrho[2] * gradrho[2]);
+ if (gradtrue) {
+ gradtrue *= p->primitives.limiter.maxr;
+ const float gradmax = p->primitives.limiter.rho[1] - p->primitives.rho;
+ const float gradmin = p->primitives.rho - p->primitives.limiter.rho[0];
+ const float gradtrue_inv = 1.f / gradtrue;
+ const float alpha =
+ min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
+ p->primitives.gradients.rho[0] *= alpha;
+ p->primitives.gradients.rho[1] *= alpha;
+ p->primitives.gradients.rho[2] *= alpha;
+ }
+
+ gradtrue = sqrtf(gradv[0][0] * gradv[0][0] + gradv[0][1] * gradv[0][1] +
+ gradv[0][2] * gradv[0][2]);
+ if (gradtrue) {
+ gradtrue *= p->primitives.limiter.maxr;
+ const float gradmax = p->primitives.limiter.v[0][1] - p->primitives.v[0];
+ const float gradmin = p->primitives.v[0] - p->primitives.limiter.v[0][0];
+ const float gradtrue_inv = 1.f / gradtrue;
+ const float alpha =
+ min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
+ p->primitives.gradients.v[0][0] *= alpha;
+ p->primitives.gradients.v[0][1] *= alpha;
+ p->primitives.gradients.v[0][2] *= alpha;
+ }
+
+ gradtrue = sqrtf(gradv[1][0] * gradv[1][0] + gradv[1][1] * gradv[1][1] +
+ gradv[1][2] * gradv[1][2]);
+ if (gradtrue) {
+ gradtrue *= p->primitives.limiter.maxr;
+ const float gradmax = p->primitives.limiter.v[1][1] - p->primitives.v[1];
+ const float gradmin = p->primitives.v[1] - p->primitives.limiter.v[1][0];
+ const float gradtrue_inv = 1.f / gradtrue;
+ const float alpha =
+ min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
+ p->primitives.gradients.v[1][0] *= alpha;
+ p->primitives.gradients.v[1][1] *= alpha;
+ p->primitives.gradients.v[1][2] *= alpha;
+ }
+
+ gradtrue = sqrtf(gradv[2][0] * gradv[2][0] + gradv[2][1] * gradv[2][1] +
+ gradv[2][2] * gradv[2][2]);
+ if (gradtrue) {
+ gradtrue *= p->primitives.limiter.maxr;
+ const float gradmax = p->primitives.limiter.v[2][1] - p->primitives.v[2];
+ const float gradmin = p->primitives.v[2] - p->primitives.limiter.v[2][0];
+ const float gradtrue_inv = 1.f / gradtrue;
+ const float alpha =
+ min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
+ p->primitives.gradients.v[2][0] *= alpha;
+ p->primitives.gradients.v[2][1] *= alpha;
+ p->primitives.gradients.v[2][2] *= alpha;
+ }
+
+ gradtrue =
+ sqrtf(gradP[0] * gradP[0] + gradP[1] * gradP[1] + gradP[2] * gradP[2]);
+ if (gradtrue) {
+ gradtrue *= p->primitives.limiter.maxr;
+ const float gradmax = p->primitives.limiter.P[1] - p->primitives.P;
+ const float gradmin = p->primitives.P - p->primitives.limiter.P[0];
+ const float gradtrue_inv = 1.f / gradtrue;
+ const float alpha =
+ min3(1.0f, gradmax * gradtrue_inv, gradmin * gradtrue_inv);
+ p->primitives.gradients.P[0] *= alpha;
+ p->primitives.gradients.P[1] *= alpha;
+ p->primitives.gradients.P[2] *= alpha;
+ }
+}
+
+#endif /* SWIFT_GIZMO_MFV_SLOPE_LIMITER_CELL_H */
diff --git a/src/hydro/GizmoMFV/hydro_slope_limiters_face.h b/src/hydro/GizmoMFV/hydro_slope_limiters_face.h
new file mode 100644
index 0000000000000000000000000000000000000000..b12f15590f51c16f494978d9501401709b2cbf59
--- /dev/null
+++ b/src/hydro/GizmoMFV/hydro_slope_limiters_face.h
@@ -0,0 +1,128 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/**
+ * @brief Slope limit a single quantity at the interface
+ *
+ * @param phi_i Value of the quantity at the particle position.
+ * @param phi_j Value of the quantity at the neighbouring particle position.
+ * @param phi_mid0 Extrapolated value of the quantity at the interface position.
+ * @param xij_norm Distance between the particle position and the interface
+ * position.
+ * @param r Distance between the particle and its neighbour.
+ * @return The slope limited difference between the quantity at the particle
+ * position and the quantity at the interface position.
+ */
+#ifndef SWIFT_GIZMO_MFV_SLOPE_LIMITER_FACE_H
+#define SWIFT_GIZMO_MFV_SLOPE_LIMITER_FACE_H
+
+/* Some standard headers. */
+#include <float.h>
+
+/* Local headers. */
+#include "minmax.h"
+#include "sign.h"
+
+__attribute__((always_inline)) INLINE static float
+hydro_slope_limit_face_quantity(float phi_i, float phi_j, float phi_mid0,
+ float xij_norm, float r_inv) {
+
+ const float psi1 = 0.5f;
+ const float psi2 = 0.25f;
+
+ const float delta1 = psi1 * fabsf(phi_i - phi_j);
+ const float delta2 = psi2 * fabsf(phi_i - phi_j);
+
+ const float phimin = min(phi_i, phi_j);
+ const float phimax = max(phi_i, phi_j);
+
+ const float phibar = phi_i + xij_norm * r_inv * (phi_j - phi_i);
+
+ float phiplus, phiminus, phi_mid;
+
+ if (same_signf(phimax + delta1, phimax))
+ phiplus = phimax + delta1;
+ else
+ phiplus = phimax / (1.0f + delta1 / (fabsf(phimax) + FLT_MIN));
+
+ if (same_signf(phimin - delta1, phimin))
+ phiminus = phimin - delta1;
+ else
+ phiminus = phimin / (1.0f + delta1 / (fabsf(phimin) + FLT_MIN));
+
+ if (phi_i < phi_j) {
+ const float temp = min(phibar + delta2, phi_mid0);
+ phi_mid = max(phiminus, temp);
+ } else {
+ const float temp = max(phibar - delta2, phi_mid0);
+ phi_mid = min(phiplus, temp);
+ }
+
+ return phi_mid - phi_i;
+}
+
+/**
+ * @brief Slope limit the slopes at the interface between two particles
+ *
+ * @param Wi Hydrodynamic variables of particle i.
+ * @param Wj Hydrodynamic variables of particle j.
+ * @param dWi Difference between the hydrodynamic variables of particle i at the
+ * position of particle i and at the interface position.
+ * @param dWj Difference between the hydrodynamic variables of particle j at the
+ * position of particle j and at the interface position.
+ * @param xij_i Relative position vector of the interface w.r.t. particle i.
+ * @param xij_j Relative position vector of the interface w.r.t. partilce j.
+ * @param r Distance between particle i and particle j.
+ */
+__attribute__((always_inline)) INLINE static void hydro_slope_limit_face(
+ float *Wi, float *Wj, float *dWi, float *dWj, const float *xij_i,
+ const float *xij_j, float r) {
+
+ const float xij_i_norm =
+ sqrtf(xij_i[0] * xij_i[0] + xij_i[1] * xij_i[1] + xij_i[2] * xij_i[2]);
+
+ const float xij_j_norm =
+ sqrtf(xij_j[0] * xij_j[0] + xij_j[1] * xij_j[1] + xij_j[2] * xij_j[2]);
+
+ const float r_inv = 1.f / r;
+
+ dWi[0] = hydro_slope_limit_face_quantity(Wi[0], Wj[0], Wi[0] + dWi[0],
+ xij_i_norm, r_inv);
+ dWi[1] = hydro_slope_limit_face_quantity(Wi[1], Wj[1], Wi[1] + dWi[1],
+ xij_i_norm, r_inv);
+ dWi[2] = hydro_slope_limit_face_quantity(Wi[2], Wj[2], Wi[2] + dWi[2],
+ xij_i_norm, r_inv);
+ dWi[3] = hydro_slope_limit_face_quantity(Wi[3], Wj[3], Wi[3] + dWi[3],
+ xij_i_norm, r_inv);
+ dWi[4] = hydro_slope_limit_face_quantity(Wi[4], Wj[4], Wi[4] + dWi[4],
+ xij_i_norm, r_inv);
+
+ dWj[0] = hydro_slope_limit_face_quantity(Wj[0], Wi[0], Wj[0] + dWj[0],
+ xij_j_norm, r_inv);
+ dWj[1] = hydro_slope_limit_face_quantity(Wj[1], Wi[1], Wj[1] + dWj[1],
+ xij_j_norm, r_inv);
+ dWj[2] = hydro_slope_limit_face_quantity(Wj[2], Wi[2], Wj[2] + dWj[2],
+ xij_j_norm, r_inv);
+ dWj[3] = hydro_slope_limit_face_quantity(Wj[3], Wi[3], Wj[3] + dWj[3],
+ xij_j_norm, r_inv);
+ dWj[4] = hydro_slope_limit_face_quantity(Wj[4], Wi[4], Wj[4] + dWj[4],
+ xij_j_norm, r_inv);
+}
+
+#endif /* SWIFT_GIZMO_MFV_SLOPE_LIMITER_FACE_H */
diff --git a/src/hydro/GizmoMFV/hydro_unphysical.h b/src/hydro/GizmoMFV/hydro_unphysical.h
new file mode 100644
index 0000000000000000000000000000000000000000..81c644e6cce00e0d871aafc39140e420e042a926
--- /dev/null
+++ b/src/hydro/GizmoMFV/hydro_unphysical.h
@@ -0,0 +1,69 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_HYDRO_UNPHYSICAL_H
+#define SWIFT_HYDRO_UNPHYSICAL_H
+
+#if defined(GIZMO_UNPHYSICAL_ERROR) || defined(GIZMO_UNPHYSICAL_RESCUE)
+
+#if defined(GIZMO_UNPHYSICAL_ERROR)
+
+/*! @brief Crash whenever an unphysical value is detected. */
+#define gizmo_unphysical_message(name, quantity) \
+ error("Unphysical " name " detected (%g)!", quantity);
+
+#elif defined(GIZMO_UNPHYSICAL_WARNING)
+
+/*! @brief Show a warning whenever an unphysical value is detected. */
+#define gizmo_unphysical_message(name, quantity) \
+ message("Unphysical " name " detected (%g), reset to 0!", quantity);
+
+#else
+
+/*! @brief Don't tell anyone an unphysical value was detected. */
+#define gizmo_unphysical_message(name, quantity)
+
+#endif
+
+#define gizmo_check_physical_quantity(name, quantity) \
+ if (quantity < 0.f) { \
+ gizmo_unphysical_message(name, quantity); \
+ quantity = 0.f; \
+ }
+
+#define gizmo_check_physical_quantities( \
+ mass_name, energy_name, mass, momentum_x, momentum_y, momentum_z, energy) \
+ gizmo_check_physical_quantity(mass_name, mass); \
+ gizmo_check_physical_quantity(energy_name, energy); \
+ /* now check for vacuum and make sure we have a real vacuum */ \
+ if (mass == 0.f || energy == 0.f) { \
+ mass = 0.f; \
+ momentum_x = 0.f; \
+ momentum_y = 0.f; \
+ momentum_z = 0.f; \
+ energy = 0.f; \
+ }
+
+#else // defined(GIZMO_UNPHYSICAL_ERROR) || defined(GIZMO_UNPHYSICAL_RESCUE)
+
+#define gizmo_check_physical_quantities( \
+ mass_name, energy_name, mass, momentum_x, momentum_y, momentum_z, energy)
+
+#endif // defined(GIZMO_UNPHYSICAL_ERROR) || defined(GIZMO_UNPHYSICAL_RESCUE)
+
+#endif /* SWIFT_HYDRO_UNPHYSICAL_H */
diff --git a/src/hydro/GizmoMFV/hydro_velocities.h b/src/hydro/GizmoMFV/hydro_velocities.h
new file mode 100644
index 0000000000000000000000000000000000000000..ea30d5a6c9c74d34ffd73fa6ab941640f37b02e4
--- /dev/null
+++ b/src/hydro/GizmoMFV/hydro_velocities.h
@@ -0,0 +1,157 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2017 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_HYDRO_VELOCITIES_H
+#define SWIFT_HYDRO_VELOCITIES_H
+
+/**
+ * @brief Initialize the GIZMO particle velocities before the start of the
+ * actual run based on the initial value of the primitive velocity.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_velocities_init(
+ struct part* restrict p, struct xpart* restrict xp) {
+
+#ifdef GIZMO_FIX_PARTICLES
+ p->v[0] = 0.f;
+ p->v[1] = 0.f;
+ p->v[2] = 0.f;
+#else
+ p->v[0] = p->primitives.v[0];
+ p->v[1] = p->primitives.v[1];
+ p->v[2] = p->primitives.v[2];
+#endif
+
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
+}
+
+/**
+ * @brief Set the particle velocity field that will be used to deboost fluid
+ * velocities during the force loop.
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_velocities_prepare_force(struct part* restrict p,
+ const struct xpart* restrict xp) {}
+
+/**
+ * @brief Set the variables that will be used to update the smoothing length
+ * during the drift (these will depend on the movement of the particles).
+ *
+ * @param p The particle to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_velocities_end_force(
+ struct part* restrict p) {
+
+#ifdef GIZMO_FIX_PARTICLES
+ /* disable the smoothing length update, since the smoothing lengths should
+ stay the same for all steps (particles don't move) */
+ p->force.h_dt = 0.0f;
+#else
+ /* Add normalization to h_dt. */
+ p->force.h_dt *= p->h * hydro_dimension_inv;
+#endif
+}
+
+/**
+ * @brief Set the velocity of a GIZMO particle, based on the values of its
+ * primitive variables and the geometry of its mesh-free "cell".
+ *
+ * @param p The particle to act upon.
+ * @param xp The extended particle data to act upon.
+ */
+__attribute__((always_inline)) INLINE static void hydro_velocities_set(
+ struct part* restrict p, struct xpart* restrict xp) {
+
+/* We first set the particle velocity. */
+#ifdef GIZMO_FIX_PARTICLES
+
+ p->v[0] = 0.f;
+ p->v[1] = 0.f;
+ p->v[2] = 0.f;
+
+#else // GIZMO_FIX_PARTICLES
+
+ if (p->conserved.mass > 0.f && p->primitives.rho > 0.f) {
+
+ const float inverse_mass = 1.f / p->conserved.mass;
+
+ /* Normal case: set particle velocity to fluid velocity. */
+ p->v[0] = p->conserved.momentum[0] * inverse_mass;
+ p->v[1] = p->conserved.momentum[1] * inverse_mass;
+ p->v[2] = p->conserved.momentum[2] * inverse_mass;
+
+#ifdef GIZMO_STEER_MOTION
+
+ /* Add a correction to the velocity to keep particle positions close enough
+ to
+ the centroid of their mesh-free "cell". */
+ /* The correction term below is the same one described in Springel (2010).
+ */
+ float ds[3];
+ ds[0] = p->geometry.centroid[0];
+ ds[1] = p->geometry.centroid[1];
+ ds[2] = p->geometry.centroid[2];
+ const float d = sqrtf(ds[0] * ds[0] + ds[1] * ds[1] + ds[2] * ds[2]);
+ const float R = get_radius_dimension_sphere(p->geometry.volume);
+ const float eta = 0.25f;
+ const float etaR = eta * R;
+ const float xi = 1.f;
+ const float soundspeed =
+ sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+ /* We only apply the correction if the offset between centroid and position
+ is too large. */
+ if (d > 0.9f * etaR) {
+ float fac = xi * soundspeed / d;
+ if (d < 1.1f * etaR) {
+ fac *= 5.f * (d - 0.9f * etaR) / etaR;
+ }
+ p->v[0] -= ds[0] * fac;
+ p->v[1] -= ds[1] * fac;
+ p->v[2] -= ds[2] * fac;
+ }
+
+#endif // GIZMO_STEER_MOTION
+ } else {
+ /* Vacuum particles have no fluid velocity. */
+ p->v[0] = 0.f;
+ p->v[1] = 0.f;
+ p->v[2] = 0.f;
+ }
+
+#endif // GIZMO_FIX_PARTICLES
+
+ /* Now make sure all velocity variables are up to date. */
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
+
+ if (p->gpart) {
+ p->gpart->v_full[0] = p->v[0];
+ p->gpart->v_full[1] = p->v[1];
+ p->gpart->v_full[2] = p->v[2];
+ }
+}
+
+#endif /* SWIFT_HYDRO_VELOCITIES_H */
diff --git a/src/hydro/MinimalMultiMat/hydro_io.h b/src/hydro/MinimalMultiMat/hydro_io.h
index 6e55497dc42111de291b8984a7dca048bb774722..2a5eeb6a54d079ae72e1591116a8984b0d7a6f38 100644
--- a/src/hydro/MinimalMultiMat/hydro_io.h
+++ b/src/hydro/MinimalMultiMat/hydro_io.h
@@ -72,6 +72,12 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
io_make_input_field("MaterialID", INT, 1, OPTIONAL, 1, parts, mat_id);
}
+void convert_S(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ ret[0] = hydro_get_comoving_entropy(p);
+}
+
void convert_P(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
@@ -146,7 +152,7 @@ void convert_part_potential(const struct engine* e, const struct part* p,
void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
struct io_props* list, int* num_fields) {
- *num_fields = 10;
+ *num_fields = 11;
/* List what we want to write */
list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
@@ -164,13 +170,16 @@ void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
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("MaterialID", INT, 1, UNIT_CONV_NO_UNITS,
+ list[7] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
+ UNIT_CONV_ENTROPY_PER_UNIT_MASS,
+ parts, xparts, convert_S);
+ list[8] = io_make_output_field("MaterialID", INT, 1, UNIT_CONV_NO_UNITS,
parts, mat_id);
- list[8] = io_make_output_field_convert_part(
+ list[9] = io_make_output_field_convert_part(
"Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, xparts, convert_P);
- list[9] = io_make_output_field_convert_part("Potential", FLOAT, 1,
- UNIT_CONV_POTENTIAL, parts,
- xparts, convert_part_potential);
+ list[10] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+ UNIT_CONV_POTENTIAL, parts,
+ xparts, convert_part_potential);
}
/**
diff --git a/src/hydro_io.h b/src/hydro_io.h
index b1d43592d504f82a2db1a9f265d30c20d6eddc35..d752bb8bc03f619fe759fc8f5de32a01b3a61abe 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -33,8 +33,10 @@
#include "./hydro/PressureEnergy/hydro_io.h"
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_io.h"
-#elif defined(GIZMO_SPH)
-#include "./hydro/Gizmo/hydro_io.h"
+#elif defined(GIZMO_MFV_SPH)
+#include "./hydro/GizmoMFV/hydro_io.h"
+#elif defined(GIZMO_MFM_SPH)
+#include "./hydro/GizmoMFM/hydro_io.h"
#elif defined(SHADOWFAX_SPH)
#include "./hydro/Shadowswift/hydro_io.h"
#elif defined(MINIMAL_MULTI_MAT_SPH)
diff --git a/src/part.h b/src/part.h
index a207e82c2ee19f6ddad6dd261167a3360e23bfdd..7a1b6dd2b0e6a8f062658c5fac41f3534baa522e 100644
--- a/src/part.h
+++ b/src/part.h
@@ -57,8 +57,12 @@
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_part.h"
#define hydro_need_extra_init_loop 0
-#elif defined(GIZMO_SPH)
-#include "./hydro/Gizmo/hydro_part.h"
+#elif defined(GIZMO_MFV_SPH)
+#include "./hydro/GizmoMFV/hydro_part.h"
+#define hydro_need_extra_init_loop 0
+#define EXTRA_HYDRO_LOOP
+#elif defined(GIZMO_MFM_SPH)
+#include "./hydro/GizmoMFM/hydro_part.h"
#define hydro_need_extra_init_loop 0
#define EXTRA_HYDRO_LOOP
#elif defined(SHADOWFAX_SPH)
diff --git a/src/riemann/riemann_exact.h b/src/riemann/riemann_exact.h
index 0c6c1c86273b42f521102e6e91d642d3dff30b27..9a671559c12a9df6d1e11016adc619c6043c83a9 100644
--- a/src/riemann/riemann_exact.h
+++ b/src/riemann/riemann_exact.h
@@ -485,6 +485,87 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
Whalf[3] += vhalf * n_unit[2];
}
+/**
+ * @brief Solve the Riemann problem between the given left and right state and
+ * return the velocity and pressure in the middle state
+ *
+ * Based on chapter 4 in Toro
+ *
+ * @param WL Left state.
+ * @param vL Left state velocity.
+ * @param WR Right state.
+ * @param vR Right state velocity.
+ * @param vM Middle state velocity.
+ * @param PM Middle state pressure.
+ */
+__attribute__((always_inline)) INLINE static void
+riemann_solver_solve_middle_state(const float* WL, const float vL,
+ const float* WR, const float vR, float* vM,
+ float* PM) {
+
+ /* sound speeds */
+ float aL, aR;
+ /* variables used for finding pstar */
+ float p, pguess, fp, fpguess;
+
+ /* calculate sound speeds */
+ aL = sqrtf(hydro_gamma * WL[4] / WL[0]);
+ aR = sqrtf(hydro_gamma * WR[4] / WR[0]);
+
+ /* vacuum */
+ /* check vacuum (generation) condition */
+ if (riemann_is_vacuum(WL, WR, vL, vR, aL, aR)) {
+ *vM = 0.0f;
+ *PM = 0.0f;
+ return;
+ }
+
+ /* values are ok: let's find pstar (riemann_f(pstar) = 0)! */
+ /* We normally use a Newton-Raphson iteration to find the zeropoint
+ of riemann_f(p), but if pstar is close to 0, we risk negative p values.
+ Since riemann_f(p) is undefined for negative pressures, we don't
+ want this to happen.
+ We therefore use Brent's method if riemann_f(0) is larger than some
+ value. -5 makes the iteration fail safe while almost never invoking
+ the expensive Brent solver. */
+ p = 0.0f;
+ /* obtain a first guess for p */
+ pguess = riemann_guess_p(WL, WR, vL, vR, aL, aR);
+ fp = riemann_f(p, WL, WR, vL, vR, aL, aR);
+ fpguess = riemann_f(pguess, WL, WR, vL, vR, aL, aR);
+ /* ok, pstar is close to 0, better use Brent's method... */
+ /* we use Newton-Raphson until we find a suitable interval */
+ if (fp * fpguess >= 0.0f) {
+ /* Newton-Raphson until convergence or until suitable interval is found
+ to use Brent's method */
+ unsigned int counter = 0;
+ while (fabs(p - pguess) > 1.e-6f * 0.5f * (p + pguess) && fpguess < 0.0f) {
+ p = pguess;
+ pguess = pguess - fpguess / riemann_fprime(pguess, WL, WR, aL, aR);
+ fpguess = riemann_f(pguess, WL, WR, vL, vR, aL, aR);
+ counter++;
+ if (counter > 1000) {
+ error("Stuck in Newton-Raphson!\n");
+ }
+ }
+ }
+ /* As soon as there is a suitable interval: use Brent's method */
+ if (1.e6 * fabs(p - pguess) > 0.5f * (p + pguess) && fpguess > 0.0f) {
+ p = 0.0f;
+ fp = riemann_f(p, WL, WR, vL, vR, aL, aR);
+ /* use Brent's method to find the zeropoint */
+ p = riemann_solve_brent(p, pguess, fp, fpguess, 1.e-6, WL, WR, vL, vR, aL,
+ aR);
+ } else {
+ p = pguess;
+ }
+
+ *PM = p;
+ /* calculate the velocity in the intermediate state */
+ *vM =
+ 0.5f * (vL + vR) + 0.5f * (riemann_fb(p, WR, aR) - riemann_fb(p, WL, aL));
+}
+
__attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
const float* Wi, const float* Wj, const float* n_unit, const float* vij,
float* totflux) {
@@ -543,4 +624,43 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
#endif
}
+__attribute__((always_inline)) INLINE static void
+riemann_solve_for_middle_state_flux(const float* Wi, const float* Wj,
+ const float* n_unit, const float* vij,
+ float* totflux) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ riemann_check_input(Wi, Wj, n_unit, vij);
+#endif
+
+ /* vacuum? */
+ if (Wi[0] == 0.0f || Wj[0] == 0.0f) {
+ totflux[0] = 0.0f;
+ totflux[1] = 0.0f;
+ totflux[2] = 0.0f;
+ totflux[3] = 0.0f;
+ totflux[4] = 0.0f;
+ return;
+ }
+
+ const float vL = Wi[1] * n_unit[0] + Wi[2] * n_unit[1] + Wi[3] * n_unit[2];
+ const float vR = Wj[1] * n_unit[0] + Wj[2] * n_unit[1] + Wj[3] * n_unit[2];
+
+ float vM, PM;
+ riemann_solver_solve_middle_state(Wi, vL, Wj, vR, &vM, &PM);
+
+ const float vface =
+ vij[0] * n_unit[0] + vij[1] * n_unit[1] + vij[2] * n_unit[2];
+
+ totflux[0] = 0.0f;
+ totflux[1] = PM * n_unit[0];
+ totflux[2] = PM * n_unit[1];
+ totflux[3] = PM * n_unit[2];
+ totflux[4] = (vM + vface) * PM;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ riemann_check_output(Wi, Wj, n_unit, vij, totflux);
+#endif
+}
+
#endif /* SWIFT_RIEMANN_EXACT_H */
diff --git a/src/riemann/riemann_hllc.h b/src/riemann/riemann_hllc.h
index 26fcd627fcf11769baf18955af53ed6e948ba513..68a92111d8e075faff5f93ab9c3d8a1ebc496cb0 100644
--- a/src/riemann/riemann_hllc.h
+++ b/src/riemann/riemann_hllc.h
@@ -185,4 +185,77 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
#endif
}
+__attribute__((always_inline)) INLINE static void
+riemann_solve_for_middle_state_flux(const float *WL, const float *WR,
+ const float *n, const float *vij,
+ float *totflux) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ riemann_check_input(WL, WR, n, vij);
+#endif
+
+ /* Handle pure vacuum */
+ if (!WL[0] && !WR[0]) {
+ totflux[0] = 0.f;
+ totflux[1] = 0.f;
+ totflux[2] = 0.f;
+ totflux[3] = 0.f;
+ totflux[4] = 0.f;
+ return;
+ }
+
+ /* STEP 0: obtain velocity in interface frame */
+ const float uL = WL[1] * n[0] + WL[2] * n[1] + WL[3] * n[2];
+ const float uR = WR[1] * n[0] + WR[2] * n[1] + WR[3] * n[2];
+ const float aL = sqrtf(hydro_gamma * WL[4] / WL[0]);
+ const float aR = sqrtf(hydro_gamma * WR[4] / WR[0]);
+
+ /* Handle vacuum: vacuum does not require iteration and is always exact */
+ if (riemann_is_vacuum(WL, WR, uL, uR, aL, aR)) {
+ totflux[0] = 0.f;
+ totflux[1] = 0.f;
+ totflux[2] = 0.f;
+ totflux[3] = 0.f;
+ totflux[4] = 0.f;
+ return;
+ }
+
+ /* STEP 1: pressure estimate */
+ const float rhobar = 0.5f * (WL[0] + WR[0]);
+ const float abar = 0.5f * (aL + aR);
+ const float pPVRS = 0.5f * (WL[4] + WR[4]) - 0.5f * (uR - uL) * rhobar * abar;
+ const float pstar = max(0.f, pPVRS);
+
+ /* STEP 2: wave speed estimates
+ all these speeds are along the interface normal, since uL and uR are */
+ float qL = 1.f;
+ if (pstar > WL[4] && WL[4] > 0.f) {
+ qL = sqrtf(1.f +
+ 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
+ (pstar / WL[4] - 1.f));
+ }
+ float qR = 1.f;
+ if (pstar > WR[4] && WR[4] > 0.f) {
+ qR = sqrtf(1.f +
+ 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
+ (pstar / WR[4] - 1.f));
+ }
+ const float SL = uL - aL * qL;
+ const float SR = uR + aR * qR;
+ const float Sstar =
+ (WR[4] - WL[4] + WL[0] * uL * (SL - uL) - WR[0] * uR * (SR - uR)) /
+ (WL[0] * (SL - uL) - WR[0] * (SR - uR));
+
+ totflux[0] = 0.0f;
+ totflux[1] = pstar * n[0];
+ totflux[2] = pstar * n[1];
+ totflux[3] = pstar * n[2];
+ const float vface = vij[0] * n[0] + vij[1] * n[1] + vij[2] * n[2];
+ totflux[4] = pstar * (Sstar + vface);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ riemann_check_output(WL, WR, n, vij, totflux);
+#endif
+}
+
#endif /* SWIFT_RIEMANN_HLLC_H */
diff --git a/src/riemann/riemann_trrs.h b/src/riemann/riemann_trrs.h
index 597f6f8edd9690dbdd34a3c9885ae228c2af734c..c19ff70f62a9018f4d121f8b689a0a13c371d86a 100644
--- a/src/riemann/riemann_trrs.h
+++ b/src/riemann/riemann_trrs.h
@@ -219,4 +219,64 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
#endif
}
+__attribute__((always_inline)) INLINE static void
+riemann_solve_for_middle_state_flux(const float* Wi, const float* Wj,
+ const float* n_unit, const float* vij,
+ float* totflux) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ riemann_check_input(Wi, Wj, n_unit, vij);
+#endif
+
+ if (Wi[0] == 0.0f || Wj[0] == 0.0f) {
+ totflux[0] = 0.0f;
+ totflux[1] = 0.0f;
+ totflux[2] = 0.0f;
+ totflux[3] = 0.0f;
+ totflux[4] = 0.0f;
+ return;
+ }
+
+ /* calculate the velocities along the interface normal */
+ const float vL = Wi[1] * n_unit[0] + Wi[2] * n_unit[1] + Wi[3] * n_unit[2];
+ const float vR = Wj[1] * n_unit[0] + Wj[2] * n_unit[1] + Wj[3] * n_unit[2];
+
+ /* calculate the sound speeds */
+ const float aL = sqrtf(hydro_gamma * Wi[4] / Wi[0]);
+ const float aR = sqrtf(hydro_gamma * Wj[4] / Wj[0]);
+
+ if (riemann_is_vacuum(Wi, Wj, vL, vR, aL, aR)) {
+ totflux[0] = 0.0f;
+ totflux[1] = 0.0f;
+ totflux[2] = 0.0f;
+ totflux[3] = 0.0f;
+ totflux[4] = 0.0f;
+ return;
+ }
+
+ /* calculate the velocity and pressure in the intermediate state */
+ const float PLR = pow_gamma_minus_one_over_two_gamma(Wi[4] / Wj[4]);
+ const float ustar = (PLR * vL / aL + vR / aR +
+ hydro_two_over_gamma_minus_one * (PLR - 1.0f)) /
+ (PLR / aL + 1.0f / aR);
+ const float pstar =
+ 0.5f *
+ (Wi[4] * pow_two_gamma_over_gamma_minus_one(
+ 1.0f + hydro_gamma_minus_one_over_two / aL * (vL - ustar)) +
+ Wj[4] * pow_two_gamma_over_gamma_minus_one(
+ 1.0f + hydro_gamma_minus_one_over_two / aR * (ustar - vR)));
+
+ totflux[0] = 0.0f;
+ totflux[1] = pstar * n_unit[0];
+ totflux[2] = pstar * n_unit[1];
+ totflux[3] = pstar * n_unit[2];
+ const float vface =
+ vij[0] * n_unit[0] + vij[1] * n_unit[1] + vij[2] * n_unit[2];
+ totflux[4] = pstar * (ustar + vface);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ riemann_check_output(Wi, Wj, n_unit, vij, totflux);
+#endif
+}
+
#endif /* SWIFT_RIEMANN_TRRS_H */
diff --git a/src/utilities.h b/src/utilities.h
new file mode 100644
index 0000000000000000000000000000000000000000..28faf6b1e59159af4ee4a7b87054ab6a36478c9c
--- /dev/null
+++ b/src/utilities.h
@@ -0,0 +1,59 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Jacob Kegerreis (jacob.kegerreis@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_UTILITIES_H
+#define SWIFT_UTILITIES_H
+
+/**
+ * @brief Search for a value in a monotonically increasing array to find the
+ * index such that array[index] < value < array[index + 1]
+ *
+ * @param x The value to find
+ * @param array The array to search
+ * @param n The length of the array
+ *
+ * Return -1 and n for x below and above the array edge values respectively.
+ */
+INLINE static int find_value_in_monot_incr_array(const float x,
+ const float *array,
+ const int n) {
+
+ int index_mid, index_low = 0, index_high = n;
+
+ // Until array[index_low] < x < array[index_high=index_low+1]
+ while (index_high - index_low > 1) {
+ index_mid = (index_high + index_low) / 2; // Middle index
+
+ // Replace the low or high index with the middle
+ if (array[index_mid] <= x)
+ index_low = index_mid;
+ else
+ index_high = index_mid;
+ }
+
+ // Set index with the found index_low or an error value if outside the array
+ if (x < array[0])
+ return -1;
+ else if (array[n - 1] <= x)
+ return n;
+ else
+ return index_low;
+}
+
+#endif /* SWIFT_UTILITIES_H */
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 2dfe760d80332f7da8de4cf11c0f15e21c231089..891eef3f518f83c17b66623e3dac1832512d31f3 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -27,7 +27,7 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
testMatrixInversion testThreadpool testDump testLogger testInteractions.sh \
testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf \
- testPotentialPair testEOS
+ testPotentialPair testEOS testUtilities
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
@@ -37,7 +37,7 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testAdiabaticIndex testRiemannExact testRiemannTRRS \
testRiemannHLLC testMatrixInversion testDump testLogger \
testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
- testGravityDerivatives testPotentialSelf testPotentialPair testEOS
+ testGravityDerivatives testPotentialSelf testPotentialPair testEOS testUtilities
# Rebuild tests when SWIFT is updated.
$(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -110,6 +110,8 @@ testPotentialPair_SOURCES = testPotentialPair.c
testEOS_SOURCES = testEOS.c
+testUtilities_SOURCES = testUtilities.c
+
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testActivePair.sh \
test27cells.sh test27cellsPerturbed.sh testParser.sh testPeriodicBC.sh \
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 5195f431c23067a732ae9d060ffc6d2b01b5feaf..a1661bb93bfdac3b6111148f9fae42b243c04742 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -122,7 +122,7 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(MINIMAL_MULTI_MAT_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
-#elif defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
+#elif defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
part->primitives.P = pressure;
#else
error("Need to define pressure here !");
@@ -222,9 +222,9 @@ void reset_particles(struct cell *c, struct hydro_space *hs,
hydro_init_part(p, hs);
-#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
+#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
float volume = p->conserved.mass / density;
-#if defined(GIZMO_SPH)
+#if defined(GIZMO_MFV_SPH)
p->geometry.volume = volume;
#else
p->cell.volume = volume;
@@ -300,7 +300,7 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
part->h = size * h / (float)n;
h_max = fmax(h_max, part->h);
-#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
+#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
part->conserved.mass = density * volume / count;
#else
part->mass = density * volume / count;
@@ -314,7 +314,7 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
part->id = ++(*partId);
part->time_bin = 1;
-#if defined(GIZMO_SPH)
+#if defined(GIZMO_MFV_SPH)
part->geometry.volume = part->conserved.mass / density;
part->primitives.rho = density;
part->primitives.v[0] = part->v[0];
@@ -406,7 +406,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->parts[pid].v[2], main_cell->parts[pid].h,
hydro_get_comoving_density(&main_cell->parts[pid]),
#if defined(MINIMAL_SPH) || defined(MINIMAL_MULTI_MAT_SPH) || \
- defined(SHADOWFAX_SPH) || defined(HOPKINS_PU_SPH)
+ defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH) || defined(HOPKINS_PU_SPH)
0.f,
#else
main_cell->parts[pid].density.div_v,
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 1934ca84b2610ec565503a319bd798f88c279cce..bf226335793060b338fadf97756d22e4abdc518f 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -150,7 +150,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
h_max = fmaxf(h_max, part->h);
part->id = ++(*partId);
-#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
+#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
part->conserved.mass = density * volume / count;
#ifdef SHADOWFAX_SPH
@@ -261,7 +261,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
main_cell->parts[pid].v[2],
hydro_get_comoving_density(&main_cell->parts[pid]),
-#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
+#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
0.f,
#elif defined(HOPKINS_PU_SPH)
main_cell->parts[pid].density.pressure_bar_dh,
@@ -300,7 +300,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
cj->parts[pjd].id, cj->parts[pjd].x[0], cj->parts[pjd].x[1],
cj->parts[pjd].x[2], cj->parts[pjd].v[0], cj->parts[pjd].v[1],
cj->parts[pjd].v[2], hydro_get_comoving_density(&cj->parts[pjd]),
-#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
+#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
0.f,
#else
main_cell->parts[pjd].density.rho_dh,
diff --git a/tests/testActivePair.c b/tests/testActivePair.c
index 4b0b63d9766dbfd313fd60610705ede3eea29141..0453f6d5896eaa53b0f44a567d353d7d8e8fb7df 100644
--- a/tests/testActivePair.c
+++ b/tests/testActivePair.c
@@ -94,7 +94,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
part->id = ++(*partId);
/* Set the mass */
-#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
+#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
part->conserved.mass = density * volume / count;
#ifdef SHADOWFAX_SPH
@@ -105,7 +105,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
#else
part->mass = density * volume / count;
-#endif /* GIZMO_SPH */
+#endif /* GIZMO_MFV_SPH */
/* Set the thermodynamic variable */
#if defined(GADGET2_SPH)
diff --git a/tests/testInteractions.c b/tests/testInteractions.c
index 8e29bf9b3ceb519d8a6a8271f290e0c9121d3362..f288ab3c9ac6c625075b0334e1a0e26a3d5a082b 100644
--- a/tests/testInteractions.c
+++ b/tests/testInteractions.c
@@ -74,7 +74,7 @@ struct part *make_particles(size_t count, double *offset, double spacing,
p->h = h;
p->id = ++(*partId);
-#if !defined(GIZMO_SPH) && !defined(SHADOWFAX_SPH)
+#if !defined(GIZMO_MFV_SPH) && !defined(SHADOWFAX_SPH)
p->mass = 1.0f;
#endif
@@ -107,8 +107,9 @@ struct part *make_particles(size_t count, double *offset, double spacing,
*/
void prepare_force(struct part *parts, size_t count) {
-#if !defined(GIZMO_SPH) && !defined(SHADOWFAX_SPH) && !defined(MINIMAL_SPH) && \
- !defined(MINIMAL_MULTI_MAT_SPH) && !defined(HOPKINS_PU_SPH)
+#if !defined(GIZMO_MFV_SPH) && !defined(SHADOWFAX_SPH) && \
+ !defined(MINIMAL_SPH) && !defined(MINIMAL_MULTI_MAT_SPH) && \
+ !defined(HOPKINS_PU_SPH)
struct part *p;
for (size_t i = 0; i < count; ++i) {
p = &parts[i];
diff --git a/tests/testPeriodicBC.c b/tests/testPeriodicBC.c
index 2aa15f0aa94acafa20c56767fa6d739798f6c7f3..385de9752f361f4f015eb64a466473324901030f 100644
--- a/tests/testPeriodicBC.c
+++ b/tests/testPeriodicBC.c
@@ -129,7 +129,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
h_max = fmax(h_max, part->h);
part->id = ++(*partId);
-#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
+#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
part->conserved.mass = density * volume / count;
#ifdef SHADOWFAX_SPH
@@ -237,7 +237,7 @@ void dump_particle_fields(char *fileName, struct cell *main_cell, int i, int j,
main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
main_cell->parts[pid].v[2],
hydro_get_comoving_density(&main_cell->parts[pid]),
-#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
+#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
0.f,
#else
main_cell->parts[pid].density.rho_dh,
diff --git a/tests/testRiemannExact.c b/tests/testRiemannExact.c
index b72b7b6ce9d7fb4c6f7449af73c6a67ce4d8d1e7..bce7c52d422f966e10d530cdcbc8f6d20431e153 100644
--- a/tests/testRiemannExact.c
+++ b/tests/testRiemannExact.c
@@ -23,6 +23,12 @@
#include <stdio.h>
#include <string.h>
+/* Force use of exact Riemann solver */
+#undef RIEMANN_SOLVER_TRRS
+#undef RIEMANN_SOLVER_HLLC
+#undef RIEMANN_SOLVER_EXACT
+#define RIEMANN_SOLVER_EXACT 1
+
/* Local headers. */
#include "riemann/riemann_exact.h"
#include "swift.h"
diff --git a/tests/testRiemannHLLC.c b/tests/testRiemannHLLC.c
index e843b29093f88d54dd931becf966d46670707091..b988825eb0535fcdc46baa1db1203d0dbac3537a 100644
--- a/tests/testRiemannHLLC.c
+++ b/tests/testRiemannHLLC.c
@@ -24,6 +24,12 @@
#include <stdio.h>
#include <string.h>
+/* Force use of the HLLC Riemann solver */
+#undef RIEMANN_SOLVER_TRRS
+#undef RIEMANN_SOLVER_EXACT
+#undef RIEMANN_SOLVER_HLLC
+#define RIEMANN_SOLVER_HLLC 1
+
/* Local headers. */
#include "riemann/riemann_hllc.h"
#include "swift.h"
diff --git a/tests/testSymmetry.c b/tests/testSymmetry.c
index 1caae903418294439a107ce3a0972f15850ad5d1..1ab493a7c149070dc667a2377ab205df7f873856 100644
--- a/tests/testSymmetry.c
+++ b/tests/testSymmetry.c
@@ -41,7 +41,7 @@ void test() {
/* voronoi_set_box(box_anchor, box_side);*/
#endif
- /* Start with some values for the cosmological paramters */
+ /* Start with some values for the cosmological parameters */
const float a = (float)random_uniform(0.8, 1.);
const float H = 1.f;
@@ -62,7 +62,7 @@ void test() {
pi.time_bin = 1;
pj.time_bin = 1;
-#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
+#if defined(GIZMO_MFV_SPH) || defined(SHADOWFAX_SPH)
/* Give the primitive variables sensible values, since the Riemann solver does
not like negative densities and pressures */
pi.primitives.rho = random_uniform(0.1f, 1.0f);
@@ -106,11 +106,12 @@ void test() {
pj.primitives.gradients.P[0] = 0.0f;
pj.primitives.gradients.P[1] = 0.0f;
pj.primitives.gradients.P[2] = 0.0f;
+
+#ifdef SHADOWFAX_SPH
/* set time step to reasonable value */
pi.force.dt = 0.001;
pj.force.dt = 0.001;
-#ifdef SHADOWFAX_SPH
voronoi_cell_init(&pi.cell, pi.x, box_anchor, box_side);
voronoi_cell_init(&pj.cell, pj.x, box_anchor, box_side);
#endif
@@ -187,33 +188,33 @@ void test() {
runner_iact_nonsym_force(r2, dx, pj2.h, pi2.h, &pj2, &pi2, a, H);
/* Check that the particles are the same */
-#if defined(GIZMO_SPH)
+#if defined(GIZMO_MFV_SPH)
i_not_ok = 0;
j_not_ok = 0;
for (size_t i = 0; i < sizeof(struct part) / sizeof(float); ++i) {
- float a = *(((float *)&pi) + i);
- float b = *(((float *)&pi2) + i);
- float c = *(((float *)&pj) + i);
- float d = *(((float *)&pj2) + i);
+ float aa = *(((float *)&pi) + i);
+ float bb = *(((float *)&pi2) + i);
+ float cc = *(((float *)&pj) + i);
+ float dd = *(((float *)&pj2) + i);
int a_is_b;
- if ((a + b)) {
- a_is_b = (fabs((a - b) / (a + b)) > 1.e-4);
+ if ((aa + bb)) {
+ a_is_b = (fabs((aa - bb) / (aa + bb)) > 1.e-4);
} else {
- a_is_b = !(a == 0.0f);
+ a_is_b = !(aa == 0.0f);
}
int c_is_d;
- if ((c + d)) {
- c_is_d = (fabs((c - d) / (c + d)) > 1.e-4);
+ if ((cc + dd)) {
+ c_is_d = (fabs((cc - dd) / (cc + dd)) > 1.e-4);
} else {
- c_is_d = !(c == 0.0f);
+ c_is_d = !(cc == 0.0f);
}
if (a_is_b) {
- message("%.8e, %.8e, %lu", a, b, i);
+ message("%.8e, %.8e, %lu", aa, bb, i);
}
if (c_is_d) {
- message("%.8e, %.8e, %lu", c, d, i);
+ message("%.8e, %.8e, %lu", cc, dd, i);
}
i_not_ok |= a_is_b;
diff --git a/tests/testUtilities.c b/tests/testUtilities.c
new file mode 100644
index 0000000000000000000000000000000000000000..b835faba9026661361c8828fce5c18beb2b80889
--- /dev/null
+++ b/tests/testUtilities.c
@@ -0,0 +1,75 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2018 Jacob Kegerreis (jacob.kegerreis@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/>.
+ *
+ ******************************************************************************/
+
+#include "swift.h"
+#include "utilities.h"
+
+/**
+ * @brief Test generic utility functions
+ */
+int main() {
+ /// Test find_value_in_monot_incr_array()
+ int n = 100;
+ float array[n];
+ int index;
+ float x;
+
+ // Initialise test array
+ for (int j = 0; j < n; j++) {
+ array[j] = j;
+ }
+
+ // Typical value
+ x = 42.42f;
+ index = find_value_in_monot_incr_array(x, array, n);
+ if (index != 42) {
+ error("Failed with a typical value ");
+ }
+
+ // Value on array element
+ x = 33.f;
+ index = find_value_in_monot_incr_array(x, array, n);
+ if (index != 33) {
+ error("Failed with an array element ");
+ }
+
+ // Value below array
+ x = -123.f;
+ index = find_value_in_monot_incr_array(x, array, n);
+ if (index != -1) {
+ error("Failed with a value below the array ");
+ }
+
+ // Value above array
+ x = 123.f;
+ index = find_value_in_monot_incr_array(x, array, n);
+ if (index != n) {
+ error("Failed with a value above the array ");
+ }
+
+ // Array slice with typical value
+ x = 9.81f;
+ n = 10;
+ index = find_value_in_monot_incr_array(x, array + 5, n);
+ if (index != 4) {
+ error("Failed with an array slice ");
+ }
+
+ return 0;
+}