diff --git a/.gitignore b/.gitignore index 5213975f72..94e155a565 100644 --- a/.gitignore +++ b/.gitignore @@ -72,6 +72,7 @@ examples/*/*/restart examples/Cooling/CoolingRates/cooling_rates examples/Cooling/CoolingRates/cooling_element_*.dat examples/Cooling/CoolingRates/cooling_output.dat +examples/Cooling/CoolingSedovBlast_3D/run/* examples/SubgridTests/StellarEvolution/StellarEvolutionSolution* examples/SubgridTests/CosmologicalStellarEvolution/StellarEvolutionSolution* examples/SmallCosmoVolume/SmallCosmoVolume_DM/power_spectra diff --git a/configure.ac b/configure.ac index 82e96c5e18..c642c47694 100644 --- a/configure.ac +++ b/configure.ac @@ -331,6 +331,18 @@ if test "$enable_debugging_checks" = "yes"; then AC_DEFINE([SWIFT_DEBUG_CHECKS],1,[Enable expensive debugging]) fi +# Check if expensive debugging is on. +AC_ARG_ENABLE([magma2-debugging-checks], + [AS_HELP_STRING([--enable-magma2-debugging-checks], + [Activate expensive MAGMA2 debugging checks @<:@yes/no@:>@] + )], + [enable_magma2_debugging_checks="$enableval"], + [enable_magma2_debugging_checks="no"] +) +if test "$enable_magma2_debugging_checks" = "yes"; then + AC_DEFINE([MAGMA2_DEBUG_CHECKS],1,[Enable expensive MAGMA2 debugging]) +fi + # Check if cell graph is on. AC_ARG_ENABLE([cell-graph], [AS_HELP_STRING([--enable-cell-graph], @@ -2223,7 +2235,7 @@ fi # Hydro scheme. AC_ARG_WITH([hydro], [AS_HELP_STRING([--with-hydro=], - [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, phantom, gizmo-mfv, gizmo-mfm, shadowswift, planetary, remix, sphenix, gasoline, anarchy-pu default: sphenix@:>@] + [Hydro dynamics to use @<:@gadget2, minimal, magma2, pressure-entropy, pressure-energy, pressure-energy-monaghan, phantom, gizmo-mfv, gizmo-mfm, shadowswift, planetary, remix, sphenix, gasoline, anarchy-pu default: sphenix@:>@] )], [with_hydro="$withval"], [with_hydro="sphenix"] @@ -2279,6 +2291,9 @@ case "$with_hydro" in gasoline) AC_DEFINE([GASOLINE_SPH], [1], [Gasoline SPH]) ;; + magma2) + AC_DEFINE([MAGMA2_SPH], [1], [MAGMA2 SPH]) + ;; anarchy-du) AC_DEFINE([SPHENIX_SPH], [1], [SPHENIX SPH]) ;; diff --git a/examples/Cooling/CoolingSedovBlast_3D/makeIC.py b/examples/Cooling/CoolingSedovBlast_3D/makeIC.py index 3449d73253..8b9d86b7ce 100644 --- a/examples/Cooling/CoolingSedovBlast_3D/makeIC.py +++ b/examples/Cooling/CoolingSedovBlast_3D/makeIC.py @@ -38,8 +38,8 @@ print(rho0, "cm s^-3") uL = 1.0e3 * pc -uM = 1.99e30 -uv = 1.0e10 +uM = 1.99e43 +uv = 1.0e5 ut = uL / uv uU = uv ** 2 print("ut:", ut / 3.154e7, "yr") diff --git a/examples/Cooling/CoolingSedovBlast_3D/sedov.yml b/examples/Cooling/CoolingSedovBlast_3D/sedov.yml index 352686068f..30ecc151f6 100644 --- a/examples/Cooling/CoolingSedovBlast_3D/sedov.yml +++ b/examples/Cooling/CoolingSedovBlast_3D/sedov.yml @@ -1,6 +1,6 @@ # Define the system of units to use internally. InternalUnitSystem: - UnitMass_in_cgs: 1.99e33 # g + UnitMass_in_cgs: 1.99e43 # g UnitLength_in_cgs: 3.086e21 # cm UnitVelocity_in_cgs: 1.e5 # Centimeters per second UnitCurrent_in_cgs: 1 # Amperes diff --git a/examples/Cosmology/ZeldovichPancake_3D/makeIC.py b/examples/Cosmology/ZeldovichPancake_3D/makeIC.py index bd1366b53c..c5854ce75c 100644 --- a/examples/Cosmology/ZeldovichPancake_3D/makeIC.py +++ b/examples/Cosmology/ZeldovichPancake_3D/makeIC.py @@ -43,7 +43,7 @@ # Some useful variables in h-full units H_0 = 1.0 / Mpc_in_m * 10 ** 5 # h s^-1 -rho_0 = 3.0 * H_0 ** 2 / (8 * math.pi * G_in_SI) # h^2 kg m^-3 +rho_0 = 3.0 * H_0 ** 2 / (8 * pi * G_in_SI) # h^2 kg m^-3 lambda_i = 64.0 / H_0 * 10 ** 5 # h^-1 m (= 64 h^-1 Mpc) x_min = -0.5 * lambda_i x_max = 0.5 * lambda_i diff --git a/examples/HydroTests/KelvinHelmholtzGrowthRate_3D/run.sh b/examples/HydroTests/KelvinHelmholtzGrowthRate_3D/run.sh index 73f81e3da1..574f55b9ea 100755 --- a/examples/HydroTests/KelvinHelmholtzGrowthRate_3D/run.sh +++ b/examples/HydroTests/KelvinHelmholtzGrowthRate_3D/run.sh @@ -14,7 +14,7 @@ then fi # Run SWIFT -../../../swift --hydro --threads=1 kelvinHelmholtzGrowthRate.yml 2>&1 | tee output.log +../../../swift --hydro --threads=16 kelvinHelmholtzGrowthRate.yml 2>&1 | tee output.log # Plot the solution python3 plotSolution.py 100 diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py b/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py index c6ca37e0c6..161a76ee95 100644 --- a/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py +++ b/examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py @@ -71,7 +71,7 @@ def project(data, m_res, property, ylim): x = x[mask] y = y[mask] - np.float64(ylim[0]) - h = data.gas.smoothing_length[mask] + h = data.gas.smoothing_lengths[mask] if property == "density": property = "masses" diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/run.sh b/examples/HydroTests/Rayleigh-Taylor_2D/run.sh old mode 100644 new mode 100755 diff --git a/src/debug.c b/src/debug.c index 2db8f69dff..32db343678 100644 --- a/src/debug.c +++ b/src/debug.c @@ -80,6 +80,8 @@ #include "./hydro/SPHENIX/hydro_debug.h" #elif defined(GASOLINE_SPH) #include "./hydro/Gasoline/hydro_debug.h" +#elif defined(MAGMA2_SPH) +#include "./hydro/MAGMA2/hydro_debug.h" #elif defined(ANARCHY_PU_SPH) #include "./hydro/AnarchyPU/hydro_debug.h" #else diff --git a/src/forcing/roberts_flow_acceleration/forcing.h b/src/forcing/roberts_flow_acceleration/forcing.h index d739e7e3a9..b07cdc666b 100644 --- a/src/forcing/roberts_flow_acceleration/forcing.h +++ b/src/forcing/roberts_flow_acceleration/forcing.h @@ -87,8 +87,14 @@ __attribute__((always_inline)) INLINE static void forcing_terms_apply( /* Effective viscosity from artificial viscosity, as in eq. 100 from * arXiv:1012.1885 */ +#ifndef MAGMA2_SPH const float nu = terms->nu * p->viscosity.alpha * c_s * p->h / (2.f * (hydro_dimension + 2.f)); +#else + /* TODO: MAGMA2 does not track alpha per particle */ + const float nu = terms->nu * const_viscosity_alpha * c_s * p->h / + (2.f * (hydro_dimension + 1.f)); +#endif const float Vz_factor = terms->Vz_factor; const double k0 = (2. * M_PI / L) * terms->kv; diff --git a/src/hydro.h b/src/hydro.h index c82c402e81..8f7b07e554 100644 --- a/src/hydro.h +++ b/src/hydro.h @@ -82,6 +82,10 @@ #include "./hydro/Gasoline/hydro.h" #include "./hydro/Gasoline/hydro_iact.h" #define SPH_IMPLEMENTATION "Gasoline-2 (Wadsley+ 2017)" +#elif defined(MAGMA2_SPH) +#include "./hydro/MAGMA2/hydro.h" +#include "./hydro/MAGMA2/hydro_iact.h" +#define SPH_IMPLEMENTATION "MAGMA2 (Rosswog 2020)" #elif defined(ANARCHY_PU_SPH) #include "./hydro/AnarchyPU/hydro.h" #include "./hydro/AnarchyPU/hydro_iact.h" diff --git a/src/hydro/MAGMA2/hydro.h b/src/hydro/MAGMA2/hydro.h new file mode 100644 index 0000000000..7ffb2790ff --- /dev/null +++ b/src/hydro/MAGMA2/hydro.h @@ -0,0 +1,2369 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) & + * Matthieu Schaller (schaller@strw.leidenuniv.nl) + * 2025 Doug Rennehan (douglas.rennehan@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 . + * + ******************************************************************************/ +#ifndef SWIFT_MAGMA2_HYDRO_H +#define SWIFT_MAGMA2_HYDRO_H + +/** + * @file MAGMA2/hydro.h + * @brief Density-Energy non-conservative implementation of SPH, + * with added MAGMA2 physics (Rosswog 2020) (Non-neighbour loop + * equations) + */ + +#include "adiabatic_index.h" +#include "approx_math.h" +#include "cosmology.h" +#include "dimension.h" +#include "entropy_floor.h" +#include "equation_of_state.h" +#include "hydro_parameters.h" +#include "hydro_properties.h" +#include "hydro_space.h" +#include "kernel_hydro.h" +#include "minmax.h" +#include "pressure_floor.h" + +#include +#include +#include +#include + +/** + * @brief Returns the comoving internal energy of a particle at the last + * time the particle was kicked. + * + * For implementations where the main thermodynamic variable + * is not internal energy, this function computes the internal + * energy from the thermodynamic variable. + * + * @param p The particle of interest + * @param xp The extended data of the particle of interest. + */ +__attribute__((always_inline)) INLINE static float +hydro_get_comoving_internal_energy(const struct part *restrict p, + const struct xpart *restrict xp) { + return xp->u_full; +} + +/** + * @brief Returns the physical internal energy of a particle at the last + * time the particle was kicked. + * + * For implementations where the main thermodynamic variable + * is not internal energy, this function computes the internal + * energy from the thermodynamic variable and converts it to + * physical coordinates. + * + * @param p The particle of interest. + * @param xp The extended data of 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 xpart *restrict xp, + const struct cosmology *cosmo) { + return xp->u_full * cosmo->a_factor_internal_energy; +} + +/** + * @brief Returns the comoving internal energy of a particle drifted to the + * current time. + * + * @param p The particle of interest + */ +__attribute__((always_inline)) INLINE static float +hydro_get_drifted_comoving_internal_energy(const struct part *restrict p) { + return p->u; +} + +/** + * @brief Returns the physical internal energy of a particle drifted to the + * current time. + * + * @param p The particle of interest. + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static float +hydro_get_drifted_physical_internal_energy(const struct part *restrict p, + const struct cosmology *cosmo) { + return p->u * cosmo->a_factor_internal_energy; +} + +/** + * @brief Returns the comoving pressure of a particle + * + * Computes the pressure based on the particle's properties. + * + * @param p The particle of interest + */ +__attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure( + const struct part *restrict p) { + return gas_pressure_from_internal_energy(p->rho, p->u); +} + +/** + * @brief Returns the physical pressure of a particle + * + * Computes the pressure based on the particle's properties and + * convert it to physical coordinates. + * + * @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 * hydro_get_comoving_pressure(p); +} + +/** + * @brief Returns the comoving entropy of a particle at the last + * time the particle was kicked. + * + * For implementations where the main thermodynamic variable + * is not entropy, this function computes the entropy from + * the thermodynamic variable. + * + * @param p The particle of interest + * @param xp The extended data of the particle of interest. + */ +__attribute__((always_inline)) INLINE static float hydro_get_comoving_entropy( + const struct part *restrict p, const struct xpart *restrict xp) { + return gas_entropy_from_internal_energy(p->rho, xp->u_full); +} + +/** + * @brief Returns the physical entropy of a particle at the last + * time the particle was kicked. + * + * For implementations where the main thermodynamic variable + * is not entropy, this function computes the entropy from + * the thermodynamic variable and converts it to + * physical coordinates. + * + * @param p The particle of interest. + * @param xp The extended data of 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 xpart *restrict xp, + const struct cosmology *cosmo) { + /* Note: no cosmological conversion required here with our choice of + * coordinates. */ + return gas_entropy_from_internal_energy(p->rho, xp->u_full); +} + +/** + * @brief Returns the comoving entropy of a particle drifted to the + * current time. + * + * @param p The particle of interest. + */ +__attribute__((always_inline)) INLINE static float +hydro_get_drifted_comoving_entropy(const struct part *restrict p) { + return gas_entropy_from_internal_energy(p->rho, p->u); +} + +/** + * @brief Returns the physical entropy of a particle drifted to the + * current time. + * + * @param p The particle of interest. + * @param cosmo The cosmological model. + */ +__attribute__((always_inline)) INLINE static float +hydro_get_drifted_physical_entropy(const struct part *restrict p, + const struct cosmology *cosmo) { + /* Note: no cosmological conversion required here with our choice of + * coordinates. */ + return gas_entropy_from_internal_energy(p->rho, p->u); +} + +/** + * @brief Returns the comoving 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) { + return gas_soundspeed_from_internal_energy(p->rho, p->u); +} + +/** + * @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 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->rho; +} + +/** + * @brief Returns the comoving 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->rho; +} + +/** + * @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->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->mass = m; +} + +/** + * @brief Returns the time derivative of internal energy of a particle + * + * We assume a constant density. + * + * @param p The particle of interest + */ +__attribute__((always_inline)) INLINE static float +hydro_get_comoving_internal_energy_dt(const struct part *restrict p) { + return p->u_dt; +} + +/** + * @brief Returns the time derivative of internal energy of a particle + * + * We assume a constant density. + * + * @param p The particle of interest + * @param cosmo Cosmology data structure + */ +__attribute__((always_inline)) INLINE static float +hydro_get_physical_internal_energy_dt(const struct part *restrict p, + const struct cosmology *cosmo) { + return p->u_dt * cosmo->a_factor_internal_energy; +} + +/** + * @brief Sets the time derivative of internal energy of a particle + * + * We assume a constant density. + * + * @param p The particle of interest. + * @param du_dt The new time derivative of the internal energy. + */ +__attribute__((always_inline)) INLINE static void +hydro_set_comoving_internal_energy_dt(struct part *restrict p, float du_dt) { + p->u_dt = du_dt; +} + +/** + * @brief Returns the time derivative of internal energy of a particle + * + * We assume a constant density. + * + * @param p The particle of interest. + * @param cosmo Cosmology data structure + * @param du_dt The new time derivative of the internal energy. + */ +__attribute__((always_inline)) INLINE static void +hydro_set_physical_internal_energy_dt(struct part *restrict p, + const struct cosmology *cosmo, + float du_dt) { + p->u_dt = du_dt / cosmo->a_factor_internal_energy; +} + +/** + * @brief Sets the physical entropy of a particle + * + * @param p The particle of interest. + * @param xp The extended particle data. + * @param cosmo Cosmology data structure + * @param entropy The physical entropy + */ +__attribute__((always_inline)) INLINE static void hydro_set_physical_entropy( + struct part *p, struct xpart *xp, const struct cosmology *cosmo, + const float entropy) { + /* Note there is no conversion from physical to comoving entropy */ + const float comoving_entropy = entropy; + xp->u_full = gas_internal_energy_from_entropy(p->rho, comoving_entropy); +} + +/** + * @brief Sets the physical internal energy of a particle + * + * @param p The particle of interest. + * @param xp The extended particle data. + * @param cosmo Cosmology data structure + * @param u The physical internal energy + */ +__attribute__((always_inline)) INLINE static void +hydro_set_physical_internal_energy(struct part *p, struct xpart *xp, + const struct cosmology *cosmo, + const float u) { + xp->u_full = u / cosmo->a_factor_internal_energy; +} + +/** + * @brief Sets the drifted physical internal energy of a particle + * + * @param p The particle of interest. + * @param cosmo Cosmology data structure + * @param pressure_floor The properties of the pressure floor. + * @param u The physical internal energy + */ +__attribute__((always_inline)) INLINE static void +hydro_set_drifted_physical_internal_energy( + struct part *p, const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor, const float u) { + /* There is no need to use the floor here as this function is called in the + * feedback, so the new value of the internal energy should be strictly + * higher than the old value. */ + + p->u = u / cosmo->a_factor_internal_energy; + + /* Now recompute the extra quantities */ + + /* Compute the sound speed */ + const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); + const float pressure_including_floor = + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); + const float soundspeed = + gas_soundspeed_from_pressure(p->rho, pressure_including_floor); + + /* Update variables. */ + p->force.soundspeed = soundspeed; + p->force.pressure = pressure_including_floor; + + /* Signal speed */ + const float v_sig = const_viscosity_alpha_prefactor * soundspeed; + + p->dt_min = min(p->dt_min, p->h_min / v_sig); +} + +/** + * @brief Correct the signal velocity of the particle partaking in + * supernova (kinetic) feedback based on the velocity kick the particle receives + * + * @param p The particle of interest. + * @param cosmo Cosmology data structure + * @param dv_phys The velocity kick received by the particle expressed in + * physical units (note that dv_phys must be positive or equal to zero) + */ +__attribute__((always_inline)) INLINE static void +hydro_set_v_sig_based_on_velocity_kick(struct part *p, + const struct cosmology *cosmo, + const float dv_phys) { + /* Compute the velocity kick in comoving coordinates */ + const float dv = dv_phys / cosmo->a_factor_sound_speed; + + /* Sound speed */ + const float soundspeed = hydro_get_comoving_soundspeed(p); + + /* Signal speed */ + const float v_sig_sound = const_viscosity_alpha_prefactor * soundspeed; + const float v_sig_kick = const_viscosity_beta_prefactor * dv; + const float v_sig = v_sig_sound + v_sig_kick; + + p->dt_min = min(p->dt_min, p->h_min / v_sig); +} + +/** + * @brief Update the value of the viscosity alpha for the scheme. + * + * @param p the particle of interest + * @param alpha the new value for the viscosity coefficient. + */ +__attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha( + struct part *restrict p, float alpha) { } + +/** + * @brief Update the value of the diffusive coefficients to the + * feedback reset value for the scheme. + * + * @param p the particle of interest + */ +__attribute__((always_inline)) INLINE static void +hydro_diffusive_feedback_reset(struct part *restrict p) { + /* Set the viscosity to the max, and the diffusion to the min */ + hydro_set_viscosity_alpha(p, const_viscosity_alpha); +} + +/** + * @brief Computes the hydro time-step of a given particle + * + * This function returns the time-step of a particle given its hydro-dynamical + * state. A typical time-step calculation would be the use of the CFL condition. + * + * @param p Pointer to the particle data + * @param xp Pointer to the extended particle data + * @param hydro_properties The SPH 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) { + + if (p->dt_min == 0.f) return FLT_MAX; + + /* Use full kernel support and physical time */ + const float conv = kernel_gamma * cosmo->a / cosmo->a_factor_sound_speed; + + /* CFL condition */ + const float dt_cfl = 2. * hydro_properties->CFL_condition * conv * p->dt_min; + + /* Do not allow more than 0.25 * |u|/|du/dt| per step */ + const float dt_u = + (p->u_dt_cond != 0.) ? 0.25 * p->u / fabs(p->u_dt_cond) : FLT_MAX; + +#ifdef MAGMA2_DEBUG_CHECKS + if (dt_u < dt_cfl) { + message("dt_u < dt_cfl for pid=%lld u=%g u_dt_cond=%g dt_min=%g conv=%g " + "dt_cfl=%g dt_u=%g", + p->id, p->u, p->u_dt_cond, p->dt_min, conv, dt_cfl, dt_u); + } +#endif + + return fmin(dt_cfl, dt_u); +} + +/** + * @brief Compute the signal velocity between two gas particles + * + * @param dx Comoving vector separating both particles (pi - pj). + * @brief pi The first #part. + * @brief pj The second #part. + * @brief mu_ij The velocity on the axis linking the particles, or zero if the + * particles are moving away from each other, + * @brief beta The non-linear viscosity constant. + */ +__attribute__((always_inline)) INLINE static float hydro_signal_velocity( + const float dx[3], const struct part *restrict pi, + const struct part *restrict pj, const float mu_ij, const float beta) { + + const float ci = pi->force.soundspeed; + const float cj = pj->force.soundspeed; + const float c_ij = 0.5f * (ci + cj); + + const float v_sig_alpha = const_viscosity_alpha_prefactor * c_ij; + const float v_sig_beta = const_viscosity_beta_prefactor * mu_ij; + const float v_sig = v_sig_alpha - fmin(v_sig_beta, 0.f); + + return v_sig; +} + +/** + * @brief Returns the physical velocity divergence. + * + * @brief p The particle + */ +__attribute__((always_inline)) INLINE static float hydro_get_div_v( + const struct part *restrict p) { + + return p->gradients.velocity_tensor[0][0] + + p->gradients.velocity_tensor[1][1] + + p->gradients.velocity_tensor[2][2]; +} + +/** + * @brief Returns the physical velocity divergence. + * + * @brief p The particle + */ +__attribute__((always_inline)) INLINE static float hydro_get_physical_div_v( + const struct part *restrict p, const struct cosmology* cosmo) { + + const float div_v = hydro_get_div_v(p); + return div_v * cosmo->a2_inv + hydro_dimension * cosmo->H; +} + +/** + * @brief Does some extra hydro operations once the actual physical time step + * for the particle is known. + * + * @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) {} + +/** + * @brief Operations performed when a particle gets removed from the + * simulation volume. + * + * @param p The particle. + * @param xp The extended particle data. + * @param time The simulation time. + */ +__attribute__((always_inline)) INLINE static void hydro_remove_part( + const struct part *p, const struct xpart *xp, const double time) {} + +/** + * @brief Prepares a particle for the density calculation. + * + * Zeroes all the relevant arrays in preparation for the sums taking place in + * the various density loop over neighbours. Typically, all fields of the + * density sub-structure of a particle get zeroed in here. + * + * @param p The particle to act upon + * @param hs #hydro_space containing hydro specific space information. + */ +__attribute__((always_inline)) INLINE static void hydro_init_part( + struct part *restrict p, const struct hydro_space *hs) { + p->density.wcount = 0.f; + p->density.wcount_dh = 0.f; + p->rho = 0.f; + p->density.rho_dh = 0.f; +#ifdef MAGMA2_DEBUG_CHECKS + p->debug.num_ngb = 0; + p->debug.v_sig_visc_max = 0.; + p->debug.v_sig_cond_max = 0.; +#endif + +#ifdef hydro_props_use_adiabatic_correction + p->gradients.adiabatic_f_numerator = 0.; +#endif + + p->gradients.wcount = 0.; + p->gradients.u_well_conditioned = 1; + p->gradients.D_well_conditioned = 1; + + p->gradients.du_min = FLT_MAX; + p->gradients.du_max = -FLT_MAX; + p->gradients.kernel_size = FLT_MIN; + + /* These must be zeroed before the density loop */ + for (int i = 0; i < 3; i++) { + p->gradients.u_aux[i] = 0.; + p->gradients.u_aux_norm[i] = 0.; + + p->gradients.dv_min[i] = FLT_MAX; + p->gradients.dv_max[i] = -FLT_MAX; + + for (int j = 0; j < 3; j++) { + p->gradients.velocity_tensor_aux[i][j] = 0.; + p->gradients.velocity_tensor_aux_norm[i][j] = 0.; + } + } +} + +/** + * @brief Computes the condition number of a matrix A + * + * + * @param A The matrix to compute the condition number of. + */ +__attribute__((always_inline)) INLINE static +double condition_number(gsl_matrix *A) { + + gsl_matrix *A_copy = gsl_matrix_alloc(3, 3); + gsl_matrix_memcpy(A_copy, A); + + gsl_vector *S = gsl_vector_alloc(3); + gsl_vector *work = gsl_vector_alloc(3); + gsl_matrix *V = gsl_matrix_alloc(3, 3); + + gsl_linalg_SV_decomp(A_copy, V, S, work); + + double s_max = gsl_vector_get(S, 0); + double s_min = gsl_vector_get(S, 2); + + gsl_matrix_free(A_copy); + gsl_matrix_free(V); + gsl_vector_free(S); + gsl_vector_free(work); + + return (s_min != 0.) ? s_max / s_min : const_condition_number_upper_limit; +} + +/** + * @brief Vector dot product of two 3D vectors. + * + * + * @param vec_a The first vector. + * @param vec_b The second vector. + * @param result The result of the dot product. + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_vec3_vec3_dot(const hydro_real_t *restrict vec_a, + const hydro_real_t *restrict vec_b) { + + return vec_a[0] * vec_b[0] + vec_a[1] * vec_b[1] + vec_a[2] * vec_b[2]; +} + +/** + * @brief Norm of two 3D vectors. + * + * + * @param vec The vector. + * @param result The result of the norm. + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_vec3_norm(const hydro_real_t *restrict vec) { + + return sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]); +} + +/** + * @brief Unit vector of the given vector. + * + * + * @param vec The vector. + * @param result The unit vector of vec + */ +__attribute__((always_inline)) INLINE static +void hydro_vec3_unit(const hydro_real_t *restrict vec, + hydro_real_t *restrict result) { + + result[0] = 0.; + result[1] = 0.; + result[2] = 0.; + + const hydro_real_t vec_norm = hydro_vec3_norm(vec); + if (vec_norm > 0.) { + result[0] = vec[0] / vec_norm; + result[1] = vec[1] / vec_norm; + result[2] = vec[2] / vec_norm; + } +} + +/** + * @brief The Frobenius inner product of two matrices. + * + * + * @param mat The matrix to contract with the vector. + * @param vec The vector to contract with the matrix. + * @param result The result of the contraction. + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_mat3x3_mat3x3_dot(const hydro_real_t (*restrict mat_a)[3], + const hydro_real_t (*restrict mat_b)[3]) { + + return mat_a[0][0] * mat_b[0][0] + + mat_a[0][1] * mat_b[0][1] + + mat_a[0][2] * mat_b[0][2] + + mat_a[1][0] * mat_b[1][0] + + mat_a[1][1] * mat_b[1][1] + + mat_a[1][2] * mat_b[1][2] + + mat_a[2][0] * mat_b[2][0] + + mat_a[2][1] * mat_b[2][1] + + mat_a[2][2] * mat_b[2][2]; +} + +/** + * @brief Contracts the last index of matrix mat with a vector vec and stores in + * result. + * + * + * @param mat The matrix to contract with the vector. + * @param vec The vector to contract with the matrix. + * @param result The result of the contraction. + */ +__attribute__((always_inline)) INLINE static void hydro_mat3x3_vec3_dot( + const hydro_real_t (*restrict mat)[3], + const hydro_real_t *restrict vec, + hydro_real_t *restrict result) { + + result[0] = mat[0][0] * vec[0] + mat[0][1] * vec[1] + mat[0][2] * vec[2]; + result[1] = mat[1][0] * vec[0] + mat[1][1] * vec[1] + mat[1][2] * vec[2]; + result[2] = mat[2][0] * vec[0] + mat[2][1] * vec[1] + mat[2][2] * vec[2]; +} + +/** + * @brief Contracts the last two indices of the tensor tensor with a matrix + * mat and stores in result. Form: mat^T * tensor * mat. + * + * + * @param tensor The tensor to contract with the matrix. + * @param mat The matrix to contract with the tensor. + * @param result The result of the contraction. + */ +__attribute__((always_inline)) INLINE static +void hydro_tensor3x3x3_matrix3x3_dot( + const hydro_real_t (*restrict tensor)[3][3], + const hydro_real_t (*restrict mat)[3], + hydro_real_t *restrict result) { + + result[0] = tensor[0][0][0] * mat[0][0] + + tensor[0][0][1] * mat[0][1] + + tensor[0][0][2] * mat[0][2] + + tensor[0][1][0] * mat[1][0] + + tensor[0][1][1] * mat[1][1] + + tensor[0][1][2] * mat[1][2] + + tensor[0][2][0] * mat[2][0] + + tensor[0][2][1] * mat[2][1] + + tensor[0][2][2] * mat[2][2]; + + result[1] = tensor[1][0][0] * mat[0][0] + + tensor[1][0][1] * mat[0][1] + + tensor[1][0][2] * mat[0][2] + + tensor[1][1][0] * mat[1][0] + + tensor[1][1][1] * mat[1][1] + + tensor[1][1][2] * mat[1][2] + + tensor[1][2][0] * mat[2][0] + + tensor[1][2][1] * mat[2][1] + + tensor[1][2][2] * mat[2][2]; + + result[2] = tensor[2][0][0] * mat[0][0] + + tensor[2][0][1] * mat[0][1] + + tensor[2][0][2] * mat[0][2] + + tensor[2][1][0] * mat[1][0] + + tensor[2][1][1] * mat[1][1] + + tensor[2][1][2] * mat[1][2] + + tensor[2][2][0] * mat[2][0] + + tensor[2][2][1] * mat[2][1] + + tensor[2][2][2] * mat[2][2]; +} + +/** + * @brief Constructs the outer product of two 3D vectors. + * + * + * @param vec_a The first vector. + * @param vec_b The second vector. + * @param result The result of the outer product. + */ +__attribute__((always_inline)) INLINE static void hydro_vec3_vec3_outer( + const hydro_real_t *restrict vec_a, + const hydro_real_t *restrict vec_b, + hydro_real_t (*restrict result)[3]) { + + result[0][0] = vec_a[0] * vec_b[0]; + result[0][1] = vec_a[0] * vec_b[1]; + result[0][2] = vec_a[0] * vec_b[2]; + + result[1][0] = vec_a[1] * vec_b[0]; + result[1][1] = vec_a[1] * vec_b[1]; + result[1][2] = vec_a[1] * vec_b[2]; + + result[2][0] = vec_a[2] * vec_b[0]; + result[2][1] = vec_a[2] * vec_b[1]; + result[2][2] = vec_a[2] * vec_b[2]; +} + +/** + * @brief Limit gradients to variation across the kernel. + * + * + * @param df_raw Raw value + * @param df_reconstructed Reconstructed value + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_scalar_minmod(const hydro_real_t df_raw, + const hydro_real_t df_reconstructed) { + + if (df_raw * df_reconstructed <= 0.) return 0.; + + return (fabs(df_raw) < fabs(df_reconstructed)) ? df_raw : df_reconstructed; +} + +/** + * @brief Limit gradients to variation across the kernel. + * + * + * @param df_reconstructed Reconstructed estimate of df + * @param df_raw Particle estimate of df + * @param df_min_i Minimum value of df_raw across the kernel. + * @param df_max_i Maximum value of df_raw across the kernel. + * @param df_min_j Minimum value of df_raw across the kernel. + * @param df_max_j Maximum value of df_raw across the kernel + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_scalar_minmod_limiter(const hydro_real_t df_reconstructed, + const hydro_real_t df_raw, + const hydro_real_t df_min_i, + const hydro_real_t df_max_i, + const hydro_real_t df_min_j, + const hydro_real_t df_max_j) { + +#ifdef hydro_props_use_strict_minmod_limiter + hydro_real_t df = hydro_scalar_minmod(df_raw, df_reconstructed); + + const hydro_real_t lo = fmax(df_min_i, -df_max_j); + const hydro_real_t hi = fmin(df_max_i, -df_min_j); + + if (lo > hi) return df_raw; + + if (df < lo) df = lo; + if (df > hi) df = hi; + + return df; +#else + return df_reconstructed; +#endif +} + +/** + * @brief Limit variations across the kernel (for 3-vectors) + * + * + * @param dv_reconstructed Reconstructed estimate of df + * @param dv_raw Particle estimate of df + * @param dv_min_i Minimum value of df_raw across the kernel. + * @param dv_max_i Maximum value of df_raw across the kernel. + * @param dv_min_j Minimum value of df_raw across the kernel. + * @param dv_max_j Maximum value of df_raw across the kernel + * @param dv_ij The vector to return. + */ +__attribute__((always_inline)) INLINE static +void hydro_vec_minmod_limiter(const hydro_real_t *restrict dv_reconstructed, + const hydro_real_t *restrict dv_raw, + const hydro_real_t *restrict dv_min_i, + const hydro_real_t *restrict dv_max_i, + const hydro_real_t *restrict dv_min_j, + const hydro_real_t *restrict dv_max_j, + hydro_real_t *restrict dv_ij) { + + dv_ij[0] = hydro_scalar_minmod_limiter(dv_reconstructed[0], dv_raw[0], + dv_min_i[0], dv_max_i[0], + dv_min_j[0], dv_max_j[0]); + dv_ij[1] = hydro_scalar_minmod_limiter(dv_reconstructed[1], dv_raw[1], + dv_min_i[1], dv_max_i[1], + dv_min_j[1], dv_max_j[1]); + dv_ij[2] = hydro_scalar_minmod_limiter(dv_reconstructed[2], dv_raw[2], + dv_min_i[2], dv_max_i[2], + dv_min_j[2], dv_max_j[2]); + + /* If any of the components are exactly zero, we return a zero vector */ + if (dv_ij[0] == 0. || dv_ij[1] == 0. || dv_ij[2] == 0.) { + dv_ij[0] = 0.; + dv_ij[1] = 0.; + dv_ij[2] = 0.; + } +} + +/** + * @brief Limit gradients to variation across the kernel. + * + * + * @param df_min Minimum value of delta f across the kernel. + * @param df_max Maximum value of delta f across the kernel. + * @param kernel_size Interaction distance across the kernel. + * @param grad The vector gradient to slope limit. + */ +__attribute__((always_inline)) INLINE static +void hydro_vec_slope_limiter(const hydro_real_t df_min, + const hydro_real_t df_max, + const hydro_real_t kernel_size, + hydro_real_t *restrict grad) { + +#ifdef hydro_props_use_extra_slope_limiter + const hydro_real_t grad_norm = hydro_vec3_norm(grad); + const hydro_real_t length = const_grad_overshoot_length * kernel_size; + + /* Nothing to do if there is no gradient or no look-ahead distance */ + if (grad_norm > 0. && length > 0.) { + const hydro_real_t df_min_abs = (df_min < 0.) ? fabs(df_min) : 0.; + const hydro_real_t df_max_abs = (df_max > 0.) ? fabs(df_max) : 0.; + const hydro_real_t df_abs_min = fmin(df_min_abs, df_max_abs); + + hydro_real_t bound = df_abs_min; + +#ifdef const_grad_overshoot_tolerance + const hydro_real_t tolerance = const_grad_overshoot_tolerance; + if (tolerance > 0.) { + const hydro_real_t df_abs_max = fmax(df_min_abs, df_max_abs); + const hydro_real_t extra = tolerance * df_abs_max; + const hydro_real_t cap = + (df_abs_min + extra < df_abs_max) ? df_abs_min + extra : df_abs_max; + bound = cap; + } +#endif + + const hydro_real_t limiter = + (bound > 0.) ? (bound / (length * grad_norm)) : 0.; + + if (limiter < 1.) { + grad[0] *= limiter; + grad[1] *= limiter; + grad[2] *= limiter; + } + } +#endif + +} + +/** + * @brief Computes Phi_ab from Rosswog 2020 21 for a field given A. This is the + * van Leer 1974 slope limiting procedure. + * + * + * @param A_ij The ratio of the gradients of the two particles. + * @param eta_i The normed smoothing length of the first particle. + * @param eta_j The normed smoothing length of the second particle. + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_van_leer_phi(const hydro_real_t A_ij, + const hydro_real_t eta_i, + const hydro_real_t eta_j) { + + hydro_real_t phi_raw = (4. * A_ij) / ((1. + A_ij) * (1. + A_ij)); + phi_raw = fmin(1., phi_raw); + phi_raw = fmax(0., phi_raw); + + /* η_ab and η_crit damping */ + const hydro_real_t eta_ij = min(eta_i, eta_j); + + hydro_real_t damping = 1.; + if (eta_ij <= const_slope_limiter_eta_crit) { + const hydro_real_t diff = + (eta_ij - const_slope_limiter_eta_crit) / const_slope_limiter_eta_fold; + damping = exp(-diff * diff); + } + + phi_raw *= damping; + + /* Handle any edge cases */ + phi_raw = fmin(1., phi_raw); + phi_raw = fmax(0., phi_raw); + + return phi_raw; +} + +/** + * @brief Computes A_ab from Rosswog 2020 Eq. 22 for a scalar field. + * + * + * @param grad_i The gradient of the quantity for the first particle. + * @param grad_j The gradient of the quantity for the second particle. + * @param dx The distance vector between the two particles ( ri - rj ). + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_scalar_van_leer_A(const hydro_real_t *restrict grad_i, + const hydro_real_t *restrict grad_j, + const hydro_real_t *restrict dx) { + + const hydro_real_t grad_dot_x_i = hydro_vec3_vec3_dot(grad_i, dx); + const hydro_real_t grad_dot_x_j = hydro_vec3_vec3_dot(grad_j, dx); + + /* Fall-back to raw estimates */ + if (grad_dot_x_i * grad_dot_x_j <= 0.) return 0.; + + /* Regularize denominator */ + if (fabs(grad_dot_x_j) < 1.e-10) return 0.; + + const hydro_real_t A_ij = grad_dot_x_i / grad_dot_x_j; + + return A_ij; +} + +/** + * @brief Computes A_ab from Rosswog 2020 Eq. 22 for a vector field. + * + * + * @param grad_i The gradient tensor for the first particle. + * @param grad_j The gradient tensor for the second particle. + * @param dx The distance vector between the two particles ( ri - rj ). + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_vector_van_leer_A(const hydro_real_t (*restrict grad_i)[3], + const hydro_real_t (*restrict grad_j)[3], + const hydro_real_t *restrict dx) { + + hydro_real_t delta_ij[3][3] = {0}; + hydro_vec3_vec3_outer(dx, dx, delta_ij); + + const hydro_real_t grad_dot_x_x_i = hydro_mat3x3_mat3x3_dot(grad_i, delta_ij); + const hydro_real_t grad_dot_x_x_j = hydro_mat3x3_mat3x3_dot(grad_j, delta_ij); + + /* Fall-back to raw estimates */ + if (grad_dot_x_x_i * grad_dot_x_x_j <= 0.) return 0.; + + /* Regularize denominator */ + if (fabs(grad_dot_x_x_j) < 1.e-10) return 0.; + + const hydro_real_t A_ij = grad_dot_x_x_i / grad_dot_x_x_j; + + return A_ij; +} + +/** + * @brief Computes Phi_ab from Rosswog 2020 21 for a scalar field. This is the + * van Leer 1974 slope limiting procedure. + * + * + * @param grad_i The gradient of the quantity for the first particle. + * @param grad_j The gradient of the quantity for the second particle. + * @param dx The distance vector between the two particles ( ri - rj ). + * @param eta_i The normed smoothing length of the first particle. + * @param eta_j The normed smoothing length of the second particle. + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_scalar_van_leer_phi(const hydro_real_t *restrict grad_i, + const hydro_real_t *restrict grad_j, + const hydro_real_t *restrict dx, + const hydro_real_t eta_i, + const hydro_real_t eta_j) { + + const hydro_real_t A_ij = hydro_scalar_van_leer_A(grad_i, grad_j, dx); + + return hydro_van_leer_phi(A_ij, eta_i, eta_j); +} + +/** + * @brief Computes Phi_ab from Rosswog 2020 21 for a vector field. This is the + * van Leer 1974 slope limiting procedure. + * + * + * @param grad_i The gradient of the quantity for the first particle. + * @param grad_j The gradient of the quantity for the second particle. + * @param dx The distance vector between the two particles ( ri - rj ). + * @param eta_i The normed smoothing length of the first particle. + * @param eta_j The normed smoothing length of the second particle. + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_vector_van_leer_phi(const hydro_real_t (*restrict grad_i)[3], + const hydro_real_t (*restrict grad_j)[3], + const hydro_real_t *restrict dx, + const hydro_real_t eta_i, + const hydro_real_t eta_j) { + + const hydro_real_t A_ij = hydro_vector_van_leer_A(grad_i, grad_j, dx); + + return hydro_van_leer_phi(A_ij, eta_i, eta_j); +} + +/** + * @brief Reconstructs the scalar field at the mid-point between to the + * two particles to second order. + * + * + * @param phi The slope limiter value from the van Leer 1974 scheme. + * @param dx The distance vector between the two particles ( ri - rj ). + * @param f The field at the particle. + * @param grad The gradient of the field at the particle. + * @param hess The Hessian of the field at the particle. + * @param f_reconstructed The reconstructed field at the mid-point (2nd order). + */ +__attribute__((always_inline)) INLINE static void +hydro_scalar_second_order_reconstruction(const hydro_real_t phi, + const hydro_real_t *restrict dx, + const hydro_real_t f, + const hydro_real_t *restrict grad, + const hydro_real_t (*restrict hess)[3], + hydro_real_t *f_reconstructed) { + + /* Midpoint from Equation 17 Rosswog 2020 */ + const hydro_real_t midpoint[3] = {0.5 * dx[0], + 0.5 * dx[1], + 0.5 * dx[2]}; + hydro_real_t delta_ij[3][3] = {0}; + hydro_vec3_vec3_outer(midpoint, midpoint, delta_ij); + + const hydro_real_t df_dx_dot_r = hydro_vec3_vec3_dot(grad, midpoint); + const hydro_real_t df_dx_dx_dot_r2 = hydro_mat3x3_mat3x3_dot(hess, delta_ij); + + /* Apply limited slope reconstruction */ + *f_reconstructed = f + phi * (df_dx_dot_r + 0.5 * df_dx_dx_dot_r2); +} + +/** + * @brief Reconstructs the vector field at the mid-point between to the + * two particles to second order. + * + * + * @param phi The slope limiter value from the van Leer 1974 scheme. + * @param dx The distance vector between the two particles ( ri - rj ). + * @param f The vector field at the particle. + * @param grad The gradient of the vector field at the particle. + * @param hess The Hessian of the vector field at the particle. + * @param f_reconstructed The reconstructed vector field at the mid-point (2nd order). + */ +__attribute__((always_inline)) INLINE static void +hydro_vector_second_order_reconstruction(const hydro_real_t phi, + const hydro_real_t *restrict dx, + const hydro_real_t *restrict f, + const hydro_real_t (*restrict grad)[3], + const hydro_real_t (*restrict hess)[3][3], + hydro_real_t *restrict f_reconstructed) { + + /* Midpoint from Equation 17 Rosswog 2020 */ + const hydro_real_t midpoint[3] = {0.5 * dx[0], + 0.5 * dx[1], + 0.5 * dx[2]}; + hydro_real_t delta_ij[3][3] = {0}; + hydro_vec3_vec3_outer(midpoint, midpoint, delta_ij); + + hydro_real_t df_dx_dot_r[3] = {0}; + hydro_real_t df_dx_dx_dot_r2[3] = {0}; + + hydro_mat3x3_vec3_dot(grad, midpoint, df_dx_dot_r); + hydro_tensor3x3x3_matrix3x3_dot(hess, delta_ij, df_dx_dx_dot_r2); + + /* Apply limited slope reconstruction */ + f_reconstructed[0] = f[0] + phi * (df_dx_dot_r[0] + 0.5 * df_dx_dx_dot_r2[0]); + f_reconstructed[1] = f[1] + phi * (df_dx_dot_r[1] + 0.5 * df_dx_dx_dot_r2[1]); + f_reconstructed[2] = f[2] + phi * (df_dx_dot_r[2] + 0.5 * df_dx_dx_dot_r2[2]); +} + +/** + * @brief Get the balsara limiter for viscosity. + * + * + * @param p Particle p + * @param cosmo The cosmology structure + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_get_balsara_limiter(const struct part *restrict p, + const struct cosmology* cosmo) { + +#ifdef hydro_props_use_balsara_limiter + const hydro_real_t fac_B = cosmo->a_factor_Balsara_eps; + hydro_real_t balsara = 1.; + + /* Can't trust velocity_tensor when having ill-conditioned matrices */ + if (p->gradients.C_well_conditioned && p->gradients.D_well_conditioned) { + const hydro_real_t div_v_phys = hydro_get_physical_div_v(p, cosmo); + + const hydro_real_t curl_v[3] = { + p->gradients.velocity_tensor[2][1] - p->gradients.velocity_tensor[1][2], + p->gradients.velocity_tensor[0][2] - p->gradients.velocity_tensor[2][0], + p->gradients.velocity_tensor[1][0] - p->gradients.velocity_tensor[0][1] + }; + + const hydro_real_t curl_v_norm_phys = + hydro_vec3_norm(curl_v) * cosmo->a2_inv; + const hydro_real_t Balsara_eps = 1.e-4 * fac_B * p->force.soundspeed / p->h; + balsara = + fabs(div_v_phys) / (fabs(div_v_phys) + curl_v_norm_phys + Balsara_eps); + balsara = min(1., balsara); + balsara = max(0., balsara); + } + + return balsara; +#else + return 1.; +#endif +} + +/** + * @brief Computes the average kernel gradient between pi and pj + * + * + * @param pi Particle pi + * @param pj Particle pj + * @param G_i Kernel gradient at pi + * @param G_j Kernel gradient at pj + * @param G_ij Kernel gradient average to fill + */ +__attribute__((always_inline)) INLINE static +void hydro_get_average_kernel_gradient(const struct part *restrict pi, + const struct part *restrict pj, + const hydro_real_t *restrict G_i, + const hydro_real_t *restrict G_j, + hydro_real_t *restrict G_ij) { + +#if (hydro_props_kernel_gradient_weighting == 0) + G_ij[0] = 0.5 * (G_i[0] + G_j[0]); + G_ij[1] = 0.5 * (G_i[1] + G_j[1]); + G_ij[2] = 0.5 * (G_i[2] + G_j[2]); +#elif (hydro_props_kernel_gradient_weighting == 1) + const hydro_real_t f_i = pi->force.f; + const hydro_real_t f_j = pj->force.f; + G_ij[0] = 0.5 * (f_i * G_i[0] + f_j * G_j[0]); + G_ij[1] = 0.5 * (f_i * G_i[1] + f_j * G_j[1]); + G_ij[2] = 0.5 * (f_i * G_i[2] + f_j * G_j[2]); +#elif (hydro_props_kernel_gradient_weighting == 2) + const hydro_real_t f_i = pi->force.f; + const hydro_real_t f_j = pj->force.f; + const hydro_real_t f_ij = 0.5 * (f_i + f_j); + G_ij[0] = 0.5 * f_ij * (G_i[0] + G_j[0]); + G_ij[1] = 0.5 * f_ij * (G_i[1] + G_j[1]); + G_ij[2] = 0.5 * f_ij * (G_i[2] + G_j[2]); +#elif (hydro_props_kernel_gradient_weighting == 3) + const hydro_real_t f_i = pi->force.f; + const hydro_real_t f_j = pj->force.f; + hydro_real_t f_ij = 0.; + const hydro_real_t f_sum = f_i + f_j; + if (f_sum > 0.) f_ij = 2. * f_i * f_j / f_sum; + + G_ij[0] = 0.5 * f_ij * (G_i[0] + G_j[0]); + G_ij[1] = 0.5 * f_ij * (G_i[1] + G_j[1]); + G_ij[2] = 0.5 * f_ij * (G_i[2] + G_j[2]); +#elif (hydro_props_kernel_gradient_weighting == 4) + const hydro_real_t f_i = pi->force.f; + const hydro_real_t f_j = pj->force.f; + const hydro_real_t f_ij = sqrt(f_i * f_j); + G_ij[0] = 0.5 * f_ij * (G_i[0] + G_j[0]); + G_ij[1] = 0.5 * f_ij * (G_i[1] + G_j[1]); + G_ij[2] = 0.5 * f_ij * (G_i[2] + G_j[2]); +#elif (hydro_props_kernel_gradient_weighting == 5) + const hydro_real_t f_i = pi->force.f; + const hydro_real_t f_j = pj->force.f; + const hydro_real_t rho_sum = pi->rho + pj->rho; + const hydro_real_t f_ij = (pi->rho * f_i + pj->rho * f_j) / rho_sum; + G_ij[0] = 0.5 * f_ij * (G_i[0] + G_j[0]); + G_ij[1] = 0.5 * f_ij * (G_i[1] + G_j[1]); + G_ij[2] = 0.5 * f_ij * (G_i[2] + G_j[2]); +#elif (hydro_props_kernel_gradient_weighting == 6) + const hydro_real_t f_i = pi->force.f; + const hydro_real_t f_j = pj->force.f; + const hydro_real_t rho_f_sum = pi->rho * f_i + pj->rho * f_j; + const hydro_real_t f_ij = 2. * pi->rho * f_i * pj->rho * f_j / rho_f_sum; + G_ij[0] = 0.5 * f_ij * (G_i[0] + G_j[0]); + G_ij[1] = 0.5 * f_ij * (G_i[1] + G_j[1]); + G_ij[2] = 0.5 * f_ij * (G_i[2] + G_j[2]); +#elif (hydro_props_kernel_gradient_weighting == 7) + const hydro_real_t f_i = pi->force.f; + const hydro_real_t f_j = pj->force.f; + const hydro_real_t P_sum = pi->force.pressure + pj->force.pressure; + const hydro_real_t f_ij = + (pi->force.pressure * f_i + pj->force.pressure * f_j) / P_sum; + G_ij[0] = 0.5 * f_ij * (G_i[0] + G_j[0]); + G_ij[1] = 0.5 * f_ij * (G_i[1] + G_j[1]); + G_ij[2] = 0.5 * f_ij * (G_i[2] + G_j[2]); +#elif (hydro_props_kernel_gradient_weighting == 8) + const hydro_real_t f_i = pi->force.f; + const hydro_real_t f_j = pj->force.f; + const hydro_real_t P_f_sum = + pi->force.pressure * f_i + pj->force.pressure * f_j; + const hydro_real_t f_ij = + 2. * pi->force.pressure * f_i * pj->force.pressure * f_j / P_f_sum; + G_ij[0] = 0.5 * f_ij * (G_i[0] + G_j[0]); + G_ij[1] = 0.5 * f_ij * (G_i[1] + G_j[1]); + G_ij[2] = 0.5 * f_ij * (G_i[2] + G_j[2]); +#elif (hydro_props_kernel_gradient_weighting == 9) + const hydro_real_t f_i = pi->force.f; + const hydro_real_t f_j = pj->force.f; + const hydro_real_t mi = hydro_get_mass(pi); + const hydro_real_t mj = hydro_get_mass(pj); + const hydro_real_t rhoi_inv = 1. / hydro_get_comoving_density(pi); + const hydro_real_t rhoj_inv = 1. / hydro_get_comoving_density(pj); + const hydro_real_t Vi = mi * rhoi_inv; + const hydro_real_t Vj = mj * rhoj_inv; + + const hydro_real_t f_ij = 0.5 * (Vi * f_i + Vj * f_j) / (Vi + Vj); + G_ij[0] = 0.5 * f_ij * (G_i[0] + G_j[0]); + G_ij[1] = 0.5 * f_ij * (G_i[1] + G_j[1]); + G_ij[2] = 0.5 * f_ij * (G_i[2] + G_j[2]); +#else + error("Invalid hydro_props_kernel_gradient_weighting value: %d", + hydro_props_kernel_gradient_weighting); +#endif +} + +/** + * @brief Computes the viscosity acceleration and mu_ij terms. + * + * + * @param pi Particle pi + * @param pj Particle pj + * @param dv Second order reconstructed velocity difference + * @param dx Distance vector + * @param r Distance vector norm + * @param fac_mu Cosmological factor for mu_ij + * @param a2_Hubble Hubble flow + * @param mu_i The velocity jump indicator at pi to fill + * @param mu_j The velocity jump indicator at pj to fill + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_get_visc_acc_term(const struct part *restrict pi, + const struct part *restrict pj, + const hydro_real_t *restrict dv, + const hydro_real_t *restrict dx, + const hydro_real_t fac_mu, + const hydro_real_t a2_Hubble) { + + const hydro_real_t rhoi = pi->rho; + const hydro_real_t rhoj = pj->rho; + + const hydro_real_t bi = pi->gradients.balsara; + const hydro_real_t bj = pj->gradients.balsara; + + /* Viscosity direction */ + hydro_real_t dx_hat[3] = {0., 0., 0.}; + hydro_vec3_unit(dx, dx_hat); + + /* Separation and velocity along the gradient vector */ + const hydro_real_t dv_Hubble[3] = {dv[0] + a2_Hubble * dx[0], + dv[1] + a2_Hubble * dx[1], + dv[2] + a2_Hubble * dx[2]}; + const hydro_real_t dv_dot_dx_hat = hydro_vec3_vec3_dot(dv_Hubble, dx_hat); + const hydro_real_t conv = (dv_dot_dx_hat < 0.) ? fac_mu : 0.; + + /* Must be a converging flow */ + if (conv == 0.) return 0.; + + const hydro_real_t mu_ij = conv * dv_dot_dx_hat; + +#ifdef hydro_props_viscosity_weighting_type +#if (hydro_props_viscosity_weighting_type == 0) + + /* Each particle gets its own Q and then weighted by density */ + const hydro_real_t rhoij_inv = 1.0 / (rhoi * rhoj); + + const hydro_real_t q_i_alpha = + -const_viscosity_alpha * pi->force.soundspeed * mu_ij; + const hydro_real_t q_i_beta = const_viscosity_beta * mu_ij * mu_ij; + const hydro_real_t Q_i = rhoi * (q_i_alpha + q_i_beta); + + const hydro_real_t q_j_alpha = + -const_viscosity_alpha * pj->force.soundspeed * mu_ij; + const hydro_real_t q_j_beta = const_viscosity_beta * mu_ij * mu_ij; + const hydro_real_t Q_j = rhoj * (q_j_alpha + q_j_beta); + + return (bi * Q_i + bj * Q_j) * rhoij_inv; + +#elif (hydro_props_viscosity_weighting_type == 1) + + /* Each particle has the same Q but is density weighted */ + const hydro_real_t b_ij = 0.5 * (bi + bj); + const hydro_real_t rhoij_inv = 1. / (rhoi * rhoj); + const hydro_real_t c_ij = 0.5 * (pi->force.soundspeed + pj->force.soundspeed); + const hydro_real_t q_ij_alpha = -const_viscosity_alpha * c_ij * mu_ij; + const hydro_real_t q_ij_beta = const_viscosity_beta * mu_ij * mu_ij; + const hydro_real_t Q_ij = (rhoi + rhoj) * (q_ij_alpha + q_ij_beta); + + return b_ij * Q_ij * rhoij_inv; + +#elif (hydro_props_viscosity_weighting_type == 2) + + /* Particles average symmetrically and arithmetically */ + const hydro_real_t b_ij = 0.5 * (bi + bj); + const hydro_real_t c_ij = 0.5 * (pi->force.soundspeed + pj->force.soundspeed); + const hydro_real_t q_ij_alpha = -const_viscosity_alpha * c_ij * mu_ij; + const hydro_real_t q_ij_beta = const_viscosity_beta * mu_ij * mu_ij; + const hydro_real_t q_ij = 2. * (q_ij_alpha + q_ij_beta) / (rhoi + rhoj); + + return b_ij * q_ij; + +#endif +#else + error("Unknown compiled hydro_props_viscosity_weighting_type value: %d\n" + "Valid values are 0, 1, or 2.\n", + (int)hydro_props_viscosity_weighting_type); +#endif +} + +/** + * @brief Gets the signal velocity for the conduction interaction based on mu_ij. + * + * + * @param dx The separation vector + * @param pi Particle pi + * @param pj Particle pj + * @param mu_i The velocity jump indicator at pi + * @param mu_j The velocity jump indicator at pj + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_get_cond_signal_velocity(const float *restrict dx, + const struct part *restrict pi, + const struct part *restrict pj, + const float mu_i, + const float mu_j, + const float beta) { + + + const float mu_ij = 0.5f * (mu_i + mu_j); + return hydro_signal_velocity(dx, pi, pj, mu_ij, beta); +} + +/** + * @brief Gets the signal velocity for the viscous interaction based on mu_ij. + * + * + * @param dx The separation vector + * @param pi Particle pi + * @param pj Particle pj + * @param mu_i The velocity jump indicator at pi + * @param mu_j The velocity jump indicator at pj + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_get_visc_signal_velocity(const float *restrict dx, + const struct part *restrict pi, + const struct part *restrict pj, + const float mu_i, + const float mu_j, + const float beta) { + +#ifdef hydro_props_viscosity_weighting_type +#if (hydro_props_viscosity_weighting_type == 0) + const float dx_ji[3] = {-dx[0], -dx[1], -dx[2]}; + const float v_sig_visc_i = + hydro_signal_velocity(dx, pi, pj, mu_i, beta); + const float v_sig_visc_j = + hydro_signal_velocity(dx_ji, pj, pi, mu_j, beta); + return 0.5 * (v_sig_visc_i + v_sig_visc_j); +#else + const float mu_ij = 0.5f * (mu_i + mu_j); + return hydro_signal_velocity(dx, pi, pj, mu_ij, beta); +#endif +#else + error("No viscosity weighting type defined at compile time."); +#endif +} + +/** + * @brief Finishes the density calculation. + * + * Multiplies the density and number of neighbours by the appropiate constants + * and add the self-contribution term. + * Additional quantities such as velocity gradients will also get the final + * terms added to them here. + * + * Also adds/multiplies the cosmological terms if need be. + * + * @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 h_inv = 1.0f / h; /* 1/h */ + const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ + const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */ + + /* Final operation on the density (add self-contribution). */ + p->rho += p->mass * kernel_root; + p->density.rho_dh -= hydro_dimension * p->mass * kernel_root; + p->density.wcount += kernel_root; + p->density.wcount_dh -= hydro_dimension * kernel_root; + + /* Finish the calculation by inserting the missing h-factors */ + p->rho *= h_inv_dim; + p->density.rho_dh *= h_inv_dim_plus_one; + p->density.wcount *= h_inv_dim; + p->density.wcount_dh *= h_inv_dim_plus_one; + + /* Need this for correct dh/dt */ + p->gradients.wcount = p->density.wcount; + + p->gradients.u_aux[0] *= h_inv_dim_plus_one; + p->gradients.u_aux[1] *= h_inv_dim_plus_one; + p->gradients.u_aux[2] *= h_inv_dim_plus_one; + + p->gradients.u_aux_norm[0] *= h_inv_dim_plus_one; + p->gradients.u_aux_norm[1] *= h_inv_dim_plus_one; + p->gradients.u_aux_norm[2] *= h_inv_dim_plus_one; + + if (p->gradients.u_aux_norm[0] != 0. && + p->gradients.u_aux_norm[1] != 0. && + p->gradients.u_aux_norm[2] != 0.) { + p->gradients.u_aux[0] /= p->gradients.u_aux_norm[0]; + p->gradients.u_aux[1] /= p->gradients.u_aux_norm[1]; + p->gradients.u_aux[2] /= p->gradients.u_aux_norm[2]; + + hydro_vec_slope_limiter(p->gradients.du_min, p->gradients.du_max, + p->gradients.kernel_size, + p->gradients.u_aux); + } + else { + p->gradients.u_well_conditioned = 0; +#ifdef MAGMA2_DEBUG_CHECKS + p->debug.u_aux[0] = p->gradients.u_aux[0]; + p->debug.u_aux[1] = p->gradients.u_aux[1]; + p->debug.u_aux[2] = p->gradients.u_aux[2]; + + p->debug.u_aux_norm[0] = p->gradients.u_aux_norm[0]; + p->debug.u_aux_norm[1] = p->gradients.u_aux_norm[1]; + p->debug.u_aux_norm[2] = p->gradients.u_aux_norm[2]; + + p->debug.u_ill_conditioned_count++; +#endif + + p->gradients.u_aux[0] = 0.; + p->gradients.u_aux[1] = 0.; + p->gradients.u_aux[2] = 0.; + + p->gradients.u_aux_norm[0] = 0.; + p->gradients.u_aux_norm[1] = 0.; + p->gradients.u_aux_norm[2] = 0.; + } + + double aux_norm[3][3] = {0}; + for (int i = 0; i < 3; i++) { + p->gradients.velocity_tensor_aux[i][0] *= h_inv_dim_plus_one; + p->gradients.velocity_tensor_aux[i][1] *= h_inv_dim_plus_one; + p->gradients.velocity_tensor_aux[i][2] *= h_inv_dim_plus_one; + p->gradients.velocity_tensor_aux_norm[i][0] *= h_inv_dim_plus_one; + p->gradients.velocity_tensor_aux_norm[i][1] *= h_inv_dim_plus_one; + p->gradients.velocity_tensor_aux_norm[i][2] *= h_inv_dim_plus_one; + aux_norm[i][0] = p->gradients.velocity_tensor_aux_norm[i][0]; + aux_norm[i][1] = p->gradients.velocity_tensor_aux_norm[i][1]; + aux_norm[i][2] = p->gradients.velocity_tensor_aux_norm[i][2]; + } + + /* Invert the aux_norm matrix */ + gsl_matrix_view A_view = gsl_matrix_view_array((double *)aux_norm, 3, 3); + gsl_matrix *A = &A_view.matrix; + + double cond = condition_number(A); + if (cond < const_condition_number_upper_limit) { + gsl_matrix *A_inv = gsl_matrix_alloc(3, 3); + gsl_permutation *p_perm = gsl_permutation_alloc(3); + int signum; + + gsl_linalg_LU_decomp(A, p_perm, &signum); + gsl_linalg_LU_invert(A, p_perm, A_inv); + + for (int i = 0; i < 3; i++) { + p->gradients.velocity_tensor_aux_norm[i][0] = gsl_matrix_get(A_inv, i, 0); + p->gradients.velocity_tensor_aux_norm[i][1] = gsl_matrix_get(A_inv, i, 1); + p->gradients.velocity_tensor_aux_norm[i][2] = gsl_matrix_get(A_inv, i, 2); + } + + gsl_matrix_free(A_inv); + gsl_permutation_free(p_perm); + + /* aux_norm is equation 20 in Rosswog 2020, finalize the gradient in 19 */ + hydro_real_t aux_matrix[3][3] = {0}; + for (int j = 0; j < 3; j++) { + for (int i = 0; i < 3; i++) { + /** + * The indices of aux_norm and velocity_gradient_aux are important. + * aux_norm j: dx vector direction. + * k: kernel gradient direction + * + * velocity_gradient_aux i: dv vector direction + * k: kernel gradient direction + * + * aux_matrix j: spatial derivative direction + * i: velocity direction + */ + for (int k = 0; k < 3; k++) { + aux_matrix[j][i] += p->gradients.velocity_tensor_aux_norm[j][k] * + p->gradients.velocity_tensor_aux[i][k]; + } + } + } + + /* Copy over the matrices for use later */ + + /* D. Rennehan: For some reason, memcpy does not work here? Could it + * be because of the union in the particle struct? */ + for (int j = 0; j < 3; j++) { + hydro_vec_slope_limiter(p->gradients.dv_min[j], p->gradients.dv_max[j], + p->gradients.kernel_size, + aux_matrix[j]); + + p->gradients.velocity_tensor_aux[j][0] = aux_matrix[j][0]; + p->gradients.velocity_tensor_aux[j][1] = aux_matrix[j][1]; + p->gradients.velocity_tensor_aux[j][2] = aux_matrix[j][2]; + } + } + else { + + p->gradients.D_well_conditioned = 0; +#ifdef MAGMA2_DEBUG_CHECKS + p->debug.D_ill_conditioned_count++; +#endif + + /* Ensure no crazy gradients later */ + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { +#ifdef MAGMA2_DEBUG_CHECKS + p->debug.velocity_tensor_aux[i][j] = + p->gradients.velocity_tensor_aux[i][j]; + p->debug.velocity_tensor_aux_norm[i][j] = + p->gradients.velocity_tensor_aux_norm[i][j]; +#endif + + p->gradients.velocity_tensor_aux[i][j] = 0.; + p->gradients.velocity_tensor_aux_norm[i][j] = 0.; + } + } + } +} + +/** + * @brief Prepare a particle for the gradient calculation. + * + * This function is 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. + * @param hydro_props Hydrodynamic properties. + * @param pressure_floor The properties of the pressure floor. + */ +__attribute__((always_inline)) INLINE static void hydro_prepare_gradient( + struct part *restrict p, struct xpart *restrict xp, + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { + +#ifdef hydro_props_use_adiabatic_correction + /* Prepare the denominator for the adiabatic correction term */ + p->gradients.adiabatic_f_denominator = 0.; +#endif + + /* Compute the sound speed */ + const float pressure = hydro_get_comoving_pressure(p); + const float pressure_including_floor = + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); + const float soundspeed = + gas_soundspeed_from_pressure(p->rho, pressure_including_floor); + + /* Update variables. */ + p->force.pressure = pressure_including_floor; + p->force.soundspeed = soundspeed; + + + /* ------ Compute the kernel correction for SPH gradients ------ */ + + + /* Compute the "grad h" term */ + const float common_factor = p->h * hydro_dimension_inv / p->density.wcount; + float grad_W_term = -1.f; + + /* Ignore changing-kernel effects when h ~= h_max */ + if (p->h > 0.9999f * hydro_props->h_max) { + warning("h ~ h_max for particle with ID %lld (h: %g)", p->id, p->h); + } + else { + grad_W_term = common_factor * p->density.wcount_dh; + + if (grad_W_term < -0.9999f) { + /* if we get here, we either had very small neighbour contributions + (which should be treated as a no neighbour case in the ghost) or + a very weird particle distribution (e.g. particles sitting on + top of each other). Either way, we cannot use the normal + expression, since that would lead to overflow or excessive round + off and cause excessively high accelerations in the force loop */ + grad_W_term = -1.f; + warning( + "grad_W_term very small for particle with ID %lld (h: %g, wcount: " + "%g, wcount_dh: %g)", + p->id, p->h, p->density.wcount, p->density.wcount_dh); + } + } + + /* Update variables. */ + p->force.f = (grad_W_term > -1.f) ? 1.f / (1.f + grad_W_term) : 1.f; + +} + +/** + * @brief Resets the variables that are required for a gradient calculation. + * + * This function is called after hydro_prepare_gradient. + * + * @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_reset_gradient( + struct part *restrict p) { + + p->gradients.C_well_conditioned = 1; + + for (int i = 0; i < 3; i++) { + p->gradients.u[i] = 0.; + p->gradients.u_hessian[i][0] = 0.; + p->gradients.u_hessian[i][1] = 0.; + p->gradients.u_hessian[i][2] = 0.; + + p->gradients.correction_matrix[i][0] = 0.; + p->gradients.correction_matrix[i][1] = 0.; + p->gradients.correction_matrix[i][2] = 0.; + + p->gradients.velocity_tensor[i][0] = 0.; + p->gradients.velocity_tensor[i][1] = 0.; + p->gradients.velocity_tensor[i][2] = 0.; + + for (int j = 0; j < 3; j++) { + p->gradients.velocity_hessian[i][j][0] = 0.; + p->gradients.velocity_hessian[i][j][1] = 0.; + p->gradients.velocity_hessian[i][j][2] = 0.; + } + } +} + +/** + * @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) { + + /* Calculate smoothing length powers */ + 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 */ + + p->gradients.u[0] *= h_inv_dim; + p->gradients.u[1] *= h_inv_dim; + p->gradients.u[2] *= h_inv_dim; + + /* Temporary double for GSL */ + double correction_matrix[3][3] = {0}; + + /* Apply correct normalisation */ + for (int k = 0; k < 3; k++) { + p->gradients.correction_matrix[k][0] *= h_inv_dim; + p->gradients.correction_matrix[k][1] *= h_inv_dim; + p->gradients.correction_matrix[k][2] *= h_inv_dim; + + correction_matrix[k][0] = p->gradients.correction_matrix[k][0]; + correction_matrix[k][1] = p->gradients.correction_matrix[k][1]; + correction_matrix[k][2] = p->gradients.correction_matrix[k][2]; + + p->gradients.u_hessian[k][0] *= h_inv_dim; + p->gradients.u_hessian[k][1] *= h_inv_dim; + p->gradients.u_hessian[k][2] *= h_inv_dim; + + p->gradients.velocity_tensor[k][0] *= h_inv_dim; + p->gradients.velocity_tensor[k][1] *= h_inv_dim; + p->gradients.velocity_tensor[k][2] *= h_inv_dim; + + for (int j = 0; j < 3; j++) { + p->gradients.velocity_hessian[k][j][0] *= h_inv_dim; + p->gradients.velocity_hessian[k][j][1] *= h_inv_dim; + p->gradients.velocity_hessian[k][j][2] *= h_inv_dim; + } + } + + /* Invert the p->gradients.correction_matrix[3][3] matrix */ + gsl_matrix_view C_view = + gsl_matrix_view_array((double *)correction_matrix, 3, 3); + gsl_matrix *C = &C_view.matrix; + + const double cond = condition_number(C); + if (cond < const_condition_number_upper_limit) { + gsl_matrix *C_inv = gsl_matrix_alloc(3, 3); + gsl_permutation *p_perm = gsl_permutation_alloc(3); + int signum; + + gsl_linalg_LU_decomp(C, p_perm, &signum); + gsl_linalg_LU_invert(C, p_perm, C_inv); + + for (int k = 0; k < 3; k++) { + p->gradients.correction_matrix[k][0] = gsl_matrix_get(C_inv, k, 0); + p->gradients.correction_matrix[k][1] = gsl_matrix_get(C_inv, k, 1); + p->gradients.correction_matrix[k][2] = gsl_matrix_get(C_inv, k, 2); + } + + gsl_matrix_free(C_inv); + gsl_permutation_free(p_perm); + } + else { +#ifdef MAGMA2_DEBUG_CHECKS + for (int k = 0; k < 3; k++) { + p->debug.correction_matrix[k][0] = p->gradients.correction_matrix[k][0]; + p->debug.correction_matrix[k][1] = p->gradients.correction_matrix[k][1]; + p->debug.correction_matrix[k][2] = p->gradients.correction_matrix[k][2]; + } + + p->debug.C_ill_conditioned_count++; +#endif + + /* Ill-condition matrix, revert back to normal SPH gradients */ + p->gradients.C_well_conditioned = 0; + for (int k = 0; k < 3; k++) { + p->gradients.correction_matrix[k][0] = 0.; + p->gradients.correction_matrix[k][1] = 0.; + p->gradients.correction_matrix[k][2] = 0.; + } + } + + /* Contract the correction matrix with the internal energy gradient and with + * the velocity tensor */ + hydro_real_t u_gradient[3] = {0}; + hydro_real_t u_hessian[3][3] = {0}; + hydro_real_t velocity_tensor[3][3] = {0}; + hydro_real_t velocity_hessian[3][3][3] = {0}; + for (int j = 0; j < 3; j++) { + for (int i = 0; i < 3; i++) { + u_gradient[j] += p->gradients.correction_matrix[j][i] * + p->gradients.u[i]; + + for (int k = 0; k < 3; k++) { + u_hessian[j][i] += p->gradients.correction_matrix[j][k] * + p->gradients.u_hessian[i][k]; + velocity_tensor[j][i] += p->gradients.correction_matrix[j][k] * + p->gradients.velocity_tensor[i][k]; + + for (int m = 0; m < 3; m++) { + velocity_hessian[j][i][k] += p->gradients.correction_matrix[j][m] * + p->gradients.velocity_hessian[j][i][m]; + } + } + } + } + + /* Slope limiter for internal energy */ + hydro_vec_slope_limiter(p->gradients.du_min, p->gradients.du_max, + p->gradients.kernel_size, + u_gradient); + + /* Copy back over to the particle for later */ + p->gradients.u[0] = u_gradient[0]; + p->gradients.u[1] = u_gradient[1]; + p->gradients.u[2] = u_gradient[2]; + + for (int j = 0; j < 3; j++) { + p->gradients.u_hessian[j][0] = u_hessian[j][0]; + p->gradients.u_hessian[j][1] = u_hessian[j][1]; + p->gradients.u_hessian[j][2] = u_hessian[j][2]; + + hydro_vec_slope_limiter(p->gradients.dv_min[j], p->gradients.dv_max[j], + p->gradients.kernel_size, + velocity_tensor[j]); + + p->gradients.velocity_tensor[j][0] = velocity_tensor[j][0]; + p->gradients.velocity_tensor[j][1] = velocity_tensor[j][1]; + p->gradients.velocity_tensor[j][2] = velocity_tensor[j][2]; + + p->gradients.velocity_hessian[j][0][0] = velocity_hessian[j][0][0]; + p->gradients.velocity_hessian[j][0][1] = velocity_hessian[j][0][1]; + p->gradients.velocity_hessian[j][0][2] = velocity_hessian[j][0][2]; + + p->gradients.velocity_hessian[j][1][0] = velocity_hessian[j][1][0]; + p->gradients.velocity_hessian[j][1][1] = velocity_hessian[j][1][1]; + p->gradients.velocity_hessian[j][1][2] = velocity_hessian[j][1][2]; + + p->gradients.velocity_hessian[j][2][0] = velocity_hessian[j][2][0]; + p->gradients.velocity_hessian[j][2][1] = velocity_hessian[j][2][1]; + p->gradients.velocity_hessian[j][2][2] = velocity_hessian[j][2][2]; + } +} + +/** + * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. + * + * In the desperate case where a particle has no neighbours (likely because + * of the h_max ceiling), set the particle fields to something sensible to avoid + * NaNs in the next calculations. + * + * @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 */ + + warning( + "Gas particle with ID %lld treated as having no neighbours (h: %g, " + "wcount: %g).", + p->id, h, p->density.wcount); + + /* Re-set problematic values */ + p->rho = p->mass * kernel_root * h_inv_dim; + p->h_min = FLT_MAX; + p->dt_min = FLT_MAX; + p->density.wcount = kernel_root * h_inv_dim; + p->density.rho_dh = 0.f; + p->density.wcount_dh = 0.f; + +#ifdef MAGMA2_DEBUG_CHECKS + p->debug.num_ngb = 0; + p->debug.v_sig_visc_max = 0; + p->debug.v_sig_cond_max = 0; +#endif + p->gradients.C_well_conditioned = 0; + p->gradients.D_well_conditioned = 0; + p->gradients.u_well_conditioned = 0; + p->gradients.du_min = 0.; + p->gradients.du_max = 0.; + p->gradients.kernel_size = (hydro_real_t)h; + + for (int i = 0; i < 3; i++) { + p->gradients.u[i] = 0.; + p->gradients.u_aux[i] = 0.; + p->gradients.u_aux_norm[i] = 0.; + + p->gradients.dv_min[i] = 0.; + p->gradients.dv_max[i] = 0.; + + for (int j = 0; j < 3; j++) { + p->gradients.correction_matrix[i][j] = 0.; + p->gradients.velocity_tensor[i][j] = 0.; + p->gradients.velocity_tensor_aux[i][j] = 0.; + p->gradients.velocity_tensor_aux_norm[i][j] = 0.; + p->gradients.u_hessian[i][j] = 0.; + + for (int k = 0; k < 3; k++) { + p->gradients.velocity_hessian[i][j][k] = 0.; + } + } + } +} + +/** + * @brief Prepare a particle for the force calculation. + * + * This function is called in the ghost task to convert some quantities coming + * from the density loop over neighbours into quantities ready to be used in the + * force loop over neighbours. Quantities are typically read from the density + * sub-structure and written to the force sub-structure. + * Examples of calculations done here include the calculation of viscosity term + * constants, thermal conduction terms, hydro conversions, etc. + * + * @param p The particle to act upon + * @param xp The extended particle data to act upon + * @param cosmo The current cosmological model. + * @param hydro_props Hydrodynamic properties. + * @param pressure_floor The properties of the pressure floor. + * @param dt_alpha The time-step used to evolve non-cosmological quantities such + * as the artificial viscosity. + * @param dt_therm The time-step used to evolve hydrodynamical quantities. + */ +__attribute__((always_inline)) INLINE static void hydro_prepare_force( + struct part *restrict p, struct xpart *restrict xp, + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor, const float dt_alpha, + const float dt_therm) { + + /* First estimates for the timestepping. Missing the kernel_gamma factors + * for now, but will be added at the end of the force loop. */ + p->h_min = FLT_MAX; + p->dt_min = FLT_MAX; + p->gradients.balsara = hydro_get_balsara_limiter(p, cosmo); + +#ifdef hydro_props_use_adiabatic_correction + const hydro_real_t prev_f = p->force.f; + const hydro_real_t rho_inv = 1. / hydro_get_comoving_density(p); + /* Finish final kernel gradient correction factor */ + if (p->gradients.adiabatic_f_denominator != 0.) { + const hydro_real_t kernel_r2 = p->gradients.adiabatic_f_numerator; + const hydro_real_t weighted_kernel_r2_inv = + 1. / p->gradients.adiabatic_f_denominator; + p->force.f = rho_inv * kernel_r2 * weighted_kernel_r2_inv; + } + else { + p->force.f = 1.; + } + +#ifdef MAGMA2_DEBUG_CHECKS + if (p->force.f > 100.) { + warning("Final force factor is very high for particle with ID %lld" + " (prev f: %g, f: %g, numerator: %g, denominator: %g, rho_inv: %g)", + p->id, prev_f, p->force.f, p->gradients.adiabatic_f_numerator, + p->gradients.adiabatic_f_denominator, rho_inv); + } +#endif +#endif +} + +/** + * @brief Reset acceleration fields of a particle + * + * Resets all hydro acceleration and time derivative fields in preparation + * for the sums taking place in the various force tasks. + * + * @param p The particle to act upon + */ +__attribute__((always_inline)) INLINE static void hydro_reset_acceleration( + struct part *restrict 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->u_dt = 0.0f; + p->u_dt_cond = 0.0f; + 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. + * @param cosmo The cosmological model. + * @param pressure_floor The properties of the pressure floor. + */ +__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( + struct part *restrict p, const struct xpart *restrict xp, + const struct cosmology *cosmo, + const struct pressure_floor_props *pressure_floor) { + /* Re-set the predicted velocities */ + p->v[0] = xp->v_full[0]; + p->v[1] = xp->v_full[1]; + p->v[2] = xp->v_full[2]; + + /* Re-set the entropy */ + p->u = xp->u_full; + + /* Compute the sound speed */ + const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); + const float pressure_including_floor = + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); + const float soundspeed = + gas_soundspeed_from_pressure(p->rho, pressure_including_floor); + + p->force.pressure = pressure_including_floor; + p->force.soundspeed = soundspeed; + + /* Signal speed */ + const float v_sig = const_viscosity_alpha_prefactor * soundspeed; + + /* Update the signal velocity, if we need to. */ + p->dt_min = min(p->dt_min, p->h_min / v_sig); +} + +/** + * @brief Predict additional particle fields forward in time when drifting + * + * Additional hydrodynamic quantites are drifted forward in time here. These + * include thermal quantities (thermal energy or total energy or entropy, ...). + * + * Note the different time-step sizes used for the different quantities as they + * include cosmological factors. + * + * @param p The particle. + * @param xp The extended data of the particle. + * @param dt_drift The drift time-step for positions. + * @param dt_therm The drift time-step for thermal quantities. + * @param dt_kick_grav The time-step for gravity quantities. + * @param cosmo The cosmological model. + * @param hydro_props The properties of the hydro scheme. + * @param floor_props The properties of the entropy floor. + * @param pressure_floor The properties of the pressure floor. + */ +__attribute__((always_inline)) INLINE static void hydro_predict_extra( + struct part *restrict p, const struct xpart *restrict xp, float dt_drift, + float dt_therm, float dt_kick_grav, const struct cosmology *cosmo, + const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props, + const struct pressure_floor_props *pressure_floor) { + + /* Predict the internal energy */ + p->u += p->u_dt * dt_therm; + + const float h_inv = 1.f / p->h; + + /* Predict smoothing length */ + const float w1 = p->force.h_dt * h_inv * dt_drift; + if (fabsf(w1) < 0.2f) + p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */ + else + p->h *= expf(w1); + + /* Predict density and weighted pressure */ + const float w2 = -hydro_dimension * w1; + if (fabsf(w2) < 0.2f) { + const float expf_approx = + approx_expf(w2); /* 4th order expansion of exp(w) */ + p->rho *= expf_approx; + } else { + const float expf_exact = expf(w2); + p->rho *= expf_exact; + } + + /* Check against entropy floor - explicitly do this after drifting the + * density as this has a density dependence. */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A); + + /* Check against absolute minimum */ + const float min_u = + hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy; + + p->u = max(p->u, floor_u); + p->u = max(p->u, min_u); + + /* Compute the new sound speed */ + const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); + const float pressure_including_floor = + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); + const float soundspeed = + gas_soundspeed_from_pressure(p->rho, pressure_including_floor); + + p->force.pressure = pressure_including_floor; + p->force.soundspeed = soundspeed; + + /* Signal speed */ + const float v_sig = const_viscosity_alpha_prefactor * soundspeed; + + /* Update signal velocity if we need to */ + p->dt_min = min(p->dt_min, p->h_min / v_sig); +} + +/** + * @brief Returns the sum term for the dh/dt calculation + * + * + * @param m The particle mass of the neighbour + * @param rho_inv The inverse density of the neighbour + * @param r_inv The inverse distance between particles + * @param w_dr The kernel gradient for this particle + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_get_h_dt_sum(const hydro_real_t dv_dot_dx, + const hydro_real_t dv_dot_G, + const hydro_real_t m, + const hydro_real_t rho_inv, + const hydro_real_t r_inv, + const hydro_real_t w_dr) { + + hydro_real_t dvdx = 0.; +#ifdef hydro_props_dh_dt_estimator_type +#if (hydro_props_dh_dt_estimator_type == 0) + const hydro_real_t grad = r_inv * w_dr; + const hydro_real_t wt = m * rho_inv; + dvdx = dv_dot_dx; +#elif (hydro_props_dh_dt_estimator_type == 1) + const hydro_real_t grad = r_inv * w_dr; + const hydro_real_t wt = 1.; + dvdx = dv_dot_dx; +#elif (hydro_props_dh_dt_estimator_type == 2) + const hydro_real_t grad = 1.; + const hydro_real_t wt = m * rho_inv; + dvdx = dv_dot_G; +#else + error("Compiled with an unknown dh/dt estimator type"); +#endif +#else + error("Must compile with hydro_props_dh_dt_estimator_type."); +#endif + + return wt * dvdx * grad; +} + +/** + * @brief Returns the normalization for the dh/dt sum + * + * + * @param p The particle + */ +__attribute__((always_inline)) INLINE static +hydro_real_t hydro_get_h_dt_norm(struct part *restrict p) { + + hydro_real_t renormalization = 1.; + +#ifdef hydro_props_dh_dt_estimator_type +#if (hydro_props_dh_dt_estimator_type == 1) + renormalization = p->force.f / p->gradients.wcount; +#endif +#else + error("Must compile with hydro_props_dh_dt_estimator_type."); +#endif + + return renormalization * p->h * hydro_dimension_inv; +} + +/** + * @brief Finishes the force calculation. + * + * Multiplies the force and accelerations by the appropiate constants + * and add the self-contribution term. In most cases, there is little + * to do here. + * + * Cosmological terms are also added/multiplied here. + * + * @param p The particle to act upon + * @param cosmo The current cosmological model. + */ +__attribute__((always_inline)) INLINE static void hydro_end_force( + struct part *restrict p, const struct cosmology *cosmo) { + + p->force.h_dt *= hydro_get_h_dt_norm(p); +} + +/** + * @brief Kick the additional variables + * + * Additional hydrodynamic quantites are kicked forward in time here. These + * include thermal quantities (thermal energy or total energy or entropy, ...). + * + * @param p The particle to act upon. + * @param xp The particle extended data to act upon. + * @param dt_therm The time-step for this kick (for thermodynamic quantities). + * @param dt_grav The time-step for this kick (for gravity quantities). + * @param dt_grav_mesh The time-step for this kick (mesh gravity). + * @param dt_hydro The time-step for this kick (for hydro quantities). + * @param dt_kick_corr The time-step for this kick (for gravity corrections). + * @param cosmo The cosmological model. + * @param hydro_props The constants used in the scheme + * @param floor_props The properties of the entropy floor. + */ +__attribute__((always_inline)) INLINE static void hydro_kick_extra( + struct part *restrict p, struct xpart *restrict xp, float dt_therm, + float dt_grav, float dt_grav_mesh, float dt_hydro, float dt_kick_corr, + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + /* Integrate the internal energy forward in time */ + const float delta_u = p->u_dt * dt_therm; + + /* Do not decrease the energy by more than a factor of 2*/ + xp->u_full = max(xp->u_full + delta_u, 0.5f * xp->u_full); + + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + const float floor_u = gas_internal_energy_from_entropy(p->rho, floor_A); + + /* Check against absolute minimum */ + const float min_u = + hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy; + + /* Take highest of both limits */ + const float energy_min = max(min_u, floor_u); + + if (xp->u_full < energy_min) { + xp->u_full = energy_min; + p->u_dt = 0.f; + } +} + +/** + * @brief Converts hydro quantity of a particle at the start of a run + * + * This function is called once at the end of the engine_init_particle() + * routine (at the start of a calculation) after the densities of + * particles have been computed. + * This can be used to convert internal energy into entropy for instance. + * + * @param p The particle to act upon + * @param xp The extended particle to act upon + * @param cosmo The cosmological model. + * @param hydro_props The constants used in the scheme. + * @param pressure_floor The properties of the pressure floor. + */ +__attribute__((always_inline)) INLINE static void hydro_convert_quantities( + struct part *restrict p, struct xpart *restrict xp, + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct pressure_floor_props *pressure_floor) { + /* Convert the physcial internal energy to the comoving one. */ + /* u' = a^(3(g-1)) u */ + const float factor = 1.f / cosmo->a_factor_internal_energy; + p->u *= factor; + xp->u_full = p->u; + + /* Apply the minimal energy limit */ + const float min_comoving_energy = + hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy; + if (xp->u_full < min_comoving_energy) { + xp->u_full = min_comoving_energy; + p->u = min_comoving_energy; + p->u_dt = 0.f; + } + + /* Set the initial value of the artificial viscosity based on the non-variable + schemes for safety */ + + const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); + const float pressure_including_floor = + pressure_floor_get_comoving_pressure(p, pressure_floor, pressure, cosmo); + const float soundspeed = + gas_soundspeed_from_pressure(p->rho, pressure_including_floor); + + p->force.pressure = pressure_including_floor; + p->force.soundspeed = soundspeed; +} + +/** + * @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 or assignments between the particle + * and extended particle fields. + * + * @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 *restrict p, struct xpart *restrict xp) { + p->time_bin = 0; + + xp->v_full[0] = p->v[0]; + xp->v_full[1] = p->v[1]; + xp->v_full[2] = p->v[2]; + xp->u_full = p->u; + + p->gradients.C_well_conditioned = 1; + p->gradients.D_well_conditioned = 1; + p->gradients.u_well_conditioned = 1; +#ifdef MAGMA2_DEBUG_CHECKS + p->debug.N_force_low_order_grad = 0; + p->debug.N_force_high_order_grad = 0; + p->debug.C_ill_conditioned_count = 0; + p->debug.D_ill_conditioned_count = 0; + p->debug.u_ill_conditioned_count = 0; + for (int i = 0; i < 3; i++) { + p->debug.correction_matrix[i][0] = 0.; + p->debug.correction_matrix[i][1] = 0.; + p->debug.correction_matrix[i][2] = 0.; + + p->debug.velocity_tensor_aux[i][0] = 0.; + p->debug.velocity_tensor_aux[i][1] = 0.; + p->debug.velocity_tensor_aux[i][2] = 0.; + + p->debug.velocity_tensor_aux_norm[i][0] = 0.; + p->debug.velocity_tensor_aux_norm[i][1] = 0.; + p->debug.velocity_tensor_aux_norm[i][2] = 0.; + } +#endif + + hydro_reset_acceleration(p); + hydro_init_part(p, NULL); +} + +/** + * @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->u = u_init; +} + +#endif /* SWIFT_MAGMA2_HYDRO_H */ diff --git a/src/hydro/MAGMA2/hydro_debug.h b/src/hydro/MAGMA2/hydro_debug.h new file mode 100644 index 0000000000..84852bbb73 --- /dev/null +++ b/src/hydro/MAGMA2/hydro_debug.h @@ -0,0 +1,66 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) & + * Matthieu Schaller (schaller@strw.leidenuniv.nl) + * 2025 Doug Rennehan (douglas.rennehan@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 . + * + ******************************************************************************/ + +#ifndef SWIFT_MAGMA2_HYDRO_DEBUG_H +#define SWIFT_MAGMA2_HYDRO_DEBUG_H + +/** + * @file MAGMA2/hydro_debug.h + * @brief Density-Energy non-conservative implementation of SPH, + * with added MAGMA2 physics (Rosswog 2020) (Debugging routines) + */ + +__attribute__((always_inline)) INLINE static void hydro_debug_particle( + const struct part* p, const struct xpart* xp) { + warning("[PID%lld] part:", p->id); + warning("[PID%lld] x=[%.3e,%.3e,%.3e]", p->id, p->x[0], p->x[1], p->x[2]); + warning("[PID%lld] v=[%.3e,%.3e,%.3e]", p->id, p->v[0], p->v[1], p->v[2]); + warning("[PID%lld] a=[%.3e,%.3e,%.3e]", p->id, p->a_hydro[0], p->a_hydro[1], + p->a_hydro[2]); + warning("[PID%lld] u=%.3e, du/dt=%.3e P=%.3e", p->id, p->u, + p->u_dt, hydro_get_comoving_pressure(p)); + warning("[PID%lld] h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e", p->id, + p->h, p->force.h_dt, (int)p->density.wcount, p->mass, + p->density.rho_dh); + warning( + "[PID%lld] time_bin=%d, rho=%.3e, velocity_gradient=[" + "[%.3e,%.3e,%.3e]," + "[%.3e,%.3e,%.3e]," + "[%.3e,%.3e,%.3e]]", + p->id, p->time_bin, p->rho, + p->gradients.velocity_tensor[0][0], + p->gradients.velocity_tensor[0][1], + p->gradients.velocity_tensor[0][2], + p->gradients.velocity_tensor[1][0], + p->gradients.velocity_tensor[1][1], + p->gradients.velocity_tensor[1][2], + p->gradients.velocity_tensor[2][0], + p->gradients.velocity_tensor[2][1], + p->gradients.velocity_tensor[2][2]); + + if (xp != NULL) { + warning("[PID%lld] xpart:", p->id); + warning("[PID%lld] v_full=[%.3e,%.3e,%.3e]", p->id, xp->v_full[0], + xp->v_full[1], xp->v_full[2]); + } +} + +#endif /* SWIFT_MAGMA2_HYDRO_DEBUG_H */ diff --git a/src/hydro/MAGMA2/hydro_iact.h b/src/hydro/MAGMA2/hydro_iact.h new file mode 100644 index 0000000000..a6255d7577 --- /dev/null +++ b/src/hydro/MAGMA2/hydro_iact.h @@ -0,0 +1,1501 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) & + * Matthieu Schaller (schaller@strw.leidenuniv.nl) + * 2025 Doug Rennehan (douglas.rennehan@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 . + * + ******************************************************************************/ +#ifndef SWIFT_MAGMA2_HYDRO_IACT_H +#define SWIFT_MAGMA2_HYDRO_IACT_H + +/** + * @file MAGMA2/hydro_iact.h + * @brief Density-Energy non-conservative implementation of SPH, + * with added MAGMA2 physics (Rosswog 2020) (interaction routines) + */ + +#include "adaptive_softening_iact.h" +#include "adiabatic_index.h" +#include "hydro_parameters.h" +#include "minmax.h" +#include "signal_velocity.h" + +/** + * @brief Density interaction between two particles. + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing-length of part*icle i. + * @param hj Comoving smoothing-length of part*icle j. + * @param pi First part*icle. + * @param pj Second part*icle. + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void runner_iact_density( + const float r2, const float dx[3], const float hi, const float hj, + struct part* restrict pi, struct part* restrict pj, const float a, + const float H) { + + /* Kernel weights to be filled */ + float wi, wj, wi_dx, wj_dx; + + const hydro_real_t r = sqrt(r2); + + /* Get the masses. */ + const hydro_real_t mi = pi->mass; + const hydro_real_t mj = pj->mass; + + /* Compute density of pi. */ + const hydro_real_t hi_inv = 1. / hi; + const float xi = r * hi_inv; + + kernel_deval(xi, &wi, &wi_dx); + + pi->rho += mj * wi; + pi->density.rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx); + + pi->density.wcount += wi; + pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx); + + adaptive_softening_add_correction_term(pi, xi, hi_inv, mj); + + /* Compute density of pj. */ + const hydro_real_t hj_inv = 1. / hj; + const float xj = r * hj_inv; + kernel_deval(xj, &wj, &wj_dx); + + pj->rho += mi * wj; + pj->density.rho_dh -= mi * (hydro_dimension * wj + xj * wj_dx); + + pj->density.wcount += wj; + pj->density.wcount_dh -= (hydro_dimension * wj + xj * wj_dx); + + adaptive_softening_add_correction_term(pj, xj, hj_inv, mi); + + /* For slope limiter */ + pi->gradients.kernel_size = fmax(r, pi->gradients.kernel_size); + pj->gradients.kernel_size = fmax(r, pj->gradients.kernel_size); + + /* Now we need to compute the derivative terms */ + const hydro_real_t r_inv = r ? 1.0 / r : 0.0; + const hydro_real_t faci = mj * wi_dx * r_inv; + const hydro_real_t facj = mi * wj_dx * r_inv; + + /* Equations 19 & 20 in Rosswog 2020. Compute the internal energy auxiliary + * vector and norm for the gradient */ + const hydro_real_t du = pi->u - pj->u; + + /* For slope limiter */ + pi->gradients.du_min = fmin(-du, pi->gradients.du_min); + pi->gradients.du_max = fmax(-du, pi->gradients.du_max); + + pj->gradients.du_min = fmin(du, pj->gradients.du_min); + pj->gradients.du_max = fmax(du, pj->gradients.du_max); + + pi->gradients.u_aux[0] += du * dx[0] * faci; + pi->gradients.u_aux[1] += du * dx[1] * faci; + pi->gradients.u_aux[2] += du * dx[2] * faci; + + pj->gradients.u_aux[0] += du * dx[0] * facj; + pj->gradients.u_aux[1] += du * dx[1] * facj; + pj->gradients.u_aux[2] += du * dx[2] * facj; + + pi->gradients.u_aux_norm[0] += dx[0] * dx[0] * faci; + pi->gradients.u_aux_norm[1] += dx[1] * dx[1] * faci; + pi->gradients.u_aux_norm[2] += dx[2] * dx[2] * faci; + + pj->gradients.u_aux_norm[0] += dx[0] * dx[0] * facj; + pj->gradients.u_aux_norm[1] += dx[1] * dx[1] * facj; + pj->gradients.u_aux_norm[2] += dx[2] * dx[2] * facj; + + const hydro_real_t dv[3] = {pi->v[0] - pj->v[0], + pi->v[1] - pj->v[1], + pi->v[2] - pj->v[2]}; + + /* Equations 19 & 20 in Rosswog 2020. Signs are all positive because + * dv * dx always results in a positive sign. */ + for (int i = 0; i < 3; i++) { + + /* For slope limiter */ + pi->gradients.dv_min[i] = fmin(-dv[i], pi->gradients.dv_min[i]); + pi->gradients.dv_max[i] = fmax(-dv[i], pi->gradients.dv_max[i]); + + pj->gradients.dv_min[i] = fmin(dv[i], pj->gradients.dv_min[i]); + pj->gradients.dv_max[i] = fmax(dv[i], pj->gradients.dv_max[i]); + + for (int k = 0; k < 3; k++) { + pi->gradients.velocity_tensor_aux[i][k] += dv[i] * dx[k] * faci; + pj->gradients.velocity_tensor_aux[i][k] += dv[i] * dx[k] * facj; + + pi->gradients.velocity_tensor_aux_norm[i][k] += dx[i] * dx[k] * faci; + pj->gradients.velocity_tensor_aux_norm[i][k] += dx[i] * dx[k] * facj; + } + } + +#ifdef MAGMA2_DEBUG_CHECKS + /* Number of neighbors */ + pi->debug.num_ngb++; + pj->debug.num_ngb++; +#endif + +#ifdef hydro_props_use_adiabatic_correction + /* Needed for the adiabatic kernel correction factor */ + pi->gradients.adiabatic_f_numerator += mj * r2 * wi; + pj->gradients.adiabatic_f_numerator += mi * r2 * wj; +#endif + +} + +/** + * @brief Density interaction between two particles (non-symmetric). + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing-length of part*icle i. + * @param hj Comoving smoothing-length of part*icle j. + * @param pi First part*icle. + * @param pj Second part*icle (not updated). + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( + const float r2, const float dx[3], const float hi, const float hj, + struct part* restrict pi, const struct part* restrict pj, const float a, + const float H) { + + /* Kernel weights to be filled */ + float wi, wi_dx; + + /* Get the masses. */ + const hydro_real_t mj = pj->mass; + + /* Get r and r inverse. */ + const hydro_real_t r = sqrt(r2); + + const hydro_real_t h_inv = 1. / hi; + const float xi = r * h_inv; + kernel_deval(xi, &wi, &wi_dx); + + pi->rho += mj * wi; + pi->density.rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx); + + pi->density.wcount += wi; + pi->density.wcount_dh -= (hydro_dimension * wi + xi * wi_dx); + + adaptive_softening_add_correction_term(pi, xi, h_inv, mj); + + /* For slope limiter */ + pi->gradients.kernel_size = fmax(r, pi->gradients.kernel_size); + + const hydro_real_t r_inv = r ? 1.0 / r : 0.0; + const hydro_real_t faci = mj * wi_dx * r_inv; + + /* Equations 19 & 20 in Rosswog 2020. Compute the internal energy auxiliary + * vector and norm for the gradient */ + const hydro_real_t du = pi->u - pj->u; + + /* For slope limiter */ + pi->gradients.du_min = fmin(-du, pi->gradients.du_min); + pi->gradients.du_max = fmax(-du, pi->gradients.du_max); + + pi->gradients.u_aux[0] += du * dx[0] * faci; + pi->gradients.u_aux[1] += du * dx[1] * faci; + pi->gradients.u_aux[2] += du * dx[2] * faci; + + pi->gradients.u_aux_norm[0] += dx[0] * dx[0] * faci; + pi->gradients.u_aux_norm[1] += dx[1] * dx[1] * faci; + pi->gradients.u_aux_norm[2] += dx[2] * dx[2] * faci; + + const hydro_real_t dv[3] = {pi->v[0] - pj->v[0], + pi->v[1] - pj->v[1], + pi->v[2] - pj->v[2]}; + + /* Equations 19 & 20 in Rosswog 2020. Signs are all positive because + * dv * dx always results in a positive sign. */ + for (int i = 0; i < 3; i++) { + + /* For slope limiter */ + pi->gradients.dv_min[i] = fmin(-dv[i], pi->gradients.dv_min[i]); + pi->gradients.dv_max[i] = fmax(-dv[i], pi->gradients.dv_max[i]); + + for (int k = 0; k < 3; k++) { + pi->gradients.velocity_tensor_aux[i][k] += dv[i] * dx[k] * faci; + pi->gradients.velocity_tensor_aux_norm[i][k] += dx[i] * dx[k] * faci; + } + } + +#ifdef MAGMA2_DEBUG_CHECKS + /* Neighbour number */ + pi->debug.num_ngb++; +#endif + +#ifdef hydro_props_use_adiabatic_correction + /* Needed for the adiabatic kernel correction factor */ + pi->gradients.adiabatic_f_numerator += mj * r2 * wi; +#endif + +} + +/** + * @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( + const float r2, const float dx[3], const float hi, const float hj, + struct part* restrict pi, struct part* restrict pj, const float a, + const float H) { + + /* Get particle properties */ + const hydro_real_t mi = hydro_get_mass(pi); + const hydro_real_t mj = hydro_get_mass(pj); + + const hydro_real_t rhoi = hydro_get_comoving_density(pi); + const hydro_real_t rhoj = hydro_get_comoving_density(pj); + + const hydro_real_t rhoi_inv = 1. / rhoi; + const hydro_real_t rhoj_inv = 1. / rhoj; + + const hydro_real_t r = sqrt(r2); + + float wi, wi_dx, wj, wj_dx; + + const float xi = r / hi; + const float xj = r / hj; + + kernel_deval(xi, &wi, &wi_dx); + kernel_deval(xj, &wj, &wj_dx); + + const hydro_real_t faci = mj * rhoj_inv * wi; + const hydro_real_t facj = mi * rhoi_inv * wj; + + /* Compute all of the first-order gradients, and second-order gradients */ + + /* Rosswog 2020 Equation 18 gradients. In the paper he uses (vj - vi) and + * (rj - ri), however this is symmetric so no sign problems. */ + + /* Internal energy gradient */ + const hydro_real_t du = pi->u - pj->u; + + pi->gradients.u[0] += du * dx[0] * faci; + pi->gradients.u[1] += du * dx[1] * faci; + pi->gradients.u[2] += du * dx[2] * faci; + + pj->gradients.u[0] += du * dx[0] * facj; + pj->gradients.u[1] += du * dx[1] * facj; + pj->gradients.u[2] += du * dx[2] * facj; + + /* Velocity gradients */ + const hydro_real_t dv[3] = {pi->v[0] - pj->v[0], + pi->v[1] - pj->v[1], + pi->v[2] - pj->v[2]}; + + + for (int k = 0; k < 3; k++) { + const hydro_real_t du_k = pi->gradients.u_aux[k] - pj->gradients.u_aux[k]; + + for (int i = 0; i < 3; i++) { + pi->gradients.u_hessian[k][i] += du_k * dx[i] * faci; + pj->gradients.u_hessian[k][i] += du_k * dx[i] * facj; + + /* dx is signed as (pi - pj), but it is symmetric so we add */ + pi->gradients.correction_matrix[k][i] += dx[k] * dx[i] * faci; + pj->gradients.correction_matrix[k][i] += dx[k] * dx[i] * facj; + + /* Indices in Rosswog 2020 are i for dv and k for dx. In this loop, + * they are swapped just because correction_matrix is computed with + * the paper indices. */ + pi->gradients.velocity_tensor[k][i] += dv[k] * dx[i] * faci; + pj->gradients.velocity_tensor[k][i] += dv[k] * dx[i] * facj; + + const hydro_real_t dv_grad_ki = pi->gradients.velocity_tensor_aux[k][i] - + pj->gradients.velocity_tensor_aux[k][i]; + + /* Equation 19 indices: + * Index i: velocity direction + * Index k: gradient direction + * + * Our indices: + * Index k: velocity direction + * Index i: gradient direction + * Index j: second derivative gradient direction + */ + for (int j = 0; j < 3; j++) { + pi->gradients.velocity_hessian[k][i][j] += dv_grad_ki * dx[j] * faci; + pj->gradients.velocity_hessian[k][i][j] += dv_grad_ki * dx[j] * facj; + } + } + } + +#ifdef hydro_props_use_adiabatic_correction + /* Correction terms for div v */ + pi->gradients.adiabatic_f_denominator += mj * rhoj_inv * r2 * wi; + pj->gradients.adiabatic_f_denominator += mi * rhoi_inv * r2 * wj; +#endif + +} + +/** + * @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( + const float r2, const float dx[3], const float hi, const float hj, + struct part* restrict pi, struct part* restrict pj, const float a, + const float H) { + + /* Get particle properties */ + const hydro_real_t mj = hydro_get_mass(pj); + const hydro_real_t rhoj = hydro_get_comoving_density(pj); + const hydro_real_t rhoj_inv = 1. / rhoj; + const hydro_real_t r = sqrt(r2); + float wi, wi_dx; + const float xi = r / hi; + kernel_deval(xi, &wi, &wi_dx); + const hydro_real_t faci = mj * rhoj_inv * wi; + + /* Compute all of the first-order gradients, and second-order gradients */ + + /* Rosswog 2020 Equation 18 gradients. In the paper he uses (vj - vi) and + * (rj - ri), however this is symmetric so no sign problems. */ + + /* Internal energy gradient */ + const hydro_real_t du = pi->u - pj->u; + + pi->gradients.u[0] += du * dx[0] * faci; + pi->gradients.u[1] += du * dx[1] * faci; + pi->gradients.u[2] += du * dx[2] * faci; + + /* Velocity gradients */ + const hydro_real_t dv[3] = {pi->v[0] - pj->v[0], + pi->v[1] - pj->v[1], + pi->v[2] - pj->v[2]}; + + + for (int k = 0; k < 3; k++) { + const hydro_real_t du_k = pi->gradients.u_aux[k] - pj->gradients.u_aux[k]; + + for (int i = 0; i < 3; i++) { + pi->gradients.u_hessian[k][i] += du_k * dx[i] * faci; + + /* dx is signed as (pi - pj), but it is symmetric so we add */ + pi->gradients.correction_matrix[k][i] += dx[k] * dx[i] * faci; + + /* Indices in Rosswog 2020 are i for dv and k for dx. In this loop, + * they are swapped just because correction_matrix is computed with + * the paper indices. */ + pi->gradients.velocity_tensor[k][i] += dv[k] * dx[i] * faci; + + const hydro_real_t dv_grad_ki = pi->gradients.velocity_tensor_aux[k][i] - + pj->gradients.velocity_tensor_aux[k][i]; + + /* Equation 19 indices: + * Index i: velocity direction + * Index k: gradient direction + * + * Our indices: + * Index k: velocity direction + * Index i: gradient direction + * Index j: second derivative gradient direction + */ + for (int j = 0; j < 3; j++) { + pi->gradients.velocity_hessian[k][i][j] += dv_grad_ki * dx[j] * faci; + } + } + } + +#ifdef hydro_props_use_adiabatic_correction + /* Correction terms for div v */ + pi->gradients.adiabatic_f_denominator += mj * rhoj_inv * r2 * wi; +#endif + +} + +/** + * @brief Force interaction between two particles. + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing-length of part*icle i. + * @param hj Comoving smoothing-length of part*icle j. + * @param pi First part*icle. + * @param pj Second part*icle. + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void runner_iact_force( + const float r2, const float dx[3], const float hi, const float hj, + struct part* restrict pi, struct part* restrict pj, const float a, + const float H) { + + /* Cosmological factors entering the EoMs */ + const hydro_real_t fac_mu = pow_three_gamma_minus_five_over_two(a); + const hydro_real_t a2_Hubble = a * a * H; + + const hydro_real_t r = sqrt(r2); + const hydro_real_t r_inv = r ? 1.0 / r : 0.0; + + /* Recover some data */ + const hydro_real_t mj = pj->mass; + const hydro_real_t mi = pi->mass; + + const hydro_real_t rhoi = pi->rho; + const hydro_real_t rhoj = pj->rho; + const hydro_real_t rhoi_inv = 1. / rhoi; + const hydro_real_t rhoj_inv = 1. / rhoj; + const hydro_real_t rhoij_inv = rhoi_inv * rhoj_inv; + + const hydro_real_t pressurei = pi->force.pressure; + const hydro_real_t pressurej = pj->force.pressure; + + /* Get the kernel for hi. */ + const hydro_real_t hi_inv = 1. / hi; + const hydro_real_t hi_inv_dim = pow_dimension(hi_inv); + const float xi = r * hi_inv; + float wi, wi_dx; + kernel_deval(xi, &wi, &wi_dx); + + /* Get the kernel for hj. */ + const hydro_real_t hj_inv = 1. / hj; + const hydro_real_t hj_inv_dim = pow_dimension(hj_inv); + const float xj = r * hj_inv; + float wj, wj_dx; + kernel_deval(xj, &wj, &wj_dx); + + /* For dh/dt and fall-back SPH */ + const hydro_real_t hi_inv_dim_plus_one = hi_inv * hi_inv_dim; /* 1/h^(d+1) */ + const hydro_real_t hj_inv_dim_plus_one = hj_inv * hj_inv_dim; + const hydro_real_t wi_dr = hi_inv_dim_plus_one * wi_dx; + const hydro_real_t wj_dr = hj_inv_dim_plus_one * wj_dx; + + /* Peculiar velocity difference vector */ + const hydro_real_t dv[3] = {pi->v[0] - pj->v[0], + pi->v[1] - pj->v[1], + pi->v[2] - pj->v[2]}; + const hydro_real_t dv_raw[3] = {dv[0], dv[1], dv[2]}; + + /* For MAGMA2, this is the full anti-symmetric gradient vector. For the + * fallback Gasoline2-style SPH, this will just be the direction vector + * between the two particles (r_i - r_j). */ + hydro_real_t G_ij[3] = {0., 0., 0.}; + hydro_real_t G_rad_ij[3] = {0., 0., 0.}; + hydro_real_t G_ij_norm = 0.; + hydro_real_t G_rad_ij_norm = 0.; + + /* Separation vectors and swapped */ + const hydro_real_t dx_ij[3] = {dx[0], dx[1], dx[2]}; + const hydro_real_t dx_ji[3] = {-dx[0], -dx[1], -dx[2]}; + hydro_real_t dx_ij_hat[3] = {0., 0., 0.}; + hydro_vec3_unit(dx_ij, dx_ij_hat); + + /* These are set whether or not we fall back onto SPH gradients */ + hydro_real_t visc_acc_term = 0.; + hydro_real_t sph_acc_term = (pressurei + pressurej) * rhoij_inv; + hydro_real_t sph_du_term_i = pressurei * rhoij_inv; + hydro_real_t sph_du_term_j = pressurej * rhoij_inv; + hydro_real_t visc_du_term = 0.; + hydro_real_t cond_du_term = 0.; + + /* Correction factor must be well-conditioned. */ + const unsigned char C_well_conditioned = + (pi->gradients.C_well_conditioned && pj->gradients.C_well_conditioned); + /* First order density-independent velocity gradients */ + const unsigned char D_well_conditioned = + (pi->gradients.D_well_conditioned && pj->gradients.D_well_conditioned); + const unsigned char u_well_conditioned = + (pi->gradients.u_well_conditioned && pj->gradients.u_well_conditioned); + + /* Flag to revert to use high order gradients */ + unsigned char high_order_gradients_flag = 0; + + /* Always use high order gradients when both pi and pj are well-conditioned */ + if (C_well_conditioned && D_well_conditioned && u_well_conditioned) { + + /* Corrected gradients */ + hydro_real_t G_i[3] = {0., 0., 0.}; + hydro_real_t G_j[3] = {0., 0., 0.}; + /* Rosswog 2020 Eqs 4 & 5 use dx_ji for both */ + hydro_mat3x3_vec3_dot(pi->gradients.correction_matrix, dx_ji, G_i); + hydro_mat3x3_vec3_dot(pj->gradients.correction_matrix, dx_ji, G_j); + + G_i[0] *= wi * hi_inv_dim; + G_i[1] *= wi * hi_inv_dim; + G_i[2] *= wi * hi_inv_dim; + + G_j[0] *= wj * hj_inv_dim; + G_j[1] *= wj * hj_inv_dim; + G_j[2] *= wj * hj_inv_dim; + + /* Averaged correction gradient. Note: antisymmetric, so only need + * a sign flip for pj */ + hydro_get_average_kernel_gradient(pi, pj, G_i, G_j, G_ij); + + /* Compute from j perspective to ensure perfectly anti-symmetric */ + G_j[0] = 0.; + G_j[1] = 0.; + G_j[2] = 0.; + + G_i[0] = 0.; + G_i[1] = 0.; + G_i[2] = 0.; + + /* Swap G_i and G_j */ + hydro_mat3x3_vec3_dot(pj->gradients.correction_matrix, dx_ij, G_j); + hydro_mat3x3_vec3_dot(pi->gradients.correction_matrix, dx_ij, G_i); + + G_j[0] *= wj * hj_inv_dim; + G_j[1] *= wj * hj_inv_dim; + G_j[2] *= wj * hj_inv_dim; + + G_i[0] *= wi * hi_inv_dim; + G_i[1] *= wi * hi_inv_dim; + G_i[2] *= wi * hi_inv_dim; + + hydro_real_t G_ji[3] = {0., 0., 0.}; + hydro_get_average_kernel_gradient(pj, pi, G_j, G_i, G_ji); + + /* Average the two estimators */ + G_ij[0] = 0.5 * (G_ij[0] - G_ji[0]); + G_ij[1] = 0.5 * (G_ij[1] - G_ji[1]); + G_ij[2] = 0.5 * (G_ij[2] - G_ji[2]); + + /* Check if G_ij is extremely misaligned with the radial direction */ + G_ij_norm = hydro_vec3_norm(G_ij); + + /* Get G_ij along the separation vector */ + hydro_real_t G_ij_dot_dx_ij_hat = hydro_vec3_vec3_dot(G_ij, dx_ij_hat); + const hydro_real_t G_ij_dot_dx_ij_hat_abs = fabs(G_ij_dot_dx_ij_hat); + + /* Find the cos(theta) term between G and dx */ + hydro_real_t cosine_G_ij_dx_ij_hat = + G_ij_dot_dx_ij_hat_abs / (G_ij_norm + 1.e-10); + + /* Handle floating point errors */ + if (cosine_G_ij_dx_ij_hat > 1.) cosine_G_ij_dx_ij_hat = 1.; + + const unsigned char G_has_large_angle = + (cosine_G_ij_dx_ij_hat < const_viscosity_cosine_limit); + const unsigned char G_in_wrong_direction = (G_ij_dot_dx_ij_hat > 0.); + + /* Good angle between separation and correct direction, good to go! */ + if (!G_has_large_angle && !G_in_wrong_direction) { + + /* Make sure we use the correct interaction */ + high_order_gradients_flag = 1; + +#ifdef hydro_props_use_radial_artificial_terms + /* G along the separation vector */ + G_rad_ij[0] = G_ij_dot_dx_ij_hat * dx_ij_hat[0]; + G_rad_ij[1] = G_ij_dot_dx_ij_hat * dx_ij_hat[1]; + G_rad_ij[2] = G_ij_dot_dx_ij_hat * dx_ij_hat[2]; +#else + G_rad_ij[0] = G_ij[0]; + G_rad_ij[1] = G_ij[1]; + G_rad_ij[2] = G_ij[2]; +#endif + } + else { + /* Revert back to standard separation vector */ + G_ij[0] = dx[0]; + G_ij[1] = dx[1]; + G_ij[2] = dx[2]; + G_ij_norm = hydro_vec3_norm(G_ij); + + /* Use the separation vector for SPH */ + G_rad_ij[0] = dx[0]; + G_rad_ij[1] = dx[1]; + G_rad_ij[2] = dx[2]; + } + } + + /* MAGMA2-style gradients (Matrix-Inversion-2 SPH) */ + if (high_order_gradients_flag) { + + /* Compute second order reconstruction of velocity between pi & pj */ + hydro_real_t vi_reconstructed[3] = {0., 0., 0.}; + hydro_real_t vj_reconstructed[3] = {0., 0., 0.}; + + /* Important: Rosswog 2020 h_i and h_j are without kernel_gamma. Therefore, + * use xi and xj in the slope limiting procedure. */ + + /* Compute global Van Leer limiter (scalar, not component-wise) */ + const hydro_real_t phi_ij_vec = + hydro_vector_van_leer_phi(pi->gradients.velocity_tensor, + pj->gradients.velocity_tensor, + dx_ij, xi, xj); + + const hydro_real_t phi_ji_vec = + hydro_vector_van_leer_phi(pj->gradients.velocity_tensor, + pi->gradients.velocity_tensor, + dx_ji, xj, xi); + + /* Make sure no floating point problems */ + hydro_real_t phi_vec_sym = fmin(phi_ij_vec, phi_ji_vec); + phi_vec_sym = fmin(1., phi_vec_sym); + phi_vec_sym = fmax(0., phi_vec_sym); + + /* Need these recast in case of switching precision */ + const hydro_real_t v_i[3] = {pi->v[0], pi->v[1], pi->v[2]}; + const hydro_real_t v_j[3] = {pj->v[0], pj->v[1], pj->v[2]}; + + /* dx_ji for particle i and dx_ij for particle j */ + hydro_vector_second_order_reconstruction(phi_vec_sym, dx_ji, v_i, + pi->gradients.velocity_tensor, + pi->gradients.velocity_hessian, + vi_reconstructed); + + hydro_vector_second_order_reconstruction(phi_vec_sym, dx_ij, v_j, + pj->gradients.velocity_tensor, + pj->gradients.velocity_hessian, + vj_reconstructed); + + const hydro_real_t dv_reconstructed[3] = { + vi_reconstructed[0] - vj_reconstructed[0], + vi_reconstructed[1] - vj_reconstructed[1], + vi_reconstructed[2] - vj_reconstructed[2] + }; + + /* Get velocity difference, but limit reconstructed values */ + hydro_real_t dv_ij[3] = {0., 0., 0.}; + hydro_vec_minmod_limiter(dv_reconstructed, dv_raw, + pi->gradients.dv_min, pi->gradients.dv_max, + pj->gradients.dv_min, pj->gradients.dv_max, + dv_ij); + + /* Artificial viscosity */ + + + /* Get the acceleration term, depends on the weighting scheme */ + visc_acc_term = + hydro_get_visc_acc_term(pi, pj, dv_ij, dx_ij, fac_mu, a2_Hubble); + + /* Split heating between the two particles */ + visc_du_term = 0.5 * visc_acc_term; + + + /* Artificial conductivity */ + + + hydro_real_t ui_reconstructed = 0.; + hydro_real_t uj_reconstructed = 0.; + + + /* Compute global Van Leer limiter (scalar, not component-wise) */ + const hydro_real_t phi_ij_scalar = + hydro_scalar_van_leer_phi(pi->gradients.u, + pj->gradients.u, + dx_ij, xi, xj); + const hydro_real_t phi_ji_scalar = + hydro_scalar_van_leer_phi(pj->gradients.u, + pi->gradients.u, + dx_ji, xj, xi); + + /* Make sure no floating point problems */ + hydro_real_t phi_scalar_sym = fmin(phi_ij_scalar, phi_ji_scalar); + phi_scalar_sym = fmin(1., phi_scalar_sym); + phi_scalar_sym = fmax(0., phi_scalar_sym); + + /* dx_ji for particle i and dx_ij for particle j */ + hydro_scalar_second_order_reconstruction(phi_scalar_sym, dx_ji, + (hydro_real_t)pi->u, + pi->gradients.u, + pi->gradients.u_hessian, + &ui_reconstructed); + + hydro_scalar_second_order_reconstruction(phi_scalar_sym, dx_ij, + (hydro_real_t)pj->u, + pj->gradients.u, + pj->gradients.u_hessian, + &uj_reconstructed); + + const hydro_real_t rho_ij = 0.5 * (rhoi + rhoj); + const hydro_real_t rho_ij_inv = 1. / rho_ij; + const hydro_real_t dv_Hubble[3] = { + dv[0] + a2_Hubble * dx_ij[0], + dv[1] + a2_Hubble * dx_ij[1], + dv[2] + a2_Hubble * dx_ij[2] + }; + + const hydro_real_t dv_Hubble_dot_dx_ij_hat = + hydro_vec3_vec3_dot(dv_Hubble, dx_ij_hat); + + /* Limit art. cond. to only when information is communicable */ + const hydro_real_t c_ij = + 0.5 * (pi->force.soundspeed + pj->force.soundspeed); + const hydro_real_t v_sig_alpha = const_viscosity_alpha_prefactor * c_ij; + + /* Must connect the particles along the LOS */ + hydro_real_t mu_ij = fac_mu * dv_Hubble_dot_dx_ij_hat; + const hydro_real_t v_sig_beta = const_viscosity_beta_prefactor * mu_ij; + + /* Skip conduction if expansion beats sound speed along LOS */ + if (v_sig_alpha > v_sig_beta) { + + const hydro_real_t dv_ij_Hubble[3] = { + dv_ij[0] + a2_Hubble * dx_ij[0], + dv_ij[1] + a2_Hubble * dx_ij[1], + dv_ij[2] + a2_Hubble * dx_ij[2] + }; + + /* Signal velocity from speed contributions */ +#ifdef hydro_props_use_radial_artificial_terms + const hydro_real_t dv_ij_Hubble_dot_dx_ij_hat = + hydro_vec3_vec3_dot(dv_ij_Hubble, dx_ij_hat); + const hydro_real_t v_sig_speed = + fac_mu * fabs(dv_ij_Hubble_dot_dx_ij_hat); +#else + const hydro_real_t v_sig_speed = fac_mu * hydro_vec3_norm(dv_ij_Hubble); +#endif + + /* Get spec. energy difference, but limit reconstructed values */ + const hydro_real_t du_raw = pi->u - pj->u; + const hydro_real_t du_reconstructed = ui_reconstructed - uj_reconstructed; + const hydro_real_t du_ij = + hydro_scalar_minmod_limiter(du_reconstructed, du_raw, + pi->gradients.du_min, pi->gradients.du_max, + pj->gradients.du_min, pj->gradients.du_max); + + const hydro_real_t alpha_cond = const_conductivity_alpha; + const hydro_real_t delta_P = fabs(pressurei - pressurej); + const hydro_real_t P_lim = delta_P / (pressurei + pressurej); + + /* Add conductivity to the specific energy */ + cond_du_term = alpha_cond * P_lim * v_sig_speed * du_ij * rho_ij_inv; + } + else { + mu_ij = 0.; + } + + + /* Finalize everything with the correct normalizations. */ + + + /* Compute dv dot G_ij, reduces to dv dot dx in regular SPH. */ + const hydro_real_t dv_dot_G_ij = hydro_vec3_vec3_dot(dv_raw, G_ij); + +#ifdef hydro_props_use_radial_artificial_terms + /* Compute Hubble flow along LOS */ + const hydro_real_t dv_Hubble_along_dx_ij[3] = { + dv_Hubble_dot_dx_ij_hat * dx_ij_hat[0], + dv_Hubble_dot_dx_ij_hat * dx_ij_hat[1], + dv_Hubble_dot_dx_ij_hat * dx_ij_hat[2] + }; + + const hydro_real_t dv_Hubble_dot_G_rad_ij = + hydro_vec3_vec3_dot(dv_Hubble_along_dx_ij, G_rad_ij); +#else + const hydro_real_t dv_Hubble_dot_G_rad_ij = + hydro_vec3_vec3_dot(dv_Hubble, G_rad_ij); +#endif + + /* Evolve the heating terms using the velocity divergence */ + sph_du_term_i *= dv_dot_G_ij; + sph_du_term_j *= dv_dot_G_ij; + cond_du_term *= -G_rad_ij_norm; /* Eq. 24 Rosswog 2020 */ + visc_du_term *= dv_Hubble_dot_G_rad_ij; + + if (visc_du_term <= 0.) { + visc_acc_term = 0.; + visc_du_term = 0.; + } + + + /* Get the time derivative for h. */ + + + /* Velocity divergence is from the SPH estimator, not the G_ij vector */ + const hydro_real_t dv_dot_dx_ij = hydro_vec3_vec3_dot(dv_raw, dx_ij); + pi->force.h_dt -= hydro_get_h_dt_sum(dv_dot_dx_ij, dv_dot_G_ij, + mj, rhoj_inv, r_inv, wi_dr); + pj->force.h_dt -= hydro_get_h_dt_sum(dv_dot_dx_ij, dv_dot_G_ij, + mi, rhoi_inv, r_inv, wj_dr); + + + /* Timestepping */ + + + /* Compute based on raw velocities */ + const hydro_real_t v_sig_visc = + signal_velocity(dx, pi, pj, mu_ij, const_viscosity_beta); + + const hydro_real_t h_ij = 0.5 * (hi + hj); + pi->h_min = fmin(pi->h_min, h_ij); + pj->h_min = fmin(pj->h_min, h_ij); + + /* New timestep estimate */ + const hydro_real_t dt_min_i = h_ij / v_sig_visc; + const hydro_real_t dt_min_j = h_ij / v_sig_visc; + pi->dt_min = fmin(pi->dt_min, dt_min_i); + pj->dt_min = fmin(pj->dt_min, dt_min_j); + +#ifdef MAGMA2_DEBUG_CHECKS + pi->debug.N_force_high_order_grad++; + pj->debug.N_force_high_order_grad++; +#endif + } + else { /* Gasoline-like SPH fallback */ + + + /* Compute dv dot dr. */ + const hydro_real_t dvdr = hydro_vec3_vec3_dot(dv_raw, G_ij); + + /* Includes the hubble flow term; not used for du/dt */ + const hydro_real_t dvdr_Hubble = dvdr + a2_Hubble * r2; + + /* Are the particles moving towards each others ? */ + const hydro_real_t omega_ij = fmin(dvdr_Hubble, 0.); + + hydro_real_t mu_full_ij = fac_mu * r_inv * dvdr_Hubble; + const hydro_real_t mu_ij = fac_mu * r_inv * omega_ij; + + /* Construct the full viscosity term */ + const hydro_real_t b_ij = + 0.5 * (pi->gradients.balsara + pj->gradients.balsara); + const hydro_real_t rho_ij = rhoi + rhoj; + const hydro_real_t cs_ij = pi->force.soundspeed + pj->force.soundspeed; + const hydro_real_t c_ij = 0.5 * cs_ij; + const hydro_real_t alpha = const_viscosity_alpha; + const hydro_real_t visc = + omega_ij < 0. + ? (-0.25 * alpha * cs_ij * + mu_ij + + const_viscosity_beta * mu_ij * mu_ij) * b_ij / + (0.5 * rho_ij) + : 0.; + + visc_acc_term = const_fallback_reduction_factor * visc; + visc_du_term = 0.5 * visc_acc_term; + + const hydro_real_t v_sig_alpha = const_viscosity_alpha_prefactor * c_ij; + const hydro_real_t v_sig_beta = const_viscosity_beta_prefactor * mu_full_ij; + + if (v_sig_alpha > v_sig_beta) { + const hydro_real_t rho_ij_inv = 2. / rho_ij; + const hydro_real_t du = pi->u - pj->u; + + const hydro_real_t alpha_cond = const_conductivity_alpha; + const hydro_real_t delta_P = fabs(pressurei - pressurej); + const hydro_real_t P_lim = delta_P / (pressurei + pressurej); + + cond_du_term = alpha_cond * P_lim * fabs(mu_full_ij) * du * rho_ij_inv; + cond_du_term *= const_fallback_reduction_factor; + } + + /* New signal velocity */ + const hydro_real_t new_v_sig_visc = + signal_velocity(dx, pi, pj, mu_ij, const_viscosity_beta); + + const hydro_real_t h_ij = 0.5 * (hi + hj); + pi->h_min = fmin(pi->h_min, h_ij); + pj->h_min = fmin(pj->h_min, h_ij); + + /* New timestep estimate */ + const hydro_real_t dt_min_i = h_ij / new_v_sig_visc; + const hydro_real_t dt_min_j = h_ij / new_v_sig_visc; + pi->dt_min = fmin(pi->dt_min, dt_min_i); + pj->dt_min = fmin(pj->dt_min, dt_min_j); + + const hydro_real_t kernel_gradient = 0.5 * (wi_dr + wj_dr) * r_inv; + + visc_acc_term *= kernel_gradient; + sph_acc_term *= kernel_gradient; + + sph_du_term_i *= dvdr * kernel_gradient; + sph_du_term_j *= dvdr * kernel_gradient; + visc_du_term *= dvdr_Hubble * kernel_gradient; + cond_du_term *= kernel_gradient; + + /* Get the time derivative for h. */ + pi->force.h_dt -= hydro_get_h_dt_sum(dvdr, 0., mj, rhoj_inv, r_inv, wi_dr); + pj->force.h_dt -= hydro_get_h_dt_sum(dvdr, 0., mi, rhoi_inv, r_inv, wj_dr); + +#ifdef MAGMA2_DEBUG_CHECKS + pi->debug.N_force_low_order_grad++; + pj->debug.N_force_low_order_grad++; +#endif + } + + + /* Get the time derivative for v. */ + + + /* Assemble the acceleration */ + const hydro_real_t acc[3] = { + sph_acc_term * G_ij[0] + visc_acc_term * G_rad_ij[0], + sph_acc_term * G_ij[1] + visc_acc_term * G_rad_ij[1], + sph_acc_term * G_ij[2] + visc_acc_term * G_rad_ij[2] + }; + + /* Use the force Luke ! */ + pi->a_hydro[0] -= mj * acc[0]; + pi->a_hydro[1] -= mj * acc[1]; + pi->a_hydro[2] -= mj * acc[2]; + + pj->a_hydro[0] += mi * acc[0]; + pj->a_hydro[1] += mi * acc[1]; + pj->a_hydro[2] += mi * acc[2]; + + + /* Get the time derivative for u. */ + + + /* Assemble the energy equation term */ + const hydro_real_t du_dt_i = sph_du_term_i + visc_du_term + cond_du_term; + const hydro_real_t du_dt_j = sph_du_term_j + visc_du_term - cond_du_term; + + /* Internal energy time derivative */ + pi->u_dt += du_dt_i * mj; + pj->u_dt += du_dt_j * mi; + + pi->u_dt_cond += cond_du_term * mj; + pj->u_dt_cond -= cond_du_term * mi; +} + +/** + * @brief Force interaction between two particles (non-symmetric). + * + * @param r2 Comoving square distance between the two particles. + * @param dx Comoving vector separating both particles (pi - pj). + * @param hi Comoving smoothing-length of part*icle i. + * @param hj Comoving smoothing-length of part*icle j. + * @param pi First part*icle. + * @param pj Second part*icle (not updated). + * @param a Current scale factor. + * @param H Current Hubble parameter. + */ +__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( + const float r2, const float dx[3], const float hi, const float hj, + struct part* restrict pi, const struct part* restrict pj, const float a, + const float H) { + + /* Cosmological factors entering the EoMs */ + const hydro_real_t fac_mu = pow_three_gamma_minus_five_over_two(a); + const hydro_real_t a2_Hubble = a * a * H; + + const hydro_real_t r = sqrt(r2); + const hydro_real_t r_inv = r ? 1.0 / r : 0.0; + + /* Recover some data */ + const hydro_real_t mj = pj->mass; + + const hydro_real_t rhoi = pi->rho; + const hydro_real_t rhoj = pj->rho; + const hydro_real_t rhoi_inv = 1. / rhoi; + const hydro_real_t rhoj_inv = 1. / rhoj; + const hydro_real_t rhoij_inv = rhoi_inv * rhoj_inv; + + const hydro_real_t pressurei = pi->force.pressure; + const hydro_real_t pressurej = pj->force.pressure; + + /* Get the kernel for hi. */ + const hydro_real_t hi_inv = 1.0 / hi; + const hydro_real_t hi_inv_dim = pow_dimension(hi_inv); + const float xi = r * hi_inv; + float wi, wi_dx; + kernel_deval(xi, &wi, &wi_dx); + + /* Get the kernel for hj. */ + const hydro_real_t hj_inv = 1.0 / hj; + const hydro_real_t hj_inv_dim = pow_dimension(hj_inv); + const float xj = r * hj_inv; + float wj, wj_dx; + kernel_deval(xj, &wj, &wj_dx); + + /* For dh/dt and fall-back SPH */ + const hydro_real_t hi_inv_dim_plus_one = hi_inv * hi_inv_dim; /* 1/h^(d+1) */ + const hydro_real_t hj_inv_dim_plus_one = hj_inv * hj_inv_dim; + const hydro_real_t wi_dr = hi_inv_dim_plus_one * wi_dx; + const hydro_real_t wj_dr = hj_inv_dim_plus_one * wj_dx; + + /* Peculiar velocity difference vector */ + const hydro_real_t dv[3] = {pi->v[0] - pj->v[0], + pi->v[1] - pj->v[1], + pi->v[2] - pj->v[2]}; + const hydro_real_t dv_raw[3] = {dv[0], dv[1], dv[2]}; + + /* For MAGMA2, this is the full anti-symmetric gradient vector. For the + * fallback Gasoline2-style SPH, this will just be the direction vector + * between the two particles (r_i - r_j). */ + hydro_real_t G_ij[3] = {0., 0., 0.}; + hydro_real_t G_rad_ij[3] = {0., 0., 0.}; + hydro_real_t G_ij_norm = 0.; + hydro_real_t G_rad_ij_norm = 0.; + + /* Separation vectors and swapped */ + const hydro_real_t dx_ij[3] = {dx[0], dx[1], dx[2]}; + const hydro_real_t dx_ji[3] = {-dx[0], -dx[1], -dx[2]}; + hydro_real_t dx_ij_hat[3] = {0., 0., 0.}; + hydro_vec3_unit(dx_ij, dx_ij_hat); + + /* These are set whether or not we fall back onto SPH gradients */ + hydro_real_t visc_acc_term = 0.; + hydro_real_t sph_acc_term = (pressurei + pressurej) * rhoij_inv; + hydro_real_t sph_du_term_i = pressurei * rhoij_inv; + hydro_real_t visc_du_term = 0.; + hydro_real_t cond_du_term = 0.; + + /* Correction factor must be well-conditioned. */ + const unsigned char C_well_conditioned = + (pi->gradients.C_well_conditioned && pj->gradients.C_well_conditioned); + /* First order density-independent velocity gradients */ + const unsigned char D_well_conditioned = + (pi->gradients.D_well_conditioned && pj->gradients.D_well_conditioned); + const unsigned char u_well_conditioned = + (pi->gradients.u_well_conditioned && pj->gradients.u_well_conditioned); + + /* Flag to use second order gradients */ + unsigned char high_order_gradients_flag = 0; + + /* Always use SPH gradients between particles if one of them has an + * ill-conditioned C matrix */ + if (C_well_conditioned && D_well_conditioned && u_well_conditioned) { + + /* Corrected gradients */ + hydro_real_t G_i[3] = {0., 0., 0.}; + hydro_real_t G_j[3] = {0., 0., 0.}; + /* Rosswog 2020 Eqs 4 & 5 use dx_ji for both */ + hydro_mat3x3_vec3_dot(pi->gradients.correction_matrix, dx_ji, G_i); + hydro_mat3x3_vec3_dot(pj->gradients.correction_matrix, dx_ji, G_j); + + G_i[0] *= wi * hi_inv_dim; + G_i[1] *= wi * hi_inv_dim; + G_i[2] *= wi * hi_inv_dim; + + G_j[0] *= wj * hj_inv_dim; + G_j[1] *= wj * hj_inv_dim; + G_j[2] *= wj * hj_inv_dim; + + /* Averaged correction gradient. Note: antisymmetric, so only need + * a sign flip for pj */ + hydro_get_average_kernel_gradient(pi, pj, G_i, G_j, G_ij); + + /* Compute from j perspective to ensure perfectly anti-symmetric */ + G_j[0] = 0.; + G_j[1] = 0.; + G_j[2] = 0.; + + G_i[0] = 0.; + G_i[1] = 0.; + G_i[2] = 0.; + + /* Swap G_i and G_j */ + hydro_mat3x3_vec3_dot(pj->gradients.correction_matrix, dx_ij, G_j); + hydro_mat3x3_vec3_dot(pi->gradients.correction_matrix, dx_ij, G_i); + + G_j[0] *= wj * hj_inv_dim; + G_j[1] *= wj * hj_inv_dim; + G_j[2] *= wj * hj_inv_dim; + + G_i[0] *= wi * hi_inv_dim; + G_i[1] *= wi * hi_inv_dim; + G_i[2] *= wi * hi_inv_dim; + + hydro_real_t G_ji[3] = {0., 0., 0.}; + hydro_get_average_kernel_gradient(pj, pi, G_j, G_i, G_ji); + + /* Average the two estimators */ + G_ij[0] = 0.5 * (G_ij[0] - G_ji[0]); + G_ij[1] = 0.5 * (G_ij[1] - G_ji[1]); + G_ij[2] = 0.5 * (G_ij[2] - G_ji[2]); + + /* Check if G_ij is extremely misaligned with the radial direction */ + G_ij_norm = hydro_vec3_norm(G_ij); + + /* Get G_ij along the separation vector */ + hydro_real_t G_ij_dot_dx_ij_hat = hydro_vec3_vec3_dot(G_ij, dx_ij_hat); + const hydro_real_t G_ij_dot_dx_ij_hat_abs = fabs(G_ij_dot_dx_ij_hat); + + /* Find the cos(theta) term between G and dx */ + hydro_real_t cosine_G_ij_dx_ij_hat = + G_ij_dot_dx_ij_hat_abs / (G_ij_norm + 1.e-10); + + /* Handle floating point errors */ + if (cosine_G_ij_dx_ij_hat > 1.) cosine_G_ij_dx_ij_hat = 1.; + + const unsigned char G_has_large_angle = + (cosine_G_ij_dx_ij_hat < const_viscosity_cosine_limit); + const unsigned char G_in_wrong_direction = (G_ij_dot_dx_ij_hat > 0.); + + /* Good angle between separation and correct direction, good to go! */ + if (!G_has_large_angle && !G_in_wrong_direction) { + + /* Make sure we use the correct gradients below */ + high_order_gradients_flag = 1; + +#ifdef hydro_props_use_radial_artificial_terms + /* G along the separation vector */ + G_rad_ij[0] = G_ij_dot_dx_ij_hat * dx_ij_hat[0]; + G_rad_ij[1] = G_ij_dot_dx_ij_hat * dx_ij_hat[1]; + G_rad_ij[2] = G_ij_dot_dx_ij_hat * dx_ij_hat[2]; +#else + G_rad_ij[0] = G_ij[0]; + G_rad_ij[1] = G_ij[1]; + G_rad_ij[2] = G_ij[2]; +#endif + } + else { + /* Revert back to standard separation vector */ + G_ij[0] = dx[0]; + G_ij[1] = dx[1]; + G_ij[2] = dx[2]; + G_ij_norm = hydro_vec3_norm(G_ij); + + /* Use the separation vector for SPH */ + G_rad_ij[0] = dx[0]; + G_rad_ij[1] = dx[1]; + G_rad_ij[2] = dx[2]; + } + + G_rad_ij_norm = hydro_vec3_norm(G_rad_ij); + } + + /* MAGMA2 style gradients (Matrix-Inversion-2 SPH) */ + if (high_order_gradients_flag) { + + /* Compute second order reconstruction of velocity between pi & pj */ + hydro_real_t vi_reconstructed[3] = {0., 0., 0.}; + hydro_real_t vj_reconstructed[3] = {0., 0., 0.}; + + /* Important: Rosswog 2020 h_i and h_j are without kernel_gamma. Therefore, + * use xi and xj in the slope limiting procedure. */ + + /* Compute global Van Leer limiter (scalar, not component-wise) */ + const hydro_real_t phi_ij_vec = + hydro_vector_van_leer_phi(pi->gradients.velocity_tensor, + pj->gradients.velocity_tensor, + dx_ij, xi, xj); + const hydro_real_t phi_ji_vec = + hydro_vector_van_leer_phi(pj->gradients.velocity_tensor, + pi->gradients.velocity_tensor, + dx_ji, xj, xi); + + /* Make sure no floating point problems */ + hydro_real_t phi_vec_sym = fmin(phi_ij_vec, phi_ji_vec); + phi_vec_sym = fmin(1., phi_vec_sym); + phi_vec_sym = fmax(0., phi_vec_sym); + + /* Need these recast in case of switching precision */ + const hydro_real_t v_i[3] = {pi->v[0], pi->v[1], pi->v[2]}; + const hydro_real_t v_j[3] = {pj->v[0], pj->v[1], pj->v[2]}; + + /* dx_ji for particle i and dx_ij for particle j */ + hydro_vector_second_order_reconstruction(phi_vec_sym, dx_ji, v_i, + pi->gradients.velocity_tensor, + pi->gradients.velocity_hessian, + vi_reconstructed); + + hydro_vector_second_order_reconstruction(phi_vec_sym, dx_ij, v_j, + pj->gradients.velocity_tensor, + pj->gradients.velocity_hessian, + vj_reconstructed); + + const hydro_real_t dv_reconstructed[3] = { + vi_reconstructed[0] - vj_reconstructed[0], + vi_reconstructed[1] - vj_reconstructed[1], + vi_reconstructed[2] - vj_reconstructed[2] + }; + + /* Get velocity difference, but limit reconstructed values */ + hydro_real_t dv_ij[3] = {0., 0., 0.}; + hydro_vec_minmod_limiter(dv_reconstructed, dv_raw, + pi->gradients.dv_min, pi->gradients.dv_max, + pj->gradients.dv_min, pj->gradients.dv_max, + dv_ij); + + + /* Artificial viscosity */ + + + /* Get the acceleration term, depends on the weighting scheme */ + visc_acc_term = + hydro_get_visc_acc_term(pi, pj, dv_ij, dx_ij, fac_mu, a2_Hubble); + + /* Split heating between the two particles */ + visc_du_term = 0.5 * visc_acc_term; + + + /* Artificial conductivity */ + + + hydro_real_t ui_reconstructed = 0.; + hydro_real_t uj_reconstructed = 0.; + + /* Compute global Van Leer limiter (scalar, not component-wise) */ + const hydro_real_t phi_ij_scalar = + hydro_scalar_van_leer_phi(pi->gradients.u, + pj->gradients.u, + dx_ij, xi, xj); + const hydro_real_t phi_ji_scalar = + hydro_scalar_van_leer_phi(pj->gradients.u, + pi->gradients.u, + dx_ji, xj, xi); + + /* Make sure no floating point problems */ + hydro_real_t phi_scalar_sym = fmin(phi_ij_scalar, phi_ji_scalar); + phi_scalar_sym = fmin(1., phi_scalar_sym); + phi_scalar_sym = fmax(0., phi_scalar_sym); + + /* dx_ji for particle i and dx_ij for particle j */ + hydro_scalar_second_order_reconstruction(phi_scalar_sym, dx_ji, + (hydro_real_t)pi->u, + pi->gradients.u, + pi->gradients.u_hessian, + &ui_reconstructed); + + hydro_scalar_second_order_reconstruction(phi_scalar_sym, dx_ij, + (hydro_real_t)pj->u, + pj->gradients.u, + pj->gradients.u_hessian, + &uj_reconstructed); + + const hydro_real_t rho_ij = 0.5 * (rhoi + rhoj); + const hydro_real_t rho_ij_inv = 1. / rho_ij; + const hydro_real_t dv_Hubble[3] = { + dv[0] + a2_Hubble * dx_ij[0], + dv[1] + a2_Hubble * dx_ij[1], + dv[2] + a2_Hubble * dx_ij[2] + }; + + const hydro_real_t dv_Hubble_dot_dx_ij_hat = + hydro_vec3_vec3_dot(dv_Hubble, dx_ij_hat); + + /* Limit art. cond. to only when information is communicable */ + const hydro_real_t c_ij = + 0.5 * (pi->force.soundspeed + pj->force.soundspeed); + const hydro_real_t v_sig_alpha = const_viscosity_alpha_prefactor * c_ij; + + /* Must connect the particles along the LOS */ + hydro_real_t mu_ij = fac_mu * dv_Hubble_dot_dx_ij_hat; + const hydro_real_t v_sig_beta = const_viscosity_beta_prefactor * mu_ij; + + /* Skip conduction if expansion beats sound speed along LOS */ + if (v_sig_alpha > v_sig_beta) { + + const hydro_real_t dv_ij_Hubble[3] = { + dv_ij[0] + a2_Hubble * dx_ij[0], + dv_ij[1] + a2_Hubble * dx_ij[1], + dv_ij[2] + a2_Hubble * dx_ij[2] + }; + + /* Signal velocity from speed contributions */ +#ifdef hydro_props_use_radial_artificial_terms + const hydro_real_t dv_ij_Hubble_dot_dx_ij_hat = + hydro_vec3_vec3_dot(dv_ij_Hubble, dx_ij_hat); + const hydro_real_t v_sig_speed = + fac_mu * fabs(dv_ij_Hubble_dot_dx_ij_hat); +#else + const hydro_real_t v_sig_speed = fac_mu * hydro_vec3_norm(dv_ij_Hubble); +#endif + + /* Get spec. energy difference, but limit reconstructed values */ + const hydro_real_t du_raw = pi->u - pj->u; + const hydro_real_t du_reconstructed = ui_reconstructed - uj_reconstructed; + const hydro_real_t du_ij = + hydro_scalar_minmod_limiter(du_reconstructed, du_raw, + pi->gradients.du_min, pi->gradients.du_max, + pj->gradients.du_min, pj->gradients.du_max); + + const hydro_real_t alpha_cond = const_conductivity_alpha; + const hydro_real_t delta_P = fabs(pressurei - pressurej); + const hydro_real_t P_lim = delta_P / (pressurei + pressurej); + + /* Add conductivity to the specific energy */ + cond_du_term = alpha_cond * P_lim * v_sig_speed * du_ij * rho_ij_inv; + } + else { + mu_ij = 0.; + } + + /* Finalize the viscosity and conductivity with correct normalizations. */ + + + /* Compute dv dot G_ij, reduces to dv dot dx in regular SPH. */ + const hydro_real_t dv_dot_G_ij = hydro_vec3_vec3_dot(dv_raw, G_ij); + +#ifdef hydro_props_use_radial_artificial_terms + /* Get Hubble flow along dx */ + const hydro_real_t dv_Hubble_along_dx_ij[3] = { + dv_Hubble_dot_dx_ij_hat * dx_ij_hat[0], + dv_Hubble_dot_dx_ij_hat * dx_ij_hat[1], + dv_Hubble_dot_dx_ij_hat * dx_ij_hat[2] + }; + + /* Get Hubble flow contribution */ + const hydro_real_t dv_Hubble_dot_G_rad_ij = + hydro_vec3_vec3_dot(dv_Hubble_along_dx_ij, G_rad_ij); +#else + const hydro_real_t dv_Hubble_dot_G_rad_ij = + hydro_vec3_vec3_dot(dv_Hubble, G_rad_ij); +#endif + + sph_du_term_i *= dv_dot_G_ij; + cond_du_term *= -G_rad_ij_norm; /* Eq. 24 Rosswog 2020 */ + visc_du_term *= dv_Hubble_dot_G_rad_ij; + + if (visc_du_term <= 0.) { + visc_acc_term = 0.; + visc_du_term = 0.; + } + + + /* Get the time derivative for h. */ + + + const hydro_real_t dv_dot_dx_ij = hydro_vec3_vec3_dot(dv_raw, dx_ij); + pi->force.h_dt -= hydro_get_h_dt_sum(dv_dot_dx_ij, dv_dot_G_ij, + mj, rhoj_inv, r_inv, wi_dr); + + + /* Timestepping */ + + + /* Compute based on raw velocities */ + const hydro_real_t v_sig_visc = + signal_velocity(dx, pi, pj, mu_ij, const_viscosity_beta); + + const hydro_real_t h_ij = 0.5 * (hi + hj); + pi->h_min = fmin(pi->h_min, h_ij); + + const hydro_real_t dt_min_i = h_ij / v_sig_visc; + pi->dt_min = fmin(pi->dt_min, dt_min_i); + +#ifdef MAGMA2_DEBUG_CHECKS + pi->debug.N_force_high_order_grad++; +#endif + } + else { /* Gasoline2 SPH fallback*/ + + + /* Compute dv dot dr. */ + const hydro_real_t dvdr = hydro_vec3_vec3_dot(dv_raw, G_ij); + + /* Includes the hubble flow term; not used for du/dt */ + const hydro_real_t dvdr_Hubble = dvdr + a2_Hubble * r2; + + /* Are the particles moving towards each others ? */ + const hydro_real_t omega_ij = fmin(dvdr_Hubble, 0.); + + hydro_real_t mu_full_ij = fac_mu * r_inv * dvdr_Hubble; + const hydro_real_t mu_ij = fac_mu * r_inv * omega_ij; + + /* Construct the full viscosity term */ + const hydro_real_t b_ij = + 0.5 * (pi->gradients.balsara + pj->gradients.balsara); + const hydro_real_t rho_ij = rhoi + rhoj; + const hydro_real_t cs_ij = pi->force.soundspeed + pj->force.soundspeed; + const hydro_real_t c_ij = 0.5 * cs_ij; + const hydro_real_t alpha = const_viscosity_alpha; + const hydro_real_t visc = + omega_ij < 0. + ? (-0.25 * alpha * cs_ij * + mu_ij + + const_viscosity_beta * mu_ij * mu_ij) * b_ij / + (0.5 * rho_ij) + : 0.; + + visc_acc_term = const_fallback_reduction_factor * visc; + visc_du_term = 0.5 * visc_acc_term; + + const hydro_real_t v_sig_alpha = const_viscosity_alpha_prefactor * c_ij; + const hydro_real_t v_sig_beta = const_viscosity_beta_prefactor * mu_full_ij; + + if (v_sig_alpha > v_sig_beta) { + const hydro_real_t rho_ij_inv = 2. / rho_ij; + const hydro_real_t du = pi->u - pj->u; + + const hydro_real_t alpha_cond = const_conductivity_alpha; + const hydro_real_t delta_P = fabs(pressurei - pressurej); + const hydro_real_t P_lim = delta_P / (pressurei + pressurej); + + cond_du_term = alpha_cond * P_lim * fabs(mu_full_ij) * du * rho_ij_inv; + cond_du_term *= const_fallback_reduction_factor; + } + + /* New signal velocity */ + const hydro_real_t new_v_sig_visc = + signal_velocity(dx, pi, pj, mu_ij, const_viscosity_beta); + + const hydro_real_t h_ij = 0.5 * (hi + hj); + pi->h_min = fmin(pi->h_min, h_ij); + + /* New time-step estimate */ + const hydro_real_t dt_min_i = h_ij / new_v_sig_visc; + pi->dt_min = fmin(pi->dt_min, dt_min_i); + + /* Variable smoothing length term */ + const hydro_real_t kernel_gradient = 0.5 * (wi_dr + wj_dr) * r_inv; + + visc_acc_term *= kernel_gradient; + sph_acc_term *= kernel_gradient; + sph_du_term_i *= dvdr * kernel_gradient; + visc_du_term *= dvdr_Hubble * kernel_gradient; + cond_du_term *= kernel_gradient; + + /* Get the time derivative for h. */ + pi->force.h_dt -= hydro_get_h_dt_sum(dvdr, 0., mj, rhoj_inv, r_inv, wi_dr); + +#ifdef MAGMA2_DEBUG_CHECKS + pi->debug.N_force_low_order_grad++; +#endif + } + + /* Assemble the acceleration */ + const hydro_real_t acc[3] = { + sph_acc_term * G_ij[0] + visc_acc_term * G_rad_ij[0], + sph_acc_term * G_ij[1] + visc_acc_term * G_rad_ij[1], + sph_acc_term * G_ij[2] + visc_acc_term * G_rad_ij[2] + }; + + /* Use the force Luke ! */ + pi->a_hydro[0] -= mj * acc[0]; + pi->a_hydro[1] -= mj * acc[1]; + pi->a_hydro[2] -= mj * acc[2]; + + /* Assemble the energy equation term */ + const hydro_real_t du_dt_i = sph_du_term_i + visc_du_term + cond_du_term; + + /* Internal energy time derivative */ + pi->u_dt += du_dt_i * mj; + + pi->u_dt_cond += cond_du_term * mj; +} + +#endif /* SWIFT_MAGMA2_HYDRO_IACT_H */ diff --git a/src/hydro/MAGMA2/hydro_io.h b/src/hydro/MAGMA2/hydro_io.h new file mode 100644 index 0000000000..473eeb5b5d --- /dev/null +++ b/src/hydro/MAGMA2/hydro_io.h @@ -0,0 +1,346 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) & + * Matthieu Schaller (schaller@strw.leidenuniv.nl) + * 2025 Doug Rennehan (douglas.rennehan@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 . + * + ******************************************************************************/ +#ifndef SWIFT_MAGMA2_HYDRO_IO_H +#define SWIFT_MAGMA2_HYDRO_IO_H + +/** + * @file MAGMA2/hydro_io.h + * @brief Density-Energy non-conservative implementation of SPH, + * with added MAGMA2 physics (Rosswog 2020) (i/o routines) + */ + +#include "adiabatic_index.h" +#include "hydro.h" +#include "hydro_parameters.h" +#include "io_properties.h" +#include "kernel_hydro.h" + +/** + * @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. + */ +INLINE static 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, 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, u); + 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, rho); +} + +INLINE static void convert_S(const struct engine* e, const struct part* p, + const struct xpart* xp, float* ret) { + ret[0] = hydro_get_comoving_entropy(p, xp); +} + +INLINE static void convert_P(const struct engine* e, const struct part* p, + const struct xpart* xp, float* ret) { + ret[0] = hydro_get_comoving_pressure(p); +} + +INLINE static void convert_div_v(const struct engine* e, const struct part* p, + const struct xpart* xp, float* ret) { + ret[0] = p->gradients.velocity_tensor[0][0] + + p->gradients.velocity_tensor[1][1] + + p->gradients.velocity_tensor[2][2]; +} + +INLINE static void convert_part_pos(const struct engine* e, + const struct part* p, + const struct xpart* xp, double* ret) { + const struct space* s = e->s; + if (s->periodic) { + ret[0] = box_wrap(p->x[0], 0.0, s->dim[0]); + ret[1] = box_wrap(p->x[1], 0.0, s->dim[1]); + ret[2] = box_wrap(p->x[2], 0.0, s->dim[2]); + } else { + ret[0] = p->x[0]; + ret[1] = p->x[1]; + ret[2] = p->x[2]; + } + if (e->snapshot_use_delta_from_edge) { + ret[0] = min(ret[0], s->dim[0] - e->snapshot_delta_from_edge); + ret[1] = min(ret[1], s->dim[1] - e->snapshot_delta_from_edge); + ret[2] = min(ret[2], s->dim[2] - e->snapshot_delta_from_edge); + } +} + +INLINE static 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 float dt_kick_grav_mesh = e->dt_kick_grav_mesh_for_io; + + 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 term)*/ + ret[0] = xp->v_full[0] + p->a_hydro[0] * dt_kick_hydro; + ret[1] = xp->v_full[1] + p->a_hydro[1] * dt_kick_hydro; + ret[2] = xp->v_full[2] + p->a_hydro[2] * dt_kick_hydro; + + /* Add the gravity term */ + if (p->gpart != NULL) { + ret[0] += p->gpart->a_grav[0] * dt_kick_grav; + ret[1] += p->gpart->a_grav[1] * dt_kick_grav; + ret[2] += p->gpart->a_grav[2] * dt_kick_grav; + } + + /* And the mesh gravity term */ + if (p->gpart != NULL) { + ret[0] += p->gpart->a_grav_mesh[0] * dt_kick_grav_mesh; + ret[1] += p->gpart->a_grav_mesh[1] * dt_kick_grav_mesh; + ret[2] += p->gpart->a_grav_mesh[2] * dt_kick_grav_mesh; + } + + /* Conversion from internal units to peculiar velocities */ + ret[0] *= cosmo->a_inv; + ret[1] *= cosmo->a_inv; + ret[2] *= cosmo->a_inv; +} + +INLINE static 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. + */ +INLINE static void hydro_write_particles(const struct part* parts, + const struct xpart* xparts, + struct io_props* list, + int* num_fields) { + *num_fields = 0; + int num = 0; + + /* List what we want to write */ + list[num] = io_make_output_field_convert_part( + "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, parts, xparts, + convert_part_pos, "Co-moving positions of the particles"); + num++; + + list[num] = io_make_output_field_convert_part( + "Velocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, parts, xparts, + convert_part_vel, + "Peculiar velocities of the stars. This is (a * dx/dt) where x is the " + "co-moving positions of the particles"); + num++; + + list[num] = io_make_output_field( + "Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f, parts, mass, + "Masses of the particles"); + num++; + + list[num] = io_make_output_field( + "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, parts, h, + "Co-moving smoothing lengths (FWHM of the kernel) of the particles"); + num++; + + list[num] = io_make_output_field( + "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, + -3.f * hydro_gamma_minus_one, parts, u, + "Co-moving thermal energies per unit mass of the particles"); + num++; + + list[num] = io_make_physical_output_field( + "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id, + /*can convert to comoving=*/0, "Unique IDs of the particles"); + num++; + + list[num] = io_make_output_field( + "Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f, parts, rho, + "Co-moving mass densities of the particles"); + num++; + + list[num] = io_make_output_field_convert_part( + "Entropies", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, 0.f, parts, + xparts, convert_S, "Co-moving entropies per unit mass of the particles"); + num++; + + list[num] = io_make_output_field_convert_part( + "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, -3.f * hydro_gamma, parts, + xparts, convert_P, "Co-moving pressures of the particles"); + num++; + + list[num] = io_make_output_field_convert_part( + "VelocityDivergences", FLOAT, 1, UNIT_CONV_FREQUENCY, 0.f, parts, xparts, + convert_div_v, + "Local velocity divergence field around the particles. Provided without " + "cosmology, as this includes the Hubble flow. To return to a peculiar " + "velocity divergence, div . v_pec = a^2 (div . v - n_D H)"); + num++; + + list[num] = io_make_output_field_convert_part( + "Potentials", FLOAT, 1, UNIT_CONV_POTENTIAL, -1.f, parts, xparts, + convert_part_potential, + "Co-moving gravitational potential at position of the particles"); + num++; + +#ifdef MAGMA2_DEBUG_CHECKS + list[num] = io_make_output_field( + "MaxViscSignalVel", FLOAT, 1, UNIT_CONV_SPEED, + -1.5 * hydro_gamma_minus_one, parts, debug.v_sig_visc_max, + "Maximum viscous signal speed encountered by the particle this step."); + num++; + + list[num] = io_make_output_field( + "MaxCondSignalVel", FLOAT, 1, UNIT_CONV_SPEED, + -1.5 * hydro_gamma_minus_one, parts, debug.v_sig_cond_max, + "Maximum conduction signal speed encountered by the particle this step."); + num++; + + list[num] = io_make_output_field( + "NumberOfNeighbours", INT, 1, UNIT_CONV_NO_UNITS, + 0.f, parts, debug.num_ngb, + "Number of neighbours."); + num++; + + list[num] = io_make_output_field( + "LowOrderGradientsCount", INT, 1, UNIT_CONV_NO_UNITS, + 0.f, parts, debug.N_force_low_order_grad, + "Cumulative number of low order gradient force interactions."); + num++; + + list[num] = io_make_output_field( + "HighOrderGradientsCount", INT, 1, UNIT_CONV_NO_UNITS, + 0.f, parts, debug.N_force_high_order_grad, + "Cumulative number of high order gradient force interactions."); + num++; + + list[num] = io_make_output_field( + "CorrectionMatrices", FLOAT, 9, UNIT_CONV_LENGTH * UNIT_CONV_LENGTH, + 2.f, parts, debug.correction_matrix, + "Co-moving correction matrices for the particles."); + num++; + + list[num] = io_make_output_field( + "CorrectionIllConditionedCounts", INT, 1, UNIT_CONV_NO_UNITS, 0.f, parts, + debug.C_ill_conditioned_count, + "Count for how many times this particle had an ill-conditioned C matrix"); + num++; + + list[num] = io_make_output_field( + "VelocityGradientNumeratorMatrices", FLOAT, 9, + UNIT_CONV_DENSITY * UNIT_CONV_SPEED / UNIT_CONV_LENGTH, + -5.f, parts, debug.velocity_tensor_aux, + "Co-moving numerator matrices for the particles."); + num++; + + list[num] = io_make_output_field( + "VelocityGradientDenominatorMatrices", FLOAT, 9, UNIT_CONV_DENSITY, + -3.f, parts, debug.velocity_tensor_aux_norm, + "Co-moving denominator matrices for the particles."); + num++; + + list[num] = io_make_output_field( + "VelocityGradientIllConditionedCounts", INT, 1, UNIT_CONV_NO_UNITS, + 0.f, parts, debug.D_ill_conditioned_count, + "Count for how many times this particle had an ill-conditioned D matrix"); + num++; + + list[num] = io_make_output_field( + "SpecificEnergyGradientNumerators", FLOAT, 3, + UNIT_CONV_DENSITY * UNIT_CONV_ENERGY_PER_UNIT_MASS / UNIT_CONV_LENGTH, + -6.f, parts, debug.u_aux, + "Co-moving specific energy numerator matrices for the particles."); + num++; + + list[num] = io_make_output_field( + "SpecificEnergyGradientDenominators", FLOAT, 3, UNIT_CONV_DENSITY, + -3.f, parts, debug.u_aux, + "Co-moving specific energy gradient numerator for the particles."); + num++; + + list[num] = io_make_output_field( + "SpecificEnergyIllConditionedCounts", INT, 1, UNIT_CONV_NO_UNITS, + 0.f, parts, debug.u_ill_conditioned_count, + "Count for how many times this particle had an ill-conditioned u_aux_norm"); + num++; +#endif + + *num_fields = num; +} + +/** + * @brief Writes the current model of SPH to the file + * @param h_grpsph The HDF5 group in which to write + */ +INLINE static void hydro_write_flavour(hid_t h_grpsph) { + /* Viscosity and thermal conduction */ + /* Nothing in this minimal model... */ + io_write_attribute_s(h_grpsph, "Thermal Conductivity Model", + "Simple treatment as in Price (2008)"); + io_write_attribute_s(h_grpsph, "Viscosity Model", + "Gingold & Monaghan (1983)"); +} + +/** + * @brief Are we writing entropy in the internal energy field ? + * + * @return 1 if entropy is in 'internal energy', 0 otherwise. + */ +INLINE static int writeEntropyFlag(void) { return 0; } + +#endif /* SWIFT_MAGMA2_HYDRO_IO_H */ diff --git a/src/hydro/MAGMA2/hydro_parameters.h b/src/hydro/MAGMA2/hydro_parameters.h new file mode 100644 index 0000000000..3e99e97207 --- /dev/null +++ b/src/hydro/MAGMA2/hydro_parameters.h @@ -0,0 +1,297 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) + * 2025 Doug Rennehan (douglas.rennehan@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 . + * + ******************************************************************************/ + +#ifndef SWIFT_MAGMA2_HYDRO_PARAMETERS_H +#define SWIFT_MAGMA2_HYDRO_PARAMETERS_H + +/* Configuration file */ +#include + +/* Global headers */ +#if defined(HAVE_HDF5) +#include +#endif + +/* Local headers */ +#include "common_io.h" +#include "error.h" +#include "inline.h" + +/** + * @file MAGMA2/hydro_parameters.h + * @brief Density-Energy non-conservative implementation of SPH, + * with added MAGMA2 physics (Rosswog 2020) (default compile-time + * parameters). + * + * This file defines a number of things that are used in + * hydro_properties.c as defaults for run-time parameters + * as well as a number of compile-time parameters. + */ + + +/* ---------- Viscosity & Conductivitiy parameters ---------- */ + + +/*! Alpha viscosity, usually =1.0. For lower N_ngb, should be higher */ +#define const_viscosity_alpha 2.0 + +/*! Alpha conductivity, usually =0.05. At lower N_ngb, should be higher */ +#define const_conductivity_alpha 0.075 + +/*! Desired number of neighbours -- CRITICAL that this matches hydro props */ +#if defined(HYDRO_DIMENSION_1D) +#define const_kernel_target_neighbours 8.0 +#elif defined(HYDRO_DIMENSION_2D) +#define const_kernel_target_neighbours 34.0 +#else +#define const_kernel_target_neighbours 114.0 +#endif + + +/* ---------- These parameters should not be changed ---------- */ + +/*! Use a Swift-like estimator for dh/dt rather than the correct formula + * 0 = Simple mass flow estimator + * 1 = Correct formula based on number density constraint + * 2 = Using v_ij dot G_ij with simple mass flow estimator + */ +#define hydro_props_dh_dt_estimator_type 0 + +/*! Flag to use Balsara limiter */ +#define hydro_props_use_balsara_limiter + +/*! Flag to use additional slope limiting procedures */ +//#define hydro_props_use_extra_slope_limiter + +/* Flag to disallow sign flip in reconstructed quantities */ +//#define hydro_props_use_strict_minmod_limiter + +/* Slope limiter length, fraction of max. distance in kernel */ +#ifdef hydro_props_use_extra_slope_limiter +#define const_grad_overshoot_length 0.25 + +/*! Slope limiter tolerance */ +#define const_grad_overshoot_tolerance 0.0 +#endif + +/* Viscosity floor when G_ij is extremely misaligned with dx_ij */ +#define const_viscosity_cosine_limit 0.1736 + +/* Viscosity weighting scheme: + * 0 = (rho_i * q_i + rho_j * q_j) / (rho_i * rho_j) + * 1 = (rho_i * q_ij + rho_j * q_ij) / (rho_i * rho_j) + * 2 = 2.0 * q_ij / (rho_i + rho_j) */ +#define hydro_props_viscosity_weighting_type 2 + +/* Flag to use radial gradients for viscosity and conductivity */ +//#define hydro_props_use_radial_artificial_terms + +/*! Use the correction terms to make the internal energy match the mass flux */ +//#define hydro_props_use_adiabatic_correction + +/* Kernel gradient weighting scheme: + * 0 = 0.5 * (G_i + G_j) + * 1 = 0.5 * (f_i * G_i + f_j * G_j) + * 2 = 0.5 * f_ij * (G_i + G_j) + * with f_ij = 0.5 * (f_i + f_j) + * 3 = 0.5 * f_ij * (G_i + G_j) + * with f_ij = 2 * f_i * f_j / (f_i + f_j) + * 4 = 0.5 * f_ij * (G_i + G_j) + * with f_ij = sqrt(f_i * f_j) + * 5 = 0.5 * f_ij * (G_i + G_j) + * with f_ij = (f_i * rho_i + f_j * rho_j) / (rho_i + rho_j) + * 6 = 0.5 * f_ij * (G_i + G_j) + * with f_ij = 2 * f_i * f_j / (f_i * rho_i + f_j * rho_j) + * 7 = 0.5 * f_ij * (G_i + G_j) + * with f_ij = (f_i * P_i + f_j * P_j) / (P_i + P_j) + * 8 = 0.5 * f_ij * (G_i + G_j) + * with f_ij = 2 * f_i * f_j / (f_i * P_i + f_j * P_j) + * 9 = 0.5 * f_ij * (G_i + G_j) + * with f_ij = (f_i * V_i + f_j * V_j) / (V_i + V_j) + */ +#define hydro_props_kernel_gradient_weighting 0 + +/*! Use double precision for all matrix/vector operations */ +//#define hydro_props_use_double_precision + +#ifdef hydro_props_use_double_precision +/*! Consider matrix inversion to be ill-conditioned above this limit */ +#define const_condition_number_upper_limit 300. +/*! Mean interparticle spacing for this kernel and neighbour number */ +#define const_kernel_mean_spacing (kernel_gamma*pow(4. * M_PI / (3. * \ + (double)const_kernel_target_neighbours), 1. / 3.)) +#else +/*! Consider matrix inversion to be ill-conditioned above this limit */ +#define const_condition_number_upper_limit 60. +/*! Mean interparticle spacing for this kernel and neighbour number */ +#define const_kernel_mean_spacing (kernel_gamma*powf(4. * M_PI / (3. * \ + (float)const_kernel_target_neighbours), 1. / 3.)) +#endif + +/*! eta_crit Rosswog 2020 Eq 23. Of order the mean interparticle spacing. */ +#define const_slope_limiter_eta_crit (const_kernel_mean_spacing) + +/*! eta_fold from Frontiere+'17 Equation 51 */ +#define const_slope_limiter_eta_fold 0.2 + +/*! Softening squared (epsilon^2) in Eq. 15 Rosswog 2020 */ +#define const_viscosity_epsilon2 0.01 + +/*! Cosmology default const_viscosity_beta=2*const_viscosity_alpha + * Beta is defined as in e.g. Price (2010) Eqn (103) */ +#define const_viscosity_beta (2.0*const_viscosity_alpha) + +/*! Prefactor for alpha term in signal velocity */ +#define const_viscosity_alpha_prefactor (1.25 * (1. + \ + 0.75 * const_viscosity_alpha)) + +/*! Prefactor for beta term in signal velocity */ +#define const_viscosity_beta_prefactor (1.25 * 0.75 * const_viscosity_beta) + +/*! Fallback multiplier for alpha/beta terms to reduce spread */ +#define const_fallback_reduction_factor 0.25 + +/* ---------- Structures for below ---------- */ + + +/*! Artificial viscosity parameters */ +struct viscosity_global_data { }; + +/*! Thermal diffusion parameters */ +struct diffusion_global_data { }; + +/* Functions for reading from parameter file */ + +/* Forward declartions */ +struct swift_params; +struct phys_const; +struct unit_system; + +/* Define float or double depending on hydro_props_use_double_precision */ +#if defined(hydro_props_use_double_precision) +typedef double hydro_real_t; +#else +typedef float hydro_real_t; +#endif + +/* Viscosity */ + +/** + * @brief Initialises the viscosity parameters in the struct from + * the parameter file, or sets them to defaults. + * + * @param params: the pointer to the swift_params file + * @param unit_system: pointer to the unit system + * @param phys_const: pointer to the physical constants system + * @param viscosity: pointer to the viscosity_global_data struct to be filled. + **/ +static INLINE void viscosity_init(struct swift_params* params, + const struct unit_system* us, + const struct phys_const* phys_const, + struct viscosity_global_data* viscosity) { } + +/** + * @brief Initialises a viscosity struct to sensible numbers for mocking + * purposes. + * + * @param viscosity: pointer to the viscosity_global_data struct to be filled. + **/ +static INLINE void viscosity_init_no_hydro( + struct viscosity_global_data* viscosity) { } + +/** + * @brief Prints out the viscosity parameters at the start of a run. + * + * @param viscosity: pointer to the viscosity_global_data struct found in + * hydro_properties + **/ +static INLINE void viscosity_print( + const struct viscosity_global_data* viscosity) { + message("Artificial viscosity alpha set to %.3f", const_viscosity_alpha); + message("Artificial viscosity beta set to %.3f", const_viscosity_beta); +} + +#if defined(HAVE_HDF5) +/** + * @brief Prints the viscosity information to the snapshot when writing. + * + * @param h_grpsph: the SPH group in the ICs to write attributes to. + * @param viscosity: pointer to the viscosity_global_data struct. + **/ +static INLINE void viscosity_print_snapshot( + hid_t h_grpsph, const struct viscosity_global_data* viscosity) { + + io_write_attribute_f(h_grpsph, "Alpha viscosity", const_viscosity_alpha); + io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta); +} +#endif + +/* Diffusion */ + +/** + * @brief Initialises the diffusion parameters in the struct from + * the parameter file, or sets them to defaults. + * + * @param params: the pointer to the swift_params file + * @param unit_system: pointer to the unit system + * @param phys_const: pointer to the physical constants system + * @param diffusion_global_data: pointer to the diffusion struct to be filled. + **/ +static INLINE void diffusion_init(struct swift_params* params, + const struct unit_system* us, + const struct phys_const* phys_const, + struct diffusion_global_data* diffusion) { } + +/** + * @brief Initialises a diffusion struct to sensible numbers for mocking + * purposes. + * + * @param diffusion: pointer to the diffusion_global_data struct to be filled. + **/ +static INLINE void diffusion_init_no_hydro( + struct diffusion_global_data* diffusion) { } + +/** + * @brief Prints out the diffusion parameters at the start of a run. + * + * @param diffusion: pointer to the diffusion_global_data struct found in + * hydro_properties + **/ +static INLINE void diffusion_print( + const struct diffusion_global_data* diffusion) { + message("Artificial conductivity alpha set to %.3f", + const_conductivity_alpha); +} + +#ifdef HAVE_HDF5 +/** + * @brief Prints the diffusion information to the snapshot when writing. + * + * @param h_grpsph: the SPH group in the ICs to write attributes to. + * @param diffusion: pointer to the diffusion_global_data struct. + **/ +static INLINE void diffusion_print_snapshot( + hid_t h_grpsph, const struct diffusion_global_data* diffusion) { + io_write_attribute_f(h_grpsph, "Conductivity alpha", + const_conductivity_alpha); +} +#endif + +#endif /* SWIFT_MAGMA2_HYDRO_PARAMETERS_H */ diff --git a/src/hydro/MAGMA2/hydro_part.h b/src/hydro/MAGMA2/hydro_part.h new file mode 100644 index 0000000000..29779f2e83 --- /dev/null +++ b/src/hydro/MAGMA2/hydro_part.h @@ -0,0 +1,348 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2019 Josh Borrow (joshua.borrow@durham.ac.uk) & + * Matthieu Schaller (schaller@strw.leidenuniv.nl) + * 2025 Doug Rennehan (douglas.rennehan@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 . + * + ******************************************************************************/ +#ifndef SWIFT_MAGMA2_HYDRO_PART_H +#define SWIFT_MAGMA2_HYDRO_PART_H + +/** + * @file MAGMA2/hydro_part.h + * @brief Density-Energy non-conservative implementation of SPH, + * with added MAGMA2 physics (Rosswog 2020) (particle definition) + */ + +#include "adaptive_softening_struct.h" +#include "black_holes_struct.h" +#include "chemistry_struct.h" +#include "cooling_struct.h" +#include "csds.h" +#include "feedback_struct.h" +#include "hydro_parameters.h" +#include "mhd_struct.h" +#include "particle_splitting_struct.h" +#include "pressure_floor_struct.h" +#include "rt_struct.h" +#include "sink_struct.h" +#include "star_formation_struct.h" +#include "timestep_limiter_struct.h" +#include "tracers_struct.h" + +/** + * @brief Particle fields not needed during the SPH loops over neighbours. + * + * This structure contains the particle fields that are not used in the + * density or force loops. Quantities should be used in the kick, drift and + * potentially ghost tasks only. + */ +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 end of the last step */ + float a_grav[3]; + + /*! Internal energy at the last full step. */ + float u_full; + + /*! Additional data used to record particle splits */ + struct particle_splitting_data split_data; + + /*! Additional data used to record cooling information */ + struct cooling_xpart_data cooling_data; + + /* Additional data used by the tracers */ + struct tracers_xpart_data tracers_data; + + /* Additional data used by the tracers */ + struct star_formation_xpart_data sf_data; + + /* Additional data used by the feedback */ + struct feedback_xpart_data feedback_data; + + /*! Additional data used by the MHD scheme */ + struct mhd_xpart_data mhd_data; + +#ifdef WITH_CSDS + /* Additional data for the particle csds */ + struct csds_part_data csds_data; +#endif + +} SWIFT_STRUCT_ALIGN; + +/** + * @brief Particle fields for the SPH particles + * + * The density and force substructures are used to contain variables only used + * within the density and force loops over neighbours. All more permanent + * variables should be declared in the main part of the part structure, + */ +struct part { + /*! Particle unique ID. */ + long long id; + + /*! Pointer to corresponding gravity part. */ + struct gpart* gpart; + + /*! Particle position. */ + double x[3]; + + /*! Particle predicted velocity. */ + float v[3]; + + /*! Particle acceleration. */ + float a_hydro[3]; + + /*! Particle mass. */ + float mass; + + /*! Particle smoothing length. */ + float h; + + /*! Particle internal energy. */ + float u; + + /*! Time derivative of the internal energy. */ + float u_dt; + + /*! Particle density. */ + float rho; + + /*! Minimum smoothing length in the kernel */ + float h_min; + + /*! Minimum time-step amongst neighbours */ + float dt_min; + + /*! Conduction du/dt */ + float u_dt_cond; + + #ifdef MAGMA2_DEBUG_CHECKS + struct { + /*! Correction matrix at the last time it was ill-conditioned */ + hydro_real_t correction_matrix[3][3]; + + /*! Velocity tensor at ill-condition time */ + hydro_real_t velocity_tensor_aux[3][3]; + + /*! Velocity tensor norm ill-conditioned */ + hydro_real_t velocity_tensor_aux_norm[3][3]; + + /*! u_aux tensor at ill-condition time */ + hydro_real_t u_aux[3]; + + /*! u_aux_norm tensor ill-conditioned */ + hydro_real_t u_aux_norm[3]; + + /*! Number of times correction_matrix was ill-conditioned */ + int C_ill_conditioned_count; + + /*! Number of times velocity_tensor_aux_norm was ill-conditioned */ + int D_ill_conditioned_count; + + /*! Number of times u_aux_norm was ill-conditioned */ + int u_ill_conditioned_count; + + /*! How many low-order SPH gradients in force interactions */ + int N_force_low_order_grad; + + /*! How many high-order SPH gradients in force interactions */ + int N_force_high_order_grad; + + /*! Number of neighbors in the kernel */ + int num_ngb; + + /*! The maximum viscous signal speed this step */ + hydro_real_t v_sig_visc_max; + + /*! The maximum conductive signal speed this step */ + hydro_real_t v_sig_cond_max; + } debug; +#endif + + /* Store gradients in a separate struct */ + struct { +#ifdef hydro_props_use_adiabatic_correction + /*! Adiabatic kernel correction factor numerator */ + hydro_real_t adiabatic_f_numerator; + + /*! Adiabatic kernel correction factor denominator */ + hydro_real_t adiabatic_f_denominator; +#endif + + /*! Sum of the kernel weights */ + hydro_real_t wcount; + + /*! Correction matrix (C^ki in Rosswog 2020) */ + hydro_real_t correction_matrix[3][3]; + + /*! Flag for whether C is ill-conditioned */ + char C_well_conditioned; + + /*! Full velocity gradient tensor */ + hydro_real_t velocity_tensor[3][3]; + + /*! Auxiliary full velocity gradient tensor */ + hydro_real_t velocity_tensor_aux[3][3]; + + /*! Flag for whether D (velocity_tensor_aux) is ill-conditioned */ + char D_well_conditioned; + + /*! Normalization for computing velocity_tensor_aux */ + hydro_real_t velocity_tensor_aux_norm[3][3]; + + /*! Hessian tensor */ + hydro_real_t velocity_hessian[3][3][3]; + + /*! Internal energy gradient */ + hydro_real_t u[3]; + + /*! Auxiliary internal energy gradient */ + hydro_real_t u_aux[3]; + + /*! Normalization for computing u_aux */ + hydro_real_t u_aux_norm[3]; + + /*! Flag for whether u_aux_norm is ill-conditioned */ + char u_well_conditioned; + + /*! Internal energy Hessian */ + hydro_real_t u_hessian[3][3]; + + /*! Minimum delta u across kernel */ + hydro_real_t du_min; + + /*! Maximum delta u across kernel */ + hydro_real_t du_max; + + /*! Minimum dv across kernel */ + hydro_real_t dv_min[3]; + + /*! Maximum dv across kernel */ + hydro_real_t dv_max[3]; + + /*! Kernel size for slope limiting */ + hydro_real_t kernel_size; + + /*! Balsara limiter for divergences */ + hydro_real_t balsara; + + } gradients; + + /* Store density/force specific stuff. */ + union { + /** + * @brief Structure for the variables only used in the density loop over + * neighbours. + * + * Quantities in this sub-structure should only be accessed in the density + * loop over neighbours and the ghost task. + */ + struct { + /*! Neighbour number count. */ + float wcount; + + /*! Derivative of the neighbour number with respect to h. */ + float wcount_dh; + + /*! Derivative of density with respect to h */ + float rho_dh; + + } density; + + /** + * @brief Structure for the variables only used in the force loop over + * neighbours. + * + * Quantities in this sub-structure should only be accessed in the force + * loop over neighbours and the ghost, drift and kick tasks. + */ + struct { + /*! "Grad h" term -- only partial in P-U */ + float f; + + /*! Particle pressure. */ + float pressure; + + /*! Particle soundspeed. */ + float soundspeed; + + /*! Time derivative of smoothing length */ + float h_dt; + + } force; + }; + + /*! Additional data used for adaptive softening */ + struct adaptive_softening_part_data adaptive_softening_data; + + /*! Additional data used by the MHD scheme */ + struct mhd_part_data mhd_data; + + /*! Chemistry information */ + struct chemistry_part_data chemistry_data; + + /*! Cooling information */ + struct cooling_part_data cooling_data; + + /*! Additional data used by the feedback */ + struct feedback_part_data feedback_data; + + /*! Black holes information (e.g. swallowing ID) */ + struct black_holes_part_data black_holes_data; + + /*! Sink information (e.g. swallowing ID) */ + struct sink_part_data sink_data; + + /*! Additional data used by the pressure floor */ + struct pressure_floor_part_data pressure_floor_data; + + /*! Additional Radiative Transfer Data */ + struct rt_part_data rt_data; + + /*! RT sub-cycling time stepping data */ + struct rt_timestepping_data rt_time_data; + + /*! Time-step length */ + timebin_t time_bin; + + /*! Tree-depth at which size / 2 <= h * gamma < size */ + char depth_h; + + /*! Time-step limiter information */ + struct timestep_limiter_data limiter_data; + +#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_MAGMA2_HYDRO_PART_H */ diff --git a/src/hydro_csds.h b/src/hydro_csds.h index a0593ca7b3..088f9f6341 100644 --- a/src/hydro_csds.h +++ b/src/hydro_csds.h @@ -53,6 +53,8 @@ #include "./hydro/SPHENIX/hydro_csds.h" #elif defined(GASOLINE_SPH) #error TODO +#elif defined(MAGMA2_SPH) +#error TODO #elif defined(ANARCHY_PU_SPH) #error TODO #else diff --git a/src/hydro_io.h b/src/hydro_io.h index 5a64a284cc..5920dd0a9b 100644 --- a/src/hydro_io.h +++ b/src/hydro_io.h @@ -49,6 +49,8 @@ #include "./hydro/SPHENIX/hydro_io.h" #elif defined(GASOLINE_SPH) #include "./hydro/Gasoline/hydro_io.h" +#elif defined(MAGMA2_SPH) +#include "./hydro/MAGMA2/hydro_io.h" #elif defined(ANARCHY_PU_SPH) #include "./hydro/AnarchyPU/hydro_io.h" #else diff --git a/src/hydro_parameters.h b/src/hydro_parameters.h index 46d93ad43a..a370c7df71 100644 --- a/src/hydro_parameters.h +++ b/src/hydro_parameters.h @@ -58,6 +58,8 @@ #include "./hydro/SPHENIX/hydro_parameters.h" #elif defined(GASOLINE_SPH) #include "./hydro/Gasoline/hydro_parameters.h" +#elif defined(MAGMA2_SPH) +#include "./hydro/MAGMA2/hydro_parameters.h" #elif defined(ANARCHY_PU_SPH) #include "./hydro/AnarchyPU/hydro_parameters.h" #else diff --git a/src/hydro_properties.c b/src/hydro_properties.c index f8b547bc56..3107ed188a 100644 --- a/src/hydro_properties.c +++ b/src/hydro_properties.c @@ -63,19 +63,54 @@ void hydro_props_init(struct hydro_props *p, /* ------ Smoothing lengths parameters ---------- */ /* Kernel properties */ - p->eta_neighbours = parser_get_param_float(params, "SPH:resolution_eta"); + p->eta_neighbours = + parser_get_opt_param_float(params, "SPH:resolution_eta", 0.f); + + /* Target number of neighbours (optional) */ + p->target_neighbours = + parser_get_opt_param_float(params, "SPH:target_neighbours", 0.f); + + if (p->eta_neighbours <= 0.f && p->target_neighbours <=0.f) { + error("You must set either SPH:resolution_eta or SPH:target_neighbours " + "in the parameter file."); + } /* Tolerance for the smoothing length Newton-Raphson scheme */ p->h_tolerance = parser_get_opt_param_float(params, "SPH:h_tolerance", hydro_props_default_h_tolerance); /* Get derived properties */ - p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm; + /* Target number of neighbours */ + if (p->eta_neighbours > 0.f) { + p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm; + } + else { + p->eta_neighbours = + powf(p->target_neighbours / kernel_norm, hydro_dimension_inv); + } + const float delta_eta = p->eta_neighbours * (1.f + p->h_tolerance); p->delta_neighbours = (pow_dimension(delta_eta) - pow_dimension(p->eta_neighbours)) * kernel_norm; +#ifdef MAGMA2_SPH +#ifndef const_kernel_target_neighbours + error("When using MAGMA2 SPH, the constant " + "const_kernel_target_neighbours must be defined in the header file " + "hydro_parameters.h. This is a compile-time constant and must be set " + "to the desired number of neighbours in the parameter file."); +#else + if (fabsf((float)const_kernel_target_neighbours - p->target_neighbours) > + 0.05f * p->target_neighbours) { + error("When using MAGMA2 SPH, the compiled constant " + "const_kernel_target_neighbours (%g) must be within 5 percent of " + "the desired number of neighbours (%g) in the parameter file.", + (float)const_kernel_target_neighbours, p->target_neighbours); + } +#endif +#endif + #ifdef SHADOWFAX_SPH /* change the meaning of target_neighbours and delta_neighbours */ p->target_neighbours = 1.0f; diff --git a/src/part.h b/src/part.h index 08fb3121f8..9a8d70e95d 100644 --- a/src/part.h +++ b/src/part.h @@ -93,6 +93,10 @@ struct threadpool; #include "./hydro/Gasoline/hydro_part.h" #define hydro_need_extra_init_loop 0 #define EXTRA_HYDRO_LOOP +#elif defined(MAGMA2_SPH) +#include "./hydro/MAGMA2/hydro_part.h" +#define hydro_need_extra_init_loop 0 +#define EXTRA_HYDRO_LOOP #elif defined(ANARCHY_PU_SPH) #include "./hydro/AnarchyPU/hydro_part.h" #define hydro_need_extra_init_loop 0 diff --git a/src/sink/GEAR/sink.h b/src/sink/GEAR/sink.h index 47baafcffc..8f46702ecb 100644 --- a/src/sink/GEAR/sink.h +++ b/src/sink/GEAR/sink.h @@ -458,7 +458,7 @@ INLINE static int sink_is_forming( const float maximal_density_threshold = sink_props->maximal_density_threshold; const float density = hydro_get_physical_density(p, cosmo); - const float div_v = sink_get_physical_div_v_from_part(p); + const float div_v = sink_get_physical_div_v_from_part(p, cosmo); const float h = p->h; const float sink_cut_off_radius = sink_props->cut_off_radius; diff --git a/src/sink/GEAR/sink_getters.h b/src/sink/GEAR/sink_getters.h index 4612555b4a..beb75dda3a 100644 --- a/src/sink/GEAR/sink_getters.h +++ b/src/sink/GEAR/sink_getters.h @@ -82,7 +82,8 @@ INLINE static double sink_compute_neighbour_rotation_energy_magnitude( * */ INLINE static float sink_get_physical_div_v_from_part( - const struct part* restrict p) { + const struct part* restrict p, + const struct cosmology* cosmo) { float div_v = 0.0; @@ -106,6 +107,9 @@ INLINE static float sink_get_physical_div_v_from_part( div_v = (1. / 3.) * (p->viscosity.velocity_gradient[0][0] + p->viscosity.velocity_gradient[1][1] + p->viscosity.velocity_gradient[2][2]); +#elif MAGMA2_SPH + /* Copy the velocity divergence */ + div_v = hydro_get_physical_div_v(p, cosmo); #elif HOPKINS_PU_SPH div_v = p->density.div_v; #else diff --git a/src/star_formation/GEAR/star_formation.h b/src/star_formation/GEAR/star_formation.h index 40b95b2bd7..f0712925d6 100644 --- a/src/star_formation/GEAR/star_formation.h +++ b/src/star_formation/GEAR/star_formation.h @@ -448,6 +448,9 @@ __attribute__((always_inline)) INLINE static void star_formation_end_density( xp->sf_data.div_v = (1. / 3.) * (p->viscosity.velocity_gradient[0][0] + p->viscosity.velocity_gradient[1][1] + p->viscosity.velocity_gradient[2][2]); +#elif MAGMA2_SPH + /* Copy the velocity divergence */ + xp->sf_data.div_v = hydro_get_physical_div_v(p, cosmo); #elif HOPKINS_PU_SPH xp->sf_data.div_v = p->density.div_v; #else diff --git a/src/timestep.h b/src/timestep.h index a600241808..75e1ae759d 100644 --- a/src/timestep.h +++ b/src/timestep.h @@ -204,6 +204,27 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep( /* Limit timestep within the allowed range */ new_dt = min(new_dt, e->dt_max); +#ifdef MAGMA2_DEBUG_CHECKS + if (new_dt < e->dt_min) { + error("part (id=%lld) wants a time-step (%e) below dt_min (%e)" + "dt_h_change = %e, " + "new_dt_hydro = %e, new_dt_mhd = %e, " + "new_dt_cooling = %e, new_dt_grav = %e, " + "new_dt_forcing = %e, new_dt_chemistry = %e" + ", dt_max_RMS_displacement = %e", + p->id, + new_dt, e->dt_min, + dt_h_change * e->cosmology->time_step_factor, + new_dt_hydro * e->cosmology->time_step_factor, + new_dt_mhd * e->cosmology->time_step_factor, + new_dt_cooling * e->cosmology->time_step_factor, + new_dt_grav * e->cosmology->time_step_factor, + new_dt_forcing * e->cosmology->time_step_factor, + new_dt_chemistry * e->cosmology->time_step_factor, + e->dt_max_RMS_displacement * e->cosmology->time_step_factor); + } +#endif + if (new_dt < e->dt_min) error("part (id=%lld) wants a time-step (%e) below dt_min (%e)", p->id, new_dt, e->dt_min);