forces.hpp

Go to the documentation of this file.

// -------------------------------------------------------//
//
// SHAMROCK code for hydrodynamics
// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
//
// -------------------------------------------------------//
 
#pragma once

 
#include "shambase/constants.hpp"
#include "shambackends/math.hpp"
#include "shambackends/sycl.hpp"
#include "shammodels/gsph/math/reconstruction.hpp"
#include "shammodels/gsph/math/riemann/iterative.hpp"
#include "shammodels/sph/math/forces.hpp"
 
namespace shammodels::gsph {
 
    // Note: For GSPH acceleration, use shamrock::sph::sph_pressure_symetric() with p_star
    // as both P_a and P_b. This provides proper handling of zero denominators via
    // sham::inv_sat_zero() and avoids code duplication.

    template<class Tvec, class Tscal>
    inline Tscal gsph_energy_rate(
        Tscal m_b,
        Tscal p_star,
        Tscal v_star,
        Tscal rho_a_sq,
        Tscal rho_b_sq,
        Tscal omega_a,
        Tscal omega_b,
        Tvec v_a,
        Tvec r_ab_unit,
        Tvec nabla_W_a,
        Tvec nabla_W_b) {
 
        // Interface velocity vector (in direction of pair axis)
        Tvec v_star_vec = v_star * r_ab_unit;
 
        // Compute symmetric force (same as momentum equation)
        // f = m_b * p* * (nabla_W_a / (rho_a^2 * Omega_a) + nabla_W_b / (rho_b^2 * Omega_b))
        Tscal sub_fact_a = rho_a_sq * omega_a;
        Tscal sub_fact_b = rho_b_sq * omega_b;
        Tvec f           = m_b * p_star
                           * (nabla_W_a * sham::inv_sat_zero(sub_fact_a)
                              + nabla_W_b * sham::inv_sat_zero(sub_fact_b));
 
        // Energy rate: -f dot (v* - v_a)
        return -sycl::dot(f, v_star_vec - v_a);
    }

    template<class Tvec, class Tscal>
    inline void add_gsph_force_contribution(
        Tscal m_b,
        Tscal p_star,
        Tscal v_star,
        Tscal rho_a,
        Tscal rho_b,
        Tscal omega_a,
        Tscal omega_b,
        Tscal Fab_a,
        Tscal Fab_b,
        Tvec r_ab_unit,
        Tvec v_a,
        Tvec &dv_dt,
        Tscal &du_dt) {
 
        const Tscal rho_a_sq = rho_a * rho_a;
        const Tscal rho_b_sq = rho_b * rho_b;
 
        // Kernel gradient vectors (pointing from a to b)
        Tvec nabla_W_a = Fab_a * r_ab_unit;
        Tvec nabla_W_b = Fab_b * r_ab_unit;
 
        // Acceleration: use sph_pressure_symetric with p_star as both P_a and P_b
        // This provides proper handling of zero denominators via sham::inv_sat_zero()
        dv_dt += shamrock::sph::sph_pressure_symetric<Tvec, Tscal>(
            m_b, rho_a_sq, rho_b_sq, p_star, p_star, omega_a, omega_b, nabla_W_a, nabla_W_b);
 
        // Energy rate (uses symmetric force, same as momentum equation)
        du_dt += gsph_energy_rate<Tvec, Tscal>(
            m_b,
            p_star,
            v_star,
            rho_a_sq,
            rho_b_sq,
            omega_a,
            omega_b,
            v_a,
            r_ab_unit,
            nabla_W_a,
            nabla_W_b);
    }
 
    // Note: For velocity projection onto pair axis, use sycl::dot(v, r_ab_unit) directly.
    // For density from smoothing length, use shamrock::sph::rho_h() from density.hpp.

    template<class Tvec, class Tscal>
    inline void add_gsph_force_contribution_inutsuka(
        Tscal m_b,
        Tscal p_star,
        Tscal v_star,
        Tscal V2_ij,
        Tvec grad_W_ij,
        Tvec r_ab_unit,
        Tvec v_a,
        Tvec &dv_dt,
        Tscal &du_dt) {
 
        dv_dt -= m_b * p_star * V2_ij * grad_W_ij;
 
        Tvec v_star_vec = v_star * r_ab_unit;
        du_dt -= m_b * p_star * V2_ij * sycl::dot(grad_W_ij, v_star_vec - v_a);
    }

    template<class Kernel, class Tvec, class Tscal>
    inline void accumulate_gsph_pair_force(
        bool use_inutsuka_v2,
        Tscal pmass,
        Tscal p_star,
        Tscal v_star,
        Tscal rho_a,
        Tscal rho_b,
        Tscal omega_a,
        Tscal omega_b,
        Tscal rab,
        Tscal rab_inv,
        Tscal h_a,
        Tscal h_b,
        Tvec r_ab_unit,
        Tvec vxyz_a,
        Tvec &sum_axyz,
        Tscal &sum_du_a) {
 
        if (use_inutsuka_v2) {
            // Effective volume/face interpolation (Inutsuka 2002), linear
            // (1st order): specific volume is 1/rho for equal-mass particles.
            const Tscal vol_a = Tscal{1} / rho_a;
            const Tscal vol_b = Tscal{1} / rho_b;
 
            auto face = lin_v2_sast_ij<Tscal>(vol_a, vol_b, h_a, h_b, rab_inv);
 
            // Pair-symmetrized kernel gradient at sqrt(2)*h (Inutsuka 2002)
            constexpr Tscal sqrt2 = shambase::constants::sqrt_2<Tscal>;
            const Tscal Fab2_a    = Kernel::dW_3d(rab, sqrt2 * h_a);
            const Tscal Fab2_b    = Kernel::dW_3d(rab, sqrt2 * h_b);
            const Tvec grad_W_ij  = (Fab2_a + Fab2_b) * r_ab_unit;
 
            add_gsph_force_contribution_inutsuka<Tvec, Tscal>(
                pmass, p_star, v_star, face.V2, grad_W_ij, r_ab_unit, vxyz_a, sum_axyz, sum_du_a);
        } else {
            const Tscal Fab_a = Kernel::dW_3d(rab, h_a);
            const Tscal Fab_b = Kernel::dW_3d(rab, h_b);
 
            add_gsph_force_contribution<Tvec, Tscal>(
                pmass,
                p_star,
                v_star,
                rho_a,
                rho_b,
                omega_a,
                omega_b,
                Fab_a,
                Fab_b,
                r_ab_unit,
                vxyz_a,
                sum_axyz,
                sum_du_a);
        }
    }
 
} // namespace shammodels::gsph