Shamrock/doxygen/gsph_2src_2modules_2UpdateDerivs_8cpp_source.html

// -------------------------------------------------------//

//

// 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

//

// -------------------------------------------------------//


#include "shammodels/gsph/modules/UpdateDerivs.hpp"

#include "shambackends/math.hpp"

#include "shammath/sphkernels.hpp"

#include "shammodels/gsph/config/FieldNames.hpp"

#include "shammodels/gsph/math/forces.hpp"

#include "shammodels/gsph/math/reconstruction.hpp"

#include "shammodels/gsph/math/riemann/exact.hpp"

#include "shammodels/gsph/math/riemann/iterative.hpp"

#include "shammodels/sph/math/density.hpp"

#include "shamsys/legacy/log.hpp"

#include "shamtree/TreeTraversal.hpp"


template<class Tvec, template<class> class SPHKernel>


void shammodels::gsph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs() {

    StackEntry stack_loc{};


    Cfg_Riemann cfg_riemann = solver_config.riemann_config;


    if (Iterative *v = std::get_if<Iterative>(&cfg_riemann.config)) {

        update_derivs_iterative(*v);

    } else if (HLLC *v = std::get_if<HLLC>(&cfg_riemann.config)) {

        update_derivs_hllc(*v);

    } else if (Exact *v = std::get_if<Exact>(&cfg_riemann.config)) {

        update_derivs_exact(*v);

    } else {

        shambase::throw_unimplemented("Riemann solver type not supported by UpdateDerivs");

    }

}


template<class Tvec, template<class> class SPHKernel>

