Commit 49f09458 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'master' into Pressure-Energy-SPH

parents 2b8d1469 8650b65f
......@@ -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
......
......@@ -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])
......
......@@ -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
......@@ -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
......@@ -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 --------------------------------------------
......
......@@ -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 \
......
......@@ -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). */
......
......@@ -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)
......
......@@ -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);
}
}
/**
......
......@@ -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;
......
......@@ -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;
}
......
......@@ -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"
......
......@@ -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 */
......@@ -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 */
......@@ -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 */
......@@ -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 */
......@@ -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 */
/*******************************************************************************
* 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) {