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 "shambackends/math.hpp"
#include "shambackends/sycl.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.
 
} // namespace shammodels::gsph