void shammodels::gsph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_iterative(

    Iterative cfg) {


    StackEntry stack_loc{};


    using namespace shamrock;

    using namespace shamrock::patch;


    PatchDataLayerLayout &pdl = scheduler().pdl_old();


    // Get field indices from the patch data layout

    const u32 ixyz   = pdl.get_field_idx<Tvec>(gsph::names::common::xyz);

    const u32 ivxyz  = pdl.get_field_idx<Tvec>(gsph::names::newtonian::vxyz);

    const u32 iaxyz  = pdl.get_field_idx<Tvec>(gsph::names::newtonian::axyz);

    const u32 ihpart = pdl.get_field_idx<Tscal>(gsph::names::common::hpart);


    // Optional internal energy fields (for adiabatic EOS)

    const bool has_uint = solver_config.has_field_uint();

    const u32 iuint     = has_uint ? pdl.get_field_idx<Tscal>(gsph::names::newtonian::uint) : 0;

    const u32 iduint    = has_uint ? pdl.get_field_idx<Tscal>(gsph::names::newtonian::duint) : 0;


    // Ghost layout for neighbor data

    shamrock::patch::PatchDataLayerLayout &ghost_layout

        = shambase::get_check_ref(storage.ghost_layout.get());

    u32 ihpart_interf   = ghost_layout.get_field_idx<Tscal>(gsph::names::common::hpart);

    u32 ivxyz_interf    = ghost_layout.get_field_idx<Tvec>(gsph::names::newtonian::vxyz);

    u32 iomega_interf   = ghost_layout.get_field_idx<Tscal>(gsph::names::newtonian::omega);

    u32 idensity_interf = ghost_layout.get_field_idx<Tscal>(gsph::names::newtonian::density);

    u32 iuint_interf

        = has_uint ? ghost_layout.get_field_idx<Tscal>(gsph::names::newtonian::uint) : 0;


    // Get merged data and caches from storage

    auto &merged_xyzh                                 = storage.merged_xyzh.get();

    shamrock::solvergraph::Field<Tscal> &omega_field  = shambase::get_check_ref(storage.omega);

    shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();


    // Get pressure and soundspeed from storage (includes ghosts)

    shamrock::solvergraph::Field<Tscal> &pressure_field = shambase::get_check_ref(storage.pressure);

    shamrock::solvergraph::Field<Tscal> &soundspeed_field

        = shambase::get_check_ref(storage.soundspeed);


    // Iterate over all non-empty patches

    scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {

        PatchDataLayer &mpdat = mpdats.get(cur_p.id_patch);


        // Get buffers for local and ghost data

        sham::DeviceBuffer<Tvec> &buf_xyz

            = merged_xyzh.get(cur_p.id_patch).template get_field_buf_ref<Tvec>(0);

        sham::DeviceBuffer<Tvec> &buf_axyz   = pdat.get_field_buf_ref<Tvec>(iaxyz);

        sham::DeviceBuffer<Tvec> &buf_vxyz   = mpdat.get_field_buf_ref<Tvec>(ivxyz_interf);

        sham::DeviceBuffer<Tscal> &buf_hpart = mpdat.get_field_buf_ref<Tscal>(ihpart_interf);

        sham::DeviceBuffer<Tscal> &buf_omega = mpdat.get_field_buf_ref<Tscal>(iomega_interf);

        sham::DeviceBuffer<Tscal> &buf_pressure

            = pressure_field.get_field(cur_p.id_patch).get_buf();

        sham::DeviceBuffer<Tscal> &buf_cs = soundspeed_field.get_field(cur_p.id_patch).get_buf();


        // Get neighbor cache for this patch

        tree::ObjectCache &pcache

            = shambase::get_check_ref(storage.neigh_cache).get_cache(cur_p.id_patch);


        // Set up SYCL queue and event tracking

        sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();

        sham::EventList depends_list;


        // Get density from merged ghost data (SPH summation density)

        sham::DeviceBuffer<Tscal> &buf_density = mpdat.get_field_buf_ref<Tscal>(idensity_interf);


        // Get buffer accessors

        auto xyz          = buf_xyz.get_read_access(depends_list);

        auto axyz         = buf_axyz.get_write_access(depends_list);

        auto vxyz         = buf_vxyz.get_read_access(depends_list);

        auto hpart        = buf_hpart.get_read_access(depends_list);

        auto omega_acc    = buf_omega.get_read_access(depends_list);

        auto density_acc  = buf_density.get_read_access(depends_list);

        auto pressure_acc = buf_pressure.get_read_access(depends_list);

        auto cs_acc       = buf_cs.get_read_access(depends_list);

        auto ploop_ptrs   = pcache.get_read_access(depends_list);


        // Optional: internal energy

        sham::DeviceBuffer<Tscal> *buf_duint_ptr = nullptr;

        Tscal *duint_acc                         = nullptr;

        if (has_uint) {

            buf_duint_ptr = &pdat.get_field_buf_ref<Tscal>(iduint);

            duint_acc     = buf_duint_ptr->get_write_access(depends_list);

        }


        auto e = q.submit(depends_list, [&](sycl::handler &cgh) {

            const Tscal pmass          = solver_config.gpart_mass;

            const Tscal gamma          = solver_config.get_eos_gamma();

            const Tscal tol            = cfg.tol;

            const u32 max_iter         = cfg.max_iter;

            const bool do_energy       = has_uint;

            const bool use_inutsuka_v2 = solver_config.is_force_inutsuka_v2();


            // Use shamrock's ObjectCacheIterator for neighbor traversal

            tree::ObjectCacheIterator particle_looper(ploop_ptrs);


            constexpr Tscal Rker2 = Kernel::Rkern * Kernel::Rkern;


            // InutsukaV2 evaluates the kernel gradient at sqrt(2)*h, so the local

            // kernel-support filter below must match the widened cache radius

            // (see start_neighbors_cache()/compute_presteps_rint() in Solver.cpp),

            // or valid cached pairs beyond h*Rkern get dropped here regardless.

            const Tscal Rker2_factor = use_inutsuka_v2 ? Tscal{2} : Tscal{1};


            shambase::parallel_for(cgh, pdat.get_obj_cnt(), "GSPH derivs iterative", [=](u64 gid) {

                u32 id_a = (u32) gid;


                using namespace shamrock::sph;


                // Initialize accumulators

                Tvec sum_axyz  = {0, 0, 0};

                Tscal sum_du_a = 0;


                // Particle a state

                const Tscal h_a     = hpart[id_a];

                const Tvec xyz_a    = xyz[id_a];

                const Tvec vxyz_a   = vxyz[id_a];

                const Tscal omega_a = omega_acc[id_a];


                // Use SPH-summation density

                const Tscal rho_a = sycl::max(density_acc[id_a], Tscal(1e-30));


                // Use pressure and soundspeed from storage

                const Tscal P_a  = sycl::max(pressure_acc[id_a], Tscal(1e-30));

                const Tscal cs_a = sycl::max(cs_acc[id_a], Tscal(1e-10));


                // Loop over neighbors

                particle_looper.for_each_object(id_a, [&](u32 id_b) {

                    if (id_a == id_b)

                        return; // Skip self


                    // Distance and kernel support check

                    const Tvec dr    = xyz_a - xyz[id_b];

                    const Tscal rab2 = sycl::dot(dr, dr);

                    const Tscal h_b  = hpart[id_b];


                    // Skip if outside kernel support

                    if (rab2 > h_a * h_a * Rker2 * Rker2_factor

                        && rab2 > h_b * h_b * Rker2 * Rker2_factor) {

                        return;

                    }


                    const Tscal rab     = sycl::sqrt(rab2);

                    const Tvec vxyz_b   = vxyz[id_b];

                    const Tscal omega_b = omega_acc[id_b];


                    // Use SPH-summation density

                    const Tscal rho_b = sycl::max(density_acc[id_b], Tscal(1e-30));


                    // Use pressure and soundspeed from storage

                    const Tscal P_b  = sycl::max(pressure_acc[id_b], Tscal(1e-30));

                    const Tscal cs_b = sycl::max(cs_acc[id_b], Tscal(1e-10));


                    // Unit vector from a to b

                    const Tscal rab_inv  = sham::inv_sat_positive(rab);

                    const Tvec r_ab_unit = dr * rab_inv;


                    // Project velocities onto pair axis for 1D Riemann problem

                    const Tscal u_a_proj = sycl::dot(vxyz_a, r_ab_unit);

                    const Tscal u_b_proj = sycl::dot(vxyz_b, r_ab_unit);


                    // Solve 1D Riemann problem using iterative solver

                    // Convention: Left state = neighbor b, Right state = current a

                    auto riemann_result = riemann::iterative_solver<Tscal>(

                        u_b_proj,

                        rho_b,

                        P_b, // Left = neighbor

                        u_a_proj,

                        rho_a,

                        P_a, // Right = current

                        gamma,

                        tol,

                        max_iter);

                    const Tscal p_star = riemann_result.p_star;

                    const Tscal v_star = riemann_result.v_star;


                    shammodels::gsph::accumulate_gsph_pair_force<Kernel, Tvec, Tscal>(

                        use_inutsuka_v2,

                        pmass,

                        p_star,

                        v_star,

                        rho_a,

                        rho_b,

                        omega_a,

                        omega_b,

                        rab,

                        rab_inv,

                        h_a,

                        h_b,

                        r_ab_unit,

                        vxyz_a,

                        sum_axyz,

                        sum_du_a);

                });


                // Write accumulated derivatives

                axyz[id_a] = sum_axyz;

                if (duint_acc != nullptr) {

                    duint_acc[id_a] = sum_du_a;

                }

            });

        });


        // Complete event states for all buffers

        buf_xyz.complete_event_state(e);

        buf_axyz.complete_event_state(e);

        buf_vxyz.complete_event_state(e);

        buf_hpart.complete_event_state(e);

        buf_omega.complete_event_state(e);

        buf_density.complete_event_state(e);

        buf_pressure.complete_event_state(e);

        buf_cs.complete_event_state(e);


        if (has_uint && buf_duint_ptr) {

            buf_duint_ptr->complete_event_state(e);

        }


        sham::EventList resulting_events;

        resulting_events.add_event(e);

        pcache.complete_event_state(resulting_events);

    });

}


