diff --git a/configure.ac b/configure.ac index 3f284a2f6a25fbf8b51aa8f27fd84d11f4e6ff56..24b103480a7bbd2001d0ee473d6f2f7751a6c85f 100644 --- a/configure.ac +++ b/configure.ac @@ -1394,7 +1394,7 @@ AM_CONDITIONAL([HAVESTANDALONEFOF],[test $enable_standalone_fof = "yes"]) # Gravity scheme. AC_ARG_WITH([gravity], [AS_HELP_STRING([--with-gravity=<scheme>], - [Gravity scheme to use @<:@default, with-potential, default: default@:>@] + [Gravity scheme to use @<:@default, with-potential, with-multi-softening default: default@:>@] )], [with_gravity="$withval"], [with_gravity="default"] @@ -1404,6 +1404,9 @@ case "$with_gravity" in with-potential) AC_DEFINE([POTENTIAL_GRAVITY], [1], [Gravity scheme with potential calculation]) ;; + with-multi-softening) + AC_DEFINE([MULTI_SOFTENING_GRAVITY], [1], [Gravity scheme with per-particle softening value calculated from their mass]) + ;; default) AC_DEFINE([DEFAULT_GRAVITY], [1], [Default gravity scheme]) ;; diff --git a/src/debug.c b/src/debug.c index f8509af805fed90cbdbdfd58849cb7b6352062e9..8e3c581163bb404b8013790ef32549503c93fa38 100644 --- a/src/debug.c +++ b/src/debug.c @@ -73,6 +73,8 @@ #include "./gravity/Default/gravity_debug.h" #elif defined(POTENTIAL_GRAVITY) #include "./gravity/Potential/gravity_debug.h" +#elif defined(MULTI_SOFTENING_GRAVITY) +#include "./gravity/MultiSoftening/gravity_debug.h" #else #error "Invalid choice of gravity variant" #endif diff --git a/src/gravity.h b/src/gravity.h index 57cad8eba5f0772f85affd031a450c3ecb4dde59..4e255d2cbaaae012ab1282c82c8f6f15410e6b4b 100644 --- a/src/gravity.h +++ b/src/gravity.h @@ -34,6 +34,9 @@ #elif defined(POTENTIAL_GRAVITY) #include "./gravity/Potential/gravity.h" #define GRAVITY_IMPLEMENTATION "With potential calculation" +#elif defined(MULTI_SOFTENING_GRAVITY) +#include "./gravity/MultiSoftening/gravity.h" +#define GRAVITY_IMPLEMENTATION "With per-particle softening" #else #error "Invalid choice of gravity variant" #endif diff --git a/src/gravity/MultiSoftening/gravity.h b/src/gravity/MultiSoftening/gravity.h new file mode 100644 index 0000000000000000000000000000000000000000..6d1270c952100f3a25202fcdb22be09f9acaa8d9 --- /dev/null +++ b/src/gravity/MultiSoftening/gravity.h @@ -0,0 +1,233 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * 2016 Tom Theuns (tom.theuns@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_DEFAULT_GRAVITY_H +#define SWIFT_DEFAULT_GRAVITY_H + +#include <float.h> + +/* Local includes. */ +#include "cosmology.h" +#include "gravity_properties.h" +#include "kernel_gravity.h" +#include "minmax.h" + +/** + * @brief Returns the mass of a particle + * + * @param gp The particle of interest + */ +__attribute__((always_inline)) INLINE static float gravity_get_mass( + const struct gpart* restrict gp) { + + return gp->mass; +} + +/** + * @brief Returns the current co-moving softening of a particle + * + * @param gp The particle of interest + * @param grav_props The global gravity properties. + */ +__attribute__((always_inline)) INLINE static float gravity_get_softening( + const struct gpart* gp, const struct gravity_props* restrict grav_props) { + + return grav_props->epsilon_cur; +} + +/** + * @brief Add a contribution to this particle's potential. + * + * Here we do nothing as this version does not accumulate potential. + * + * @param gp The particle. + * @param pot The contribution to add. + */ +__attribute__((always_inline)) INLINE static void +gravity_add_comoving_potential(struct gpart* restrict gp, float pot) {} + +/** + * @brief Returns the comoving potential of a particle. + * + * This returns 0 as this flavour of gravity does not compute the + * particles' potential. + * + * @param gp The particle of interest + */ +__attribute__((always_inline)) INLINE static float +gravity_get_comoving_potential(const struct gpart* restrict gp) { + + return 0.f; +} + +/** + * @brief Returns the physical potential of a particle + * + * This returns 0 as this flavour of gravity does not compute the + * particles' potential. + * + * @param gp The particle of interest. + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static float +gravity_get_physical_potential(const struct gpart* restrict gp, + const struct cosmology* cosmo) { + + return 0.f; +} + +/** + * @brief Computes the gravity time-step of a given particle due to self-gravity + * + * We use Gadget-2's type 0 time-step criterion. + * + * @param gp Pointer to the g-particle data. + * @param a_hydro The accelerations coming from the hydro scheme (can be 0). + * @param grav_props Constants used in the gravity scheme. + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static float +gravity_compute_timestep_self(const struct gpart* const gp, + const float a_hydro[3], + const struct gravity_props* restrict grav_props, + const struct cosmology* cosmo) { + + /* Get physical acceleration (gravity contribution) */ + float a_phys_x = gp->a_grav[0] * cosmo->a_factor_grav_accel; + float a_phys_y = gp->a_grav[1] * cosmo->a_factor_grav_accel; + float a_phys_z = gp->a_grav[2] * cosmo->a_factor_grav_accel; + + /* Get physical acceleration (hydro contribution) */ + a_phys_x += a_hydro[0] * cosmo->a_factor_hydro_accel; + a_phys_y += a_hydro[1] * cosmo->a_factor_hydro_accel; + a_phys_z += a_hydro[2] * cosmo->a_factor_hydro_accel; + + const float ac2 = + a_phys_x * a_phys_x + a_phys_y * a_phys_y + a_phys_z * a_phys_z; + + const float ac_inv = (ac2 > 0.f) ? 1.f / sqrtf(ac2) : FLT_MAX; + + const float epsilon = gravity_get_softening(gp, grav_props); + + const float dt = sqrtf(2. * kernel_gravity_softening_plummer_equivalent_inv * + cosmo->a * grav_props->eta * epsilon * ac_inv); + + return dt; +} + +/** + * @brief Prepares a g-particle for the gravity calculation + * + * Zeroes all the relevant arrays in preparation for the sums taking place in + * the variaous tasks + * + * @param gp The particle to act upon + */ +__attribute__((always_inline)) INLINE static void gravity_init_gpart( + struct gpart* gp) { + + /* Zero the acceleration */ + gp->a_grav[0] = 0.f; + gp->a_grav[1] = 0.f; + gp->a_grav[2] = 0.f; + +#ifdef SWIFT_GRAVITY_FORCE_CHECKS + gp->potential_PM = 0.f; + gp->a_grav_PM[0] = 0.f; + gp->a_grav_PM[1] = 0.f; + gp->a_grav_PM[2] = 0.f; +#endif + +#ifdef SWIFT_DEBUG_CHECKS + gp->num_interacted = 0; + gp->initialised = 1; +#endif +} + +/** + * @brief Finishes the gravity calculation. + * + * Multiplies the forces and accelerations by the appropiate constants. + * Applies cosmological correction for periodic BCs. + * + * No need to apply the potential normalisation correction for periodic + * BCs here since we do not compute the potential. + * + * @param gp The particle to act upon + * @param const_G Newton's constant in internal units. + * @param potential_normalisation Term to be added to all the particles. + * @param periodic Are we using periodic BCs? + */ +__attribute__((always_inline)) INLINE static void gravity_end_force( + struct gpart* gp, float const_G, const float potential_normalisation, + const int periodic) { + + /* Let's get physical... */ + gp->a_grav[0] *= const_G; + gp->a_grav[1] *= const_G; + gp->a_grav[2] *= const_G; + +#ifdef SWIFT_GRAVITY_FORCE_CHECKS + gp->potential_PM *= const_G; + gp->a_grav_PM[0] *= const_G; + gp->a_grav_PM[1] *= const_G; + gp->a_grav_PM[2] *= const_G; +#endif + +#ifdef SWIFT_DEBUG_CHECKS + gp->initialised = 0; /* Ready for next step */ +#endif +} + +/** + * @brief Kick the additional variables + * + * @param gp The particle to act upon + * @param dt The time-step for this kick + */ +__attribute__((always_inline)) INLINE static void gravity_kick_extra( + struct gpart* gp, float dt) {} + +/** + * @brief Sets the values to be predicted in the drifts to their values at a + * kick time + * + * @param gp The particle. + */ +__attribute__((always_inline)) INLINE static void +gravity_reset_predicted_values(struct gpart* gp) {} + +/** + * @brief Initialises the g-particles for the first time + * + * This function is called only once just after the ICs have been + * read in to do some conversions. + * + * @param gp The particle to act upon + * @param grav_props The global properties of the gravity calculation. + */ +__attribute__((always_inline)) INLINE static void gravity_first_init_gpart( + struct gpart* gp, const struct gravity_props* grav_props) { + + gp->time_bin = 0; + + gravity_init_gpart(gp); +} + +#endif /* SWIFT_DEFAULT_GRAVITY_H */ diff --git a/src/gravity/MultiSoftening/gravity_debug.h b/src/gravity/MultiSoftening/gravity_debug.h new file mode 100644 index 0000000000000000000000000000000000000000..dce038c58e1769446861bdf6c9a2a44415642c68 --- /dev/null +++ b/src/gravity/MultiSoftening/gravity_debug.h @@ -0,0 +1,35 @@ +/******************************************************************************* + * 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_DEFAULT_GRAVITY_DEBUG_H +#define SWIFT_DEFAULT_GRAVITY_DEBUG_H + +__attribute__((always_inline)) INLINE static void gravity_debug_particle( + const struct gpart* p) { + printf( + "mass=%.3e time_bin=%d\n" + "x=[%.5e,%.5e,%.5e], v_full=[%.5e,%.5e,%.5e], a=[%.5e,%.5e,%.5e]\n", + p->mass, p->time_bin, p->x[0], p->x[1], p->x[2], p->v_full[0], + p->v_full[1], p->v_full[2], p->a_grav[0], p->a_grav[1], p->a_grav[2]); +#ifdef SWIFT_DEBUG_CHECKS + printf("num_interacted=%lld ti_drift=%lld ti_kick=%lld\n", p->num_interacted, + p->ti_drift, p->ti_kick); +#endif +} + +#endif /* SWIFT_DEFAULT_GRAVITY_DEBUG_H */ diff --git a/src/gravity/MultiSoftening/gravity_iact.h b/src/gravity/MultiSoftening/gravity_iact.h new file mode 100644 index 0000000000000000000000000000000000000000..6fce3ddd512018e9ea4be21111c75904c77cb925 --- /dev/null +++ b/src/gravity/MultiSoftening/gravity_iact.h @@ -0,0 +1,333 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk) + * 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_DEFAULT_GRAVITY_IACT_H +#define SWIFT_DEFAULT_GRAVITY_IACT_H + +/* Includes. */ +#include "kernel_gravity.h" +#include "kernel_long_gravity.h" +#include "multipole.h" + +/* Standard headers */ +#include <float.h> + +/** + * @brief Computes the intensity of the force at a point generated by a + * point-mass. + * + * The returned quantity needs to be multiplied by the distance vector to obtain + * the force vector. + * + * @param r2 Square of the distance to the point-mass. + * @param h2 Square of the softening length. + * @param h_inv Inverse of the softening length. + * @param h_inv3 Cube of the inverse of the softening length. + * @param mass Mass of the point-mass. + * @param f_ij (return) The force intensity. + * @param pot_ij (return) The potential. + */ +__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full( + const float r2, const float h2, const float h_inv, const float h_inv3, + const float mass, float *f_ij, float *pot_ij) { + + /* Get the inverse distance */ + const float r_inv = 1.f / sqrtf(r2 + FLT_MIN); + + /* Should we soften ? */ + if (r2 >= h2) { + + /* Get Newtonian gravity */ + *f_ij = mass * r_inv * r_inv * r_inv; + + } else { + + const float r = r2 * r_inv; + const float ui = r * h_inv; + + float W_f_ij; + kernel_grav_force_eval(ui, &W_f_ij); + + /* Get softened gravity */ + *f_ij = mass * h_inv3 * W_f_ij; + } + + /* No potential calculation */ + *pot_ij = 0.f; +} + +/** + * @brief Computes the intensity of the force at a point generated by a + * point-mass truncated for long-distance periodicity. + * + * The returned quantity needs to be multiplied by the distance vector to obtain + * the force vector. + * + * @param r2 Square of the distance to the point-mass. + * @param h2 Square of the softening length. + * @param h_inv Inverse of the softening length. + * @param h_inv3 Cube of the inverse of the softening length. + * @param mass Mass of the point-mass. + * @param r_s_inv Inverse of the mesh smoothing scale. + * @param f_ij (return) The force intensity. + * @param pot_ij (return) The potential. + */ +__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated( + const float r2, const float h2, const float h_inv, const float h_inv3, + const float mass, const float r_s_inv, float *f_ij, float *pot_ij) { + + /* Get the inverse distance */ + const float r_inv = 1.f / sqrtf(r2 + FLT_MIN); + const float r = r2 * r_inv; + + /* Should we soften ? */ + if (r2 >= h2) { + + /* Get Newtonian gravity */ + *f_ij = mass * r_inv * r_inv * r_inv; + + } else { + + const float ui = r * h_inv; + float W_f_ij; + + kernel_grav_force_eval(ui, &W_f_ij); + + /* Get softened gravity */ + *f_ij = mass * h_inv3 * W_f_ij; + } + + /* Get long-range correction */ + const float u_lr = r * r_s_inv; + float corr_f_lr; + kernel_long_grav_force_eval(u_lr, &corr_f_lr); + *f_ij *= corr_f_lr; + + /* No potential calculation */ + *pot_ij = 0.f; +} + +/** + * @brief Computes the forces at a point generated by a multipole. + * + * This assumes M_100 == M_010 == M_001 == 0. + * This uses the quadrupole and trace of the octupole terms only and defaults to + * the monopole if the code is compiled with low-order gravity only. + * + * @param r_x x-component of the distance vector to the multipole. + * @param r_y y-component of the distance vector to the multipole. + * @param r_z z-component of the distance vector to the multipole. + * @param r2 Square of the distance vector to the multipole. + * @param h The softening length. + * @param h_inv Inverse of the softening length. + * @param m The multipole. + * @param f_x (return) The x-component of the acceleration. + * @param f_y (return) The y-component of the acceleration. + * @param f_z (return) The z-component of the acceleration. + * @param pot (return) The potential. + */ +__attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full( + const float r_x, const float r_y, const float r_z, const float r2, + const float h, const float h_inv, const struct multipole *m, float *f_x, + float *f_y, float *f_z, float *pot) { + +/* In the case where the order is < 3, then there is only a monopole term left. + * We can default to the normal P-P interaction with the mass of the multipole + * and its CoM as the "particle" property */ +#if SELF_GRAVITY_MULTIPOLE_ORDER < 3 + + float f_ij, pot_ij; + runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000, + &f_ij, &pot_ij); + *f_x = f_ij * r_x; + *f_y = f_ij * r_y; + *f_z = f_ij * r_z; + +#else + + /* Get the inverse distance */ + const float r_inv = 1.f / sqrtf(r2); + + /* Compute the derivatives of the potential */ + struct potential_derivatives_M2P d; + potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, 0.f, + &d); + + /* 0th order contributions */ + *f_x = m->M_000 * d.D_100; + *f_y = m->M_000 * d.D_010; + *f_z = m->M_000 * d.D_001; + +#if SELF_GRAVITY_MULTIPOLE_ORDER > 0 + + /* 1st order contributions */ + + /* 1st order contributions are all 0 since the dipole is 0 */ + + /* *f_x = m->M_001 * d.D_101 + m->M_010 * d.D_110 + m->M_100 * d.D_200 ; */ + /* *f_y = m->M_001 * d.D_011 + m->M_010 * d.D_020 + m->M_100 * d.D_110 ; */ + /* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */ + +#endif + +#if SELF_GRAVITY_MULTIPOLE_ORDER > 2 + + /* 2nd order contributions */ + *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 + + m->M_101 * d.D_201 + m->M_110 * d.D_210 + m->M_200 * d.D_300; + *f_y += m->M_002 * d.D_012 + m->M_011 * d.D_021 + m->M_020 * d.D_030 + + m->M_101 * d.D_111 + m->M_110 * d.D_120 + m->M_200 * d.D_210; + *f_z += m->M_002 * d.D_003 + m->M_011 * d.D_012 + m->M_020 * d.D_021 + + m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201; + +#endif + +#if SELF_GRAVITY_MULTIPOLE_ORDER > 3 + + /* 3rd order contributions */ + *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 + + m->M_030 * d.D_130 + m->M_102 * d.D_202 + m->M_111 * d.D_211 + + m->M_120 * d.D_220 + m->M_201 * d.D_301 + m->M_210 * d.D_310 + + m->M_300 * d.D_400; + *f_y += m->M_003 * d.D_013 + m->M_012 * d.D_022 + m->M_021 * d.D_031 + + m->M_030 * d.D_040 + m->M_102 * d.D_112 + m->M_111 * d.D_121 + + m->M_120 * d.D_130 + m->M_201 * d.D_211 + m->M_210 * d.D_220 + + m->M_300 * d.D_310; + *f_z += m->M_003 * d.D_004 + m->M_012 * d.D_013 + m->M_021 * d.D_022 + + m->M_030 * d.D_031 + m->M_102 * d.D_103 + m->M_111 * d.D_112 + + m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 + + m->M_300 * d.D_301; + +#endif + + /* Take care of the the sign convention */ + *f_x *= -1.f; + *f_y *= -1.f; + *f_z *= -1.f; +#endif + + /* No potential calculation */ + *pot = 0.f; +} + +/** + * @brief Computes the forces at a point generated by a multipole, truncated for + * long-range periodicity. + * + * This assumes M_100 == M_010 == M_001 == 0. + * This uses the quadrupole term and trace of the octupole terms only and + * defaults to the monopole if the code is compiled with low-order gravity only. + * + * @param r_x x-component of the distance vector to the multipole. + * @param r_y y-component of the distance vector to the multipole. + * @param r_z z-component of the distance vector to the multipole. + * @param r2 Square of the distance vector to the multipole. + * @param h The softening length. + * @param h_inv Inverse of the softening length. + * @param r_s_inv The inverse of the gravity mesh-smoothing scale. + * @param m The multipole. + * @param f_x (return) The x-component of the acceleration. + * @param f_y (return) The y-component of the acceleration. + * @param f_z (return) The z-component of the acceleration. + * @param pot (return) The potential. + */ +__attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated( + const float r_x, const float r_y, const float r_z, const float r2, + const float h, const float h_inv, const float r_s_inv, + const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) { + +/* In the case where the order is < 3, then there is only a monopole term left. + * We can default to the normal P-P interaction with the mass of the multipole + * and its CoM as the "particle" property */ +#if SELF_GRAVITY_MULTIPOLE_ORDER < 3 + + float f_ij, pot_ij; + runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv, + m->M_000, r_s_inv, &f_ij, &pot_ij); + *f_x = f_ij * r_x; + *f_y = f_ij * r_y; + *f_z = f_ij * r_z; + +#else + + /* Get the inverse distance */ + const float r_inv = 1.f / sqrtf(r2); + + /* Compute the derivatives of the potential */ + struct potential_derivatives_M2P d; + potential_derivatives_compute_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1, + r_s_inv, &d); + + /* 0th order contributions */ + *f_x = m->M_000 * d.D_100; + *f_y = m->M_000 * d.D_010; + *f_z = m->M_000 * d.D_001; + +#if SELF_GRAVITY_MULTIPOLE_ORDER > 0 + + /* 1st order contributions */ + + /* 1st order contributions are all 0 since the dipole is 0 */ + + /* *f_x = m->M_001 * d.D_101 + m->M_010 * d.D_110 + m->M_100 * d.D_200 ; */ + /* *f_y = m->M_001 * d.D_011 + m->M_010 * d.D_020 + m->M_100 * d.D_110 ; */ + /* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */ + +#endif + +#if SELF_GRAVITY_MULTIPOLE_ORDER > 2 + + /* 2nd order contributions */ + *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 + + m->M_101 * d.D_201 + m->M_110 * d.D_210 + m->M_200 * d.D_300; + *f_y += m->M_002 * d.D_012 + m->M_011 * d.D_021 + m->M_020 * d.D_030 + + m->M_101 * d.D_111 + m->M_110 * d.D_120 + m->M_200 * d.D_210; + *f_z += m->M_002 * d.D_003 + m->M_011 * d.D_012 + m->M_020 * d.D_021 + + m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201; + +#endif + +#if SELF_GRAVITY_MULTIPOLE_ORDER > 3 + + /* 3rd order contributions */ + *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 + + m->M_030 * d.D_130 + m->M_102 * d.D_202 + m->M_111 * d.D_211 + + m->M_120 * d.D_220 + m->M_201 * d.D_301 + m->M_210 * d.D_310 + + m->M_300 * d.D_400; + *f_y += m->M_003 * d.D_013 + m->M_012 * d.D_022 + m->M_021 * d.D_031 + + m->M_030 * d.D_040 + m->M_102 * d.D_112 + m->M_111 * d.D_121 + + m->M_120 * d.D_130 + m->M_201 * d.D_211 + m->M_210 * d.D_220 + + m->M_300 * d.D_310; + *f_z += m->M_003 * d.D_004 + m->M_012 * d.D_013 + m->M_021 * d.D_022 + + m->M_030 * d.D_031 + m->M_102 * d.D_103 + m->M_111 * d.D_112 + + m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 + + m->M_300 * d.D_301; +#endif + + /* Take care of the the sign convention */ + *f_x *= -1.f; + *f_y *= -1.f; + *f_z *= -1.f; +#endif + + /* No potential calculation */ + *pot = 0.f; +} + +#endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */ diff --git a/src/gravity/MultiSoftening/gravity_io.h b/src/gravity/MultiSoftening/gravity_io.h new file mode 100644 index 0000000000000000000000000000000000000000..2e443e8fdc2479f3f2feff30281dccad9a7b6519 --- /dev/null +++ b/src/gravity/MultiSoftening/gravity_io.h @@ -0,0 +1,122 @@ +/******************************************************************************* + * 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_DEFAULT_GRAVITY_IO_H +#define SWIFT_DEFAULT_GRAVITY_IO_H + +#include "io_properties.h" + +INLINE static void convert_gpart_pos(const struct engine* e, + const struct gpart* gp, double* ret) { + + if (e->s->periodic) { + ret[0] = box_wrap(gp->x[0], 0.0, e->s->dim[0]); + ret[1] = box_wrap(gp->x[1], 0.0, e->s->dim[1]); + ret[2] = box_wrap(gp->x[2], 0.0, e->s->dim[2]); + } else { + ret[0] = gp->x[0]; + ret[1] = gp->x[1]; + ret[2] = gp->x[2]; + } +} + +INLINE static void convert_gpart_vel(const struct engine* e, + const struct gpart* gp, 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, gp->time_bin); + const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin); + + /* Get time-step since the last kick */ + float dt_kick_grav; + 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); + } else { + dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base; + } + + /* Extrapolate the velocites to the current time */ + ret[0] = gp->v_full[0] + gp->a_grav[0] * dt_kick_grav; + ret[1] = gp->v_full[1] + gp->a_grav[1] * dt_kick_grav; + ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav; + + /* Conversion from internal units to peculiar velocities */ + ret[0] *= cosmo->a_inv; + ret[1] *= cosmo->a_inv; + ret[2] *= cosmo->a_inv; +} + +/** + * @brief Specifies which g-particle fields to read from a dataset + * + * @param gparts The g-particle array. + * @param list The list of i/o properties to read. + * @param num_fields The number of i/o fields to read. + */ +INLINE static void darkmatter_read_particles(struct gpart* gparts, + struct io_props* list, + int* num_fields) { + + /* Say how much we want to read */ + *num_fields = 4; + + /* List what we want to read */ + list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY, + UNIT_CONV_LENGTH, gparts, x); + list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY, + UNIT_CONV_SPEED, gparts, v_full); + list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS, + gparts, mass); + list[3] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY, + UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset); +} + +/** + * @brief Specifies which g-particle fields to write to a dataset + * + * @param gparts The g-particle array. + * @param list The list of i/o properties to write. + * @param num_fields The number of i/o fields to write. + */ +INLINE static void darkmatter_write_particles(const struct gpart* gparts, + struct io_props* list, + int* num_fields) { + + /* Say how much we want to write */ + *num_fields = 5; + + /* List what we want to write */ + list[0] = io_make_output_field_convert_gpart( + "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, gparts, convert_gpart_pos); + list[1] = io_make_output_field_convert_gpart( + "Velocities", FLOAT, 3, UNIT_CONV_SPEED, gparts, convert_gpart_vel); + list[2] = + io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, gparts, mass); + list[3] = io_make_output_field("ParticleIDs", ULONGLONG, 1, + UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset); + list[4] = io_make_output_field("GroupIDs", INT, 1, UNIT_CONV_NO_UNITS, gparts, + group_id); +} + +#endif /* SWIFT_DEFAULT_GRAVITY_IO_H */ diff --git a/src/gravity/MultiSoftening/gravity_part.h b/src/gravity/MultiSoftening/gravity_part.h new file mode 100644 index 0000000000000000000000000000000000000000..2c16a34a11d3a7674f5a6c122da780aeaff4a9dc --- /dev/null +++ b/src/gravity/MultiSoftening/gravity_part.h @@ -0,0 +1,83 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@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_DEFAULT_GRAVITY_PART_H +#define SWIFT_DEFAULT_GRAVITY_PART_H + +/* Gravity particle. */ +struct gpart { + + /*! Particle ID. If negative, it is the negative offset of the #part with + which this gpart is linked. */ + long long id_or_neg_offset; + + /*! Particle position. */ + double x[3]; + + /*! Particle velocity. */ + float v_full[3]; + + /*! Particle acceleration. */ + float a_grav[3]; + + /*! Particle mass. */ + float mass; + + /*! Time-step length */ + timebin_t time_bin; + + /*! Type of the #gpart (DM, gas, star, ...) */ + enum part_type type; + + /* Particle group ID and size in the FOF. */ + size_t group_id, group_size; + +#ifdef SWIFT_DEBUG_CHECKS + + /* Numer of gparts this gpart interacted with */ + long long num_interacted; + + /* Time of the last drift */ + integertime_t ti_drift; + + /* Time of the last kick */ + integertime_t ti_kick; + + /* Has this particle been initialised? */ + int initialised; + +#endif + +#ifdef SWIFT_GRAVITY_FORCE_CHECKS + + /*! Acceleration taken from the mesh only */ + float a_grav_PM[3]; + + /*! Potential taken from the mesh only */ + float potential_PM; + + /* Brute-force particle acceleration. */ + double a_grav_exact[3]; + + /* Brute-force particle potential. */ + double potential_exact; +#endif + +} SWIFT_STRUCT_ALIGN; + +#endif /* SWIFT_DEFAULT_GRAVITY_PART_H */ diff --git a/src/gravity_iact.h b/src/gravity_iact.h index 2aaf7219d0851058cba92d4686ed989572dbcaa4..4c8d707fe33d0955809f8cb394d75c3a2f580fb7 100644 --- a/src/gravity_iact.h +++ b/src/gravity_iact.h @@ -32,6 +32,8 @@ #include "./gravity/Default/gravity_iact.h" #elif defined(POTENTIAL_GRAVITY) #include "./gravity/Potential/gravity_iact.h" +#elif defined(MULTI_SOFTENING_GRAVITY) +#include "./gravity/MultiSoftening/gravity_iact.h" #else #error "Invalid choice of gravity variant" #endif diff --git a/src/gravity_io.h b/src/gravity_io.h index 752ee906081d706865e62bae7c8f505e9ca64347..44f392ee06d43228f7aa46faa5c7d9015fde790b 100644 --- a/src/gravity_io.h +++ b/src/gravity_io.h @@ -30,6 +30,8 @@ #include "./gravity/Default/gravity_io.h" #elif defined(POTENTIAL_GRAVITY) #include "./gravity/Potential/gravity_io.h" +#elif defined(MULTI_SOFTENING_GRAVITY) +#include "./gravity/MultiSoftening/gravity_io.h" #else #error "Invalid choice of gravity variant" #endif diff --git a/src/part.h b/src/part.h index 820834d652cf53b3dce06192d7c552ec7681f2c7..440e5e1eb3415e65d7d9cc8b611bb98de107a1b2 100644 --- a/src/part.h +++ b/src/part.h @@ -95,6 +95,8 @@ #include "./gravity/Default/gravity_part.h" #elif defined(POTENTIAL_GRAVITY) #include "./gravity/Potential/gravity_part.h" +#elif defined(MULTI_SOFTENING_GRAVITY) +#include "./gravity/MultiSoftening/gravity_part.h" #else #error "Invalid choice of gravity variant" #endif diff --git a/src/part_type.c b/src/part_type.c index 1f96d4ef1db4b35a92d133e91498ea10ce472c70..75a3f270d06576831833cf75c9bed9b7fde2c68d 100644 --- a/src/part_type.c +++ b/src/part_type.c @@ -20,5 +20,5 @@ /* This object's header. */ #include "part_type.h" -const char* part_type_names[swift_type_count] = {"Gas", "DM", "Dummy", - "Dummy", "Stars", "BH"}; +const char* part_type_names[swift_type_count] = { + "Gas", "DM", "DM_boundary", "Dummy", "Stars", "BH"}; diff --git a/src/part_type.h b/src/part_type.h index 901f47193fa0e72b362c8dce5199a1d0a20526c9..edab07a17080f90bee99507cd5ee3d6ed08054aa 100644 --- a/src/part_type.h +++ b/src/part_type.h @@ -27,6 +27,7 @@ enum part_type { swift_type_gas = 0, swift_type_dark_matter = 1, + swift_type_dark_matter_boundary = 2, swift_type_stars = 4, swift_type_black_hole = 5, swift_type_count