/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2024 Darwin Roduit (darwin.roduit@alumni.epfl.ch)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
*******************************************************************************/
#ifndef SWIFT_GEAR_SINK_GETTERS_H
#define SWIFT_GEAR_SINK_GETTERS_H
#include "cosmology.h"
#include "sink_part.h"
/**
* @file src/sink/GEAR/sink_getters.h
* @brief Getter functions for GEAR sink scheme to avoid exposing
* implementation details to the outer world. Keep the code clean and lean.
*/
/**
* @brief Get the sink age in interal units.
*
* @param sink The #sink.
* @param with_cosmology If we run with cosmology.
* @param cosmo The co Birth scale-factor of the star.
* @param with_cosmology If we run with cosmology.
*/
__attribute__((always_inline)) INLINE double sink_get_sink_age(
const struct sink* restrict sink, const int with_cosmology,
const struct cosmology* cosmo, const double time) {
double sink_age;
if (with_cosmology) {
/* Deal with rounding issues */
if (sink->birth_scale_factor >= cosmo->a) {
sink_age = 0.;
} else {
sink_age = cosmology_get_delta_time_from_scale_factors(
cosmo, sink->birth_scale_factor, cosmo->a);
}
} else {
sink_age = time - sink->birth_time;
}
return sink_age;
}
/**
* @brief Compute the rotational energy of the neighbouring gas particles.
*
* Note: This function must be used after having computed the rotational energy
* per components, i.e. after sink_prepare_part_sink_formation().
*
* @param p The gas particle.
*
*/
INLINE static double sink_compute_neighbour_rotation_energy_magnitude(
const struct part* restrict p) {
double E_rot_x = p->sink_data.E_rot_neighbours[0];
double E_rot_y = p->sink_data.E_rot_neighbours[1];
double E_rot_z = p->sink_data.E_rot_neighbours[2];
double E_rot =
sqrtf(E_rot_x * E_rot_x + E_rot_y * E_rot_y + E_rot_z * E_rot_z);
return E_rot;
}
/**
* @brief Retrieve the physical velocity divergence from the gas particle.
*
* @param p The gas particles.
*
*/
INLINE static float sink_get_physical_div_v_from_part(
const struct part* restrict p) {
float div_v = 0.0;
/* The implementation of div_v depends on the Hydro scheme. Furthermore, some
add a Hubble flow term, some do not. We need to take care of this */
#ifdef SPHENIX_SPH
/* SPHENIX is already including the Hubble flow. */
div_v = hydro_get_div_v(p);
#elif GADGET2_SPH
div_v = p->density.div_v;
/* Add the missing term */
div_v += hydro_dimension * cosmo->H;
#elif MINIMAL_SPH
div_v = hydro_get_div_v(p);
/* Add the missing term */
div_v += hydro_dimension * cosmo->H;
#elif GASOLINE_SPH
/* Copy the velocity divergence */
div_v = (1. / 3.) * (p->viscosity.velocity_gradient[0][0] +
p->viscosity.velocity_gradient[1][1] +
p->viscosity.velocity_gradient[2][2]);
#elif HOPKINS_PU_SPH
div_v = p->density.div_v;
#else
#error \
"This scheme is not implemented. Note that Different scheme apply the Hubble flow in different places. Be careful about it."
#endif
return div_v;
}
/**
* @brief Compute the angular momentum-based criterion for sink-sink
* interaction.
*
* This function calculates the angular momentum of a sink particle relative to
* another particle (sink or gas) and evaluates the Keplerian angular momentum.
*
* @param dx Comoving vector separating the two particles (pi - pj).
* @param dv_plus_H_flow Comoving relative velocity including the Hubble flow.
* @param r Comoving distance between the two particles.
* @param r_cut_i Comoving cut-off radius of particle i.
* @param mass_i Mass of particle i.
* @param cosmo The cosmological parameters and properties
* @param grav_props The gravity scheme parameters and properties
* @param L2_kepler (return) Keplerian angular momentum squared of particle i.
* @param L2_j (return) Specific angular momentum squared relative to particle
* j.
*/
__attribute__((always_inline)) INLINE static void
sink_compute_angular_momenta_criterion(
const float dx[3], const float dv_plus_H_flow[3], const float r,
const float r_cut_i, const float mass_i, const struct cosmology* cosmo,
const struct gravity_props* grav_props, float* L2_kepler, float* L2_j) {
/* Compute the physical relative velocity between the particles */
const float dv_physical[3] = {dv_plus_H_flow[0] * cosmo->a_inv,
dv_plus_H_flow[1] * cosmo->a_inv,
dv_plus_H_flow[2] * cosmo->a_inv};
/* Compute the physical distance between the particles */
const float dx_physical[3] = {dx[0] * cosmo->a, dx[1] * cosmo->a,
dx[2] * cosmo->a};
const float r_physical = r * cosmo->a;
/* Momentum check------------------------------------------------------- */
/* Relative momentum of the gas */
const float specific_angular_momentum[3] = {
dx_physical[1] * dv_physical[2] - dx_physical[2] * dv_physical[1],
dx_physical[2] * dv_physical[0] - dx_physical[0] * dv_physical[2],
dx_physical[0] * dv_physical[1] - dx_physical[1] * dv_physical[0]};
*L2_j = specific_angular_momentum[0] * specific_angular_momentum[0] +
specific_angular_momentum[1] * specific_angular_momentum[1] +
specific_angular_momentum[2] * specific_angular_momentum[2];
/* Keplerian angular speed squared */
const float omega_acc_2 =
grav_props->G_Newton * mass_i / (r_physical * r_physical * r_physical);
/*Keplerian angular momentum squared */
*L2_kepler = (r_cut_i * r_cut_i * r_cut_i * r_cut_i) * omega_acc_2;
}
#endif /* SWIFT_GEAR_SINK_GETTERS_H */