template<class Tvec, template<class> class SPHKernel>

void shammodels::gsph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_exact(Exact cfg) {


    StackEntry stack_loc{};


    using namespace shamrock;

    using namespace shamrock::patch;


    PatchDataLayerLayout &pdl = scheduler().pdl_old();


    // Get field indices from the patch data layout

    const u32 ixyz   = pdl.get_field_idx<Tvec>(gsph::names::common::xyz);

    const u32 ivxyz  = pdl.get_field_idx<Tvec>(gsph::names::newtonian::vxyz);

    const u32 iaxyz  = pdl.get_field_idx<Tvec>(gsph::names::newtonian::axyz);

    const u32 ihpart = pdl.get_field_idx<Tscal>(gsph::names::common::hpart);


    // Optional internal energy fields (for adiabatic EOS)

    const bool has_uint = solver_config.has_field_uint();

    const u32 iuint     = has_uint ? pdl.get_field_idx<Tscal>(gsph::names::newtonian::uint) : 0;

    const u32 iduint    = has_uint ? pdl.get_field_idx<Tscal>(gsph::names::newtonian::duint) : 0;


    // Ghost layout for neighbor data

    shamrock::patch::PatchDataLayerLayout &ghost_layout

        = shambase::get_check_ref(storage.ghost_layout.get());

    u32 ihpart_interf   = ghost_layout.get_field_idx<Tscal>(gsph::names::common::hpart);

    u32 ivxyz_interf    = ghost_layout.get_field_idx<Tvec>(gsph::names::newtonian::vxyz);

    u32 iomega_interf   = ghost_layout.get_field_idx<Tscal>(gsph::names::newtonian::omega);

    u32 idensity_interf = ghost_layout.get_field_idx<Tscal>(gsph::names::newtonian::density);

    u32 iuint_interf

        = has_uint ? ghost_layout.get_field_idx<Tscal>(gsph::names::newtonian::uint) : 0;


    // Get merged data and caches from storage

    auto &merged_xyzh                                 = storage.merged_xyzh.get();

    shamrock::solvergraph::Field<Tscal> &omega_field  = shambase::get_check_ref(storage.omega);

    shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();


    // Get pressure and soundspeed from storage (includes ghosts)

    shamrock::solvergraph::Field<Tscal> &pressure_field = shambase::get_check_ref(storage.pressure);

    shamrock::solvergraph::Field<Tscal> &soundspeed_field

        = shambase::get_check_ref(storage.soundspeed);


    // Iterate over all non-empty patches

    scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {

        PatchDataLayer &mpdat = mpdats.get(cur_p.id_patch);


        // Get buffers for local and ghost data

        sham::DeviceBuffer<Tvec> &buf_xyz

            = merged_xyzh.get(cur_p.id_patch).template get_field_buf_ref<Tvec>(0);

        sham::DeviceBuffer<Tvec> &buf_axyz   = pdat.get_field_buf_ref<Tvec>(iaxyz);

        sham::DeviceBuffer<Tvec> &buf_vxyz   = mpdat.get_field_buf_ref<Tvec>(ivxyz_interf);

        sham::DeviceBuffer<Tscal> &buf_hpart = mpdat.get_field_buf_ref<Tscal>(ihpart_interf);

        sham::DeviceBuffer<Tscal> &buf_omega = mpdat.get_field_buf_ref<Tscal>(iomega_interf);

        sham::DeviceBuffer<Tscal> &buf_pressure

            = pressure_field.get_field(cur_p.id_patch).get_buf();

        sham::DeviceBuffer<Tscal> &buf_cs = soundspeed_field.get_field(cur_p.id_patch).get_buf();


        // Get neighbor cache for this patch

        tree::ObjectCache &pcache

            = shambase::get_check_ref(storage.neigh_cache).get_cache(cur_p.id_patch);


        // Set up SYCL queue and event tracking

        sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();

        sham::EventList depends_list;


        // Get density from merged ghost data (SPH summation density)

        sham::DeviceBuffer<Tscal> &buf_density = mpdat.get_field_buf_ref<Tscal>(idensity_interf);


        // Get buffer accessors

        auto xyz          = buf_xyz.get_read_access(depends_list);

        auto axyz         = buf_axyz.get_write_access(depends_list);

        auto vxyz         = buf_vxyz.get_read_access(depends_list);

        auto hpart        = buf_hpart.get_read_access(depends_list);

        auto omega_acc    = buf_omega.get_read_access(depends_list);

        auto density_acc  = buf_density.get_read_access(depends_list);

        auto pressure_acc = buf_pressure.get_read_access(depends_list);

        auto cs_acc       = buf_cs.get_read_access(depends_list);

        auto ploop_ptrs   = pcache.get_read_access(depends_list);


        // Optional: internal energy

        sham::DeviceBuffer<Tscal> *buf_duint_ptr = nullptr;

        Tscal *duint_acc                         = nullptr;

        if (has_uint) {

            buf_duint_ptr = &pdat.get_field_buf_ref<Tscal>(iduint);

            duint_acc     = buf_duint_ptr->get_write_access(depends_list);

        }


        auto e = q.submit(depends_list, [&](sycl::handler &cgh) {

            const Tscal pmass          = solver_config.gpart_mass;

            const Tscal gamma          = solver_config.get_eos_gamma();

            const Tscal tol            = cfg.tol;

            const u32 max_iter         = cfg.max_iter;

            const bool do_energy       = has_uint;

            const bool use_inutsuka_v2 = solver_config.is_force_inutsuka_v2();


            // Use shamrock's ObjectCacheIterator for neighbor traversal

            tree::ObjectCacheIterator particle_looper(ploop_ptrs);


            constexpr Tscal Rker2 = Kernel::Rkern * Kernel::Rkern;


            // InutsukaV2 evaluates the kernel gradient at sqrt(2)*h, so the local

            // kernel-support filter below must match the widened cache radius

            // (see start_neighbors_cache()/compute_presteps_rint() in Solver.cpp),

            // or valid cached pairs beyond h*Rkern get dropped here regardless.

            const Tscal Rker2_factor = use_inutsuka_v2 ? Tscal{2} : Tscal{1};


            shambase::parallel_for(cgh, pdat.get_obj_cnt(), "GSPH derivs exact", [=](u64 gid) {

                u32 id_a = (u32) gid;


                using namespace shamrock::sph;


                // Initialize accumulators

                Tvec sum_axyz  = {0, 0, 0};

                Tscal sum_du_a = 0;


                // Particle a state

                const Tscal h_a     = hpart[id_a];

                const Tvec xyz_a    = xyz[id_a];

                const Tvec vxyz_a   = vxyz[id_a];

                const Tscal omega_a = omega_acc[id_a];


                // Use SPH-summation density

                const Tscal rho_a = sycl::max(density_acc[id_a], Tscal(1e-30));


                // Use pressure and soundspeed from storage

                const Tscal P_a  = sycl::max(pressure_acc[id_a], Tscal(1e-30));

                const Tscal cs_a = sycl::max(cs_acc[id_a], Tscal(1e-10));


                // Loop over neighbors

                particle_looper.for_each_object(id_a, [&](u32 id_b) {

                    if (id_a == id_b)

                        return; // Skip self


                    // Distance and kernel support check

                    const Tvec dr    = xyz_a - xyz[id_b];

                    const Tscal rab2 = sycl::dot(dr, dr);

                    const Tscal h_b  = hpart[id_b];


                    // Skip if outside kernel support

                    if (rab2 > h_a * h_a * Rker2 * Rker2_factor

                        && rab2 > h_b * h_b * Rker2 * Rker2_factor) {

                        return;

                    }


                    const Tscal rab     = sycl::sqrt(rab2);

                    const Tvec vxyz_b   = vxyz[id_b];

                    const Tscal omega_b = omega_acc[id_b];


                    // Use SPH-summation density

                    const Tscal rho_b = sycl::max(density_acc[id_b], Tscal(1e-30));


                    // Use pressure and soundspeed from storage

                    const Tscal P_b  = sycl::max(pressure_acc[id_b], Tscal(1e-30));

                    const Tscal cs_b = sycl::max(cs_acc[id_b], Tscal(1e-10));


                    // Unit vector from a to b

                    const Tscal rab_inv  = sham::inv_sat_positive(rab);

                    const Tvec r_ab_unit = dr * rab_inv;


                    // Project velocities onto pair axis for 1D Riemann problem

                    const Tscal u_a_proj = sycl::dot(vxyz_a, r_ab_unit);

                    const Tscal u_b_proj = sycl::dot(vxyz_b, r_ab_unit);


                    // Solve 1D Riemann problem exactly (Toro)

                    // Convention: Left state = neighbor b, Right state = current a

                    auto riemann_result = riemann::exact_solver<Tscal>(

                        u_b_proj,

                        rho_b,

                        P_b, // Left = neighbor

                        u_a_proj,

                        rho_a,

                        P_a, // Right = current

                        gamma,

                        tol,

                        max_iter);

                    const Tscal p_star = riemann_result.p_star;

                    const Tscal v_star = riemann_result.v_star;


                    shammodels::gsph::accumulate_gsph_pair_force<Kernel, Tvec, Tscal>(

                        use_inutsuka_v2,

                        pmass,

                        p_star,

                        v_star,

                        rho_a,

                        rho_b,

                        omega_a,

                        omega_b,

                        rab,

                        rab_inv,

                        h_a,

                        h_b,

                        r_ab_unit,

                        vxyz_a,

                        sum_axyz,

                        sum_du_a);

                });


                // Write accumulated derivatives

                axyz[id_a] = sum_axyz;

                if (duint_acc != nullptr) {

                    duint_acc[id_a] = sum_du_a;

                }

            });

        });


        // Complete event states for all buffers

        buf_xyz.complete_event_state(e);

        buf_axyz.complete_event_state(e);

        buf_vxyz.complete_event_state(e);

        buf_hpart.complete_event_state(e);

        buf_omega.complete_event_state(e);

        buf_density.complete_event_state(e);

        buf_pressure.complete_event_state(e);

        buf_cs.complete_event_state(e);


        if (has_uint && buf_duint_ptr) {

            buf_duint_ptr->complete_event_state(e);

        }


        sham::EventList resulting_events;

        resulting_events.add_event(e);

        pcache.complete_event_state(resulting_events);

    });

}


template<class Tvec, template<class> class SPHKernel>

void shammodels::gsph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_hllc(HLLC cfg) {


    StackEntry stack_loc{};


    using namespace shamrock;

    using namespace shamrock::patch;


    PatchDataLayerLayout &pdl = scheduler().pdl_old();


    // Get field indices

    const u32 ixyz   = pdl.get_field_idx<Tvec>(gsph::names::common::xyz);

    const u32 ivxyz  = pdl.get_field_idx<Tvec>(gsph::names::newtonian::vxyz);

    const u32 iaxyz  = pdl.get_field_idx<Tvec>(gsph::names::newtonian::axyz);

    const u32 ihpart = pdl.get_field_idx<Tscal>(gsph::names::common::hpart);


    const bool has_uint = solver_config.has_field_uint();

    const u32 iuint     = has_uint ? pdl.get_field_idx<Tscal>(gsph::names::newtonian::uint) : 0;

    const u32 iduint    = has_uint ? pdl.get_field_idx<Tscal>(gsph::names::newtonian::duint) : 0;


    // Ghost layout

    shamrock::patch::PatchDataLayerLayout &ghost_layout

        = shambase::get_check_ref(storage.ghost_layout.get());

    u32 ihpart_interf   = ghost_layout.get_field_idx<Tscal>(gsph::names::common::hpart);

    u32 ivxyz_interf    = ghost_layout.get_field_idx<Tvec>(gsph::names::newtonian::vxyz);

    u32 iomega_interf   = ghost_layout.get_field_idx<Tscal>(gsph::names::newtonian::omega);

    u32 idensity_interf = ghost_layout.get_field_idx<Tscal>(gsph::names::newtonian::density);

    u32 iuint_interf

        = has_uint ? ghost_layout.get_field_idx<Tscal>(gsph::names::newtonian::uint) : 0;


    auto &merged_xyzh                                 = storage.merged_xyzh.get();

    shamrock::solvergraph::Field<Tscal> &omega_field  = shambase::get_check_ref(storage.omega);

    shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();


    // Get pressure and soundspeed from storage (includes ghosts)

    shamrock::solvergraph::Field<Tscal> &pressure_field = shambase::get_check_ref(storage.pressure);

    shamrock::solvergraph::Field<Tscal> &soundspeed_field

        = shambase::get_check_ref(storage.soundspeed);


    scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {

        PatchDataLayer &mpdat = mpdats.get(cur_p.id_patch);


        sham::DeviceBuffer<Tvec> &buf_xyz

            = merged_xyzh.get(cur_p.id_patch).template get_field_buf_ref<Tvec>(0);

        sham::DeviceBuffer<Tvec> &buf_axyz   = pdat.get_field_buf_ref<Tvec>(iaxyz);

        sham::DeviceBuffer<Tvec> &buf_vxyz   = mpdat.get_field_buf_ref<Tvec>(ivxyz_interf);

        sham::DeviceBuffer<Tscal> &buf_hpart = mpdat.get_field_buf_ref<Tscal>(ihpart_interf);

        sham::DeviceBuffer<Tscal> &buf_omega = mpdat.get_field_buf_ref<Tscal>(iomega_interf);

        sham::DeviceBuffer<Tscal> &buf_pressure

            = pressure_field.get_field(cur_p.id_patch).get_buf();

        sham::DeviceBuffer<Tscal> &buf_cs = soundspeed_field.get_field(cur_p.id_patch).get_buf();


        tree::ObjectCache &pcache

            = shambase::get_check_ref(storage.neigh_cache).get_cache(cur_p.id_patch);


        sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();

        sham::EventList depends_list;


        // Get density from merged ghost data

        sham::DeviceBuffer<Tscal> &buf_density = mpdat.get_field_buf_ref<Tscal>(idensity_interf);


        auto xyz          = buf_xyz.get_read_access(depends_list);

        auto axyz         = buf_axyz.get_write_access(depends_list);

        auto vxyz         = buf_vxyz.get_read_access(depends_list);

        auto hpart        = buf_hpart.get_read_access(depends_list);

        auto omega_acc    = buf_omega.get_read_access(depends_list);

        auto density_acc  = buf_density.get_read_access(depends_list);

        auto pressure_acc = buf_pressure.get_read_access(depends_list);

        auto cs_acc       = buf_cs.get_read_access(depends_list);

        auto ploop_ptrs   = pcache.get_read_access(depends_list);


        sham::DeviceBuffer<Tscal> *buf_duint_ptr = nullptr;

        Tscal *duint_acc                         = nullptr;

        if (has_uint) {

            buf_duint_ptr = &pdat.get_field_buf_ref<Tscal>(iduint);

            duint_acc     = buf_duint_ptr->get_write_access(depends_list);

        }


        auto e = q.submit(depends_list, [&](sycl::handler &cgh) {

            const Tscal pmass    = solver_config.gpart_mass;

            const Tscal gamma    = solver_config.get_eos_gamma();

            const bool do_energy = has_uint;


            tree::ObjectCacheIterator particle_looper(ploop_ptrs);


            constexpr Tscal Rker2 = Kernel::Rkern * Kernel::Rkern;


            shambase::parallel_for(cgh, pdat.get_obj_cnt(), "GSPH derivs HLLC", [=](u64 gid) {

                u32 id_a = (u32) gid;


                using namespace shamrock::sph;


                Tvec sum_axyz  = {0, 0, 0};

                Tscal sum_du_a = 0;


                const Tscal h_a     = hpart[id_a];

                const Tvec xyz_a    = xyz[id_a];

                const Tvec vxyz_a   = vxyz[id_a];

                const Tscal omega_a = omega_acc[id_a];


                // Use SPH-summation density

                const Tscal rho_a = sycl::max(density_acc[id_a], Tscal(1e-30));


                // Use pressure and soundspeed from storage

                const Tscal P_a  = sycl::max(pressure_acc[id_a], Tscal(1e-30));

                const Tscal cs_a = sycl::max(cs_acc[id_a], Tscal(1e-10));


                particle_looper.for_each_object(id_a, [&](u32 id_b) {

                    if (id_a == id_b)

                        return;


                    const Tvec dr    = xyz_a - xyz[id_b];

                    const Tscal rab2 = sycl::dot(dr, dr);

                    const Tscal h_b  = hpart[id_b];


                    if (rab2 > h_a * h_a * Rker2 && rab2 > h_b * h_b * Rker2) {

                        return;

                    }


                    const Tscal rab     = sycl::sqrt(rab2);

                    const Tvec vxyz_b   = vxyz[id_b];

                    const Tscal omega_b = omega_acc[id_b];


                    // Use SPH-summation density

                    const Tscal rho_b = sycl::max(density_acc[id_b], Tscal(1e-30));


                    // Use pressure and soundspeed from storage

                    const Tscal P_b  = sycl::max(pressure_acc[id_b], Tscal(1e-30));

                    const Tscal cs_b = sycl::max(cs_acc[id_b], Tscal(1e-10));


                    const Tscal rab_inv  = sham::inv_sat_positive(rab);

                    const Tvec r_ab_unit = dr * rab_inv;


                    // Project velocities onto pair axis for 1D Riemann problem

                    const Tscal u_a_proj = sycl::dot(vxyz_a, r_ab_unit);

                    const Tscal u_b_proj = sycl::dot(vxyz_b, r_ab_unit);


                    // Use HLLC approximate Riemann solver

                    // Convention: Left state = neighbor b, Right state = current a

                    auto riemann_result = riemann::hllc_solver<Tscal>(

                        u_b_proj,

                        rho_b,

                        P_b, // Left = neighbor

                        u_a_proj,

                        rho_a,

                        P_a, // Right = current

                        gamma);

                    const Tscal p_star = riemann_result.p_star;

                    const Tscal v_star = riemann_result.v_star;


                    const Tscal Fab_a = Kernel::dW_3d(rab, h_a);

                    const Tscal Fab_b = Kernel::dW_3d(rab, h_b);


                    // GSPH force contribution

                    shammodels::gsph::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);

                });


                axyz[id_a] = sum_axyz;

                if (duint_acc != nullptr) {

                    duint_acc[id_a] = sum_du_a;

                }

            });

        });


        buf_xyz.complete_event_state(e);

        buf_axyz.complete_event_state(e);

        buf_vxyz.complete_event_state(e);

        buf_hpart.complete_event_state(e);

        buf_omega.complete_event_state(e);

        buf_density.complete_event_state(e);

        buf_pressure.complete_event_state(e);

        buf_cs.complete_event_state(e);


        if (has_uint && buf_duint_ptr) {

            buf_duint_ptr->complete_event_state(e);

        }


        sham::EventList resulting_events;

        resulting_events.add_event(e);

        pcache.complete_event_state(resulting_events);

    });

}


// Explicit template instantiations

// M-spline kernels (Monaghan)

using namespace shammath;

template class shammodels::gsph::modules::UpdateDerivs<f64_3, M4>;

template class shammodels::gsph::modules::UpdateDerivs<f64_3, M6>;

template class shammodels::gsph::modules::UpdateDerivs<f64_3, M8>;


// Wendland kernels (C2, C4, C6)

template class shammodels::gsph::modules::UpdateDerivs<f64_3, C2>;

template class shammodels::gsph::modules::UpdateDerivs<f64_3, C4>;

template class shammodels::gsph::modules::UpdateDerivs<f64_3, C6>;

FieldNames.hpp
Constants for field names in GSPH solver, organized by physics mode.

shammodels::gsph::names::newtonian::duint
constexpr const char * duint
Time derivative of internal energy du/dt.
Definition FieldNames.hpp:52

shammodels::gsph::names::newtonian::axyz
constexpr const char * axyz
3-acceleration field
Definition FieldNames.hpp:46

shammodels::gsph::names::newtonian::uint
constexpr const char * uint
Specific internal energy u.
Definition FieldNames.hpp:49

shammodels::gsph::names::newtonian::vxyz
constexpr const char * vxyz
3-velocity field
Definition FieldNames.hpp:43

shammodels::gsph::names::common::xyz
constexpr const char * xyz
Position field (3D coordinates).
Definition FieldNames.hpp:33

shammodels::gsph::names::newtonian::density
constexpr const char * density
Density \rho (derived from h).
Definition FieldNames.hpp:55

shammodels::gsph::names::common::hpart
constexpr const char * hpart
Smoothing length field.
Definition FieldNames.hpp:34

shammodels::gsph::names::newtonian::omega
constexpr const char * omega
Grad-h correction factor \Omega.
Definition FieldNames.hpp:64

TreeTraversal.hpp

u32
std::uint32_t u32
32 bit unsigned integer
Definition aliases_int.hpp:27

u64
std::uint64_t u64
64 bit unsigned integer
Definition aliases_int.hpp:26

sham::DeviceBuffer
A buffer allocated in USM (Unified Shared Memory).
Definition DeviceBuffer.hpp:106

sham::DeviceBuffer::complete_event_state
void complete_event_state(sycl::event e) const
Complete the event state of the buffer.
Definition DeviceBuffer.hpp:368

sham::DeviceBuffer::get_write_access
T * get_write_access(sham::EventList &depends_list, SourceLocation src_loc=SourceLocation{})
Get a read-write pointer to the buffer's data.
Definition DeviceBuffer.hpp:349

sham::DeviceBuffer::get_read_access
const T * get_read_access(sham::EventList &depends_list, SourceLocation src_loc=SourceLocation{}) const
Get a read-only pointer to the buffer's data.
Definition DeviceBuffer.hpp:331

sham::DeviceQueue
A SYCL queue associated with a device and a context.
Definition DeviceQueue.hpp:32

sham::DeviceQueue::submit
sycl::event submit(Fct &&fct)
Submits a kernel to the SYCL queue.
Definition DeviceQueue.hpp:101

sham::EventList
Class to manage a list of SYCL events.
Definition EventList.hpp:31

sham::EventList::add_event
void add_event(sycl::event e)
Add an event to the list of events.
Definition EventList.hpp:87

shambase::DistributedData
Represents a collection of objects distributed across patches identified by a u64 id.
Definition DistributedData.hpp:43

shambase::DistributedData::get
T & get(u64 id)
Returns a reference to an object in the collection.
Definition DistributedData.hpp:145

shammodels::gsph::modules::UpdateDerivs
GSPH derivative update module.
Definition UpdateDerivs.hpp:52

shammodels::gsph::modules::UpdateDerivs::update_derivs
void update_derivs()
Update all derivatives using GSPH Riemann solver approach.
Definition UpdateDerivs.cpp:42

shamrock::patch::PatchDataLayerLayout
Definition PatchDataLayerLayout.hpp:66

shamrock::patch::PatchDataLayerLayout::get_field_idx
u32 get_field_idx(const std::string &field_name) const
Get the field id if matching name & type.
Definition PatchDataLayerLayout.hpp:349

shamrock::patch::PatchDataLayer
PatchDataLayer container class, the layout is described in patchdata_layout.
Definition PatchDataLayer.hpp:35

shamrock::solvergraph::Field
Definition Field.hpp:27

shamrock::solvergraph::IFieldRefs::get_field
PatchDataField< T > & get_field(u64 id) const
Get the underlying PatchDataField at the given id.
Definition IFieldRefs.hpp:54

shamrock.tree::ObjectCacheIterator
Definition TreeTraversal.hpp:487

density.hpp

exact.hpp
Exact Riemann solver for GSPH (Toro 2009).

forces.hpp
GSPH force computation using Riemann solver results.

shammodels::gsph::accumulate_gsph_pair_force
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)
Dispatch a single neighbor pair's force contribution to ChaWhitworth or InutsukaV2,...
Definition forces.hpp:232

shammodels::gsph::add_gsph_force_contribution
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)
Add GSPH force contribution from a single neighbor pair.
Definition forces.hpp:121

UpdateDerivs.hpp
GSPH derivative update module.

iterative.hpp
Iterative Riemann solver for GSPH (van Leer 1997).

log.hpp

math.hpp

sham::inv_sat_positive
T inv_sat_positive(T v, T minvsat=T{1e-9}, T satval=T{0.}) noexcept
inverse saturated (positive numbers only)
Definition math.hpp:841

shambase::get_check_ref
T & get_check_ref(const std::unique_ptr< T > &ptr, SourceLocation loc=SourceLocation())
Takes a std::unique_ptr and returns a reference to the object it holds. It throws a std::runtime_erro...
Definition memory.hpp:110

shambase::throw_unimplemented
void throw_unimplemented(SourceLocation loc=SourceLocation{})
Throw a std::runtime_error saying that the function is unimplemented.
Definition exception.hpp:153

shammath
namespace for math utility
Definition AABB.hpp:26

shamrock
namespace for the main framework
Definition __init__.py:1

reconstruction.hpp
Effective face (volume element) interpolation for the Inutsuka (2002) GSPH formulation.

sphkernels.hpp
sph kernels

StackEntry
shambase::details::BasicStackEntry StackEntry
Alias for shambase::details::BasicStackEntry.
Definition stacktrace.hpp:242

shamrock::patch::Patch
Patch object that contain generic patch information.
Definition Patch.hpp:33

shamrock::patch::Patch::id_patch
u64 id_patch
unique key that identify the patch
Definition Patch.hpp:86

shamrock.tree::ObjectCache
Definition TreeTraversal.hpp:375