api_device_sph.hpp | Ansys Developer Portal

#pragma once
 
#include <rocky20/api/rocky_api_sph_functions.h>
#include <rocky20/api/rocky_model_api.hpp>
#include <rocky20/api/sph/rocky_sph_element_api.hpp>
#include <rocky20/api/sph/rocky_sph_interaction_api.hpp>
#include <rocky20/api/sph/rocky_sph_model_api.hpp>
 
namespace device {
 
template <class Backend>
struct DefaultElementsOperator
{
    template <typename Callable>
    __host__ __device__ void operator()(
        Callable function_ptr,
        IRockyDeviceModel& device_model_api,
        IRockySPHDeviceModel& sph_model,
        IRockySPHElementInteraction& sph_interaction_api,
        void* udf_data)
    {
        function_ptr(device_model_api, sph_model, sph_interaction_api, udf_data);
    }
};
 
template <class Backend>
struct DefaultTrianglesOperator
{
    template <typename Callable>
    __host__ __device__ void operator()(
        Callable function_ptr,
        IRockyDeviceModel& device_model_api,
        IRockySPHDeviceModel& sph_model,
        IRockySPHTriangleInteraction& sph_interaction_api,
        void* udf_data)
    {
        function_ptr(device_model_api, sph_model, sph_interaction_api, udf_data);
    }
};
 
template <class Backend>
struct DefaultOperator
{
    template <typename Callable>
    __host__ __device__ void operator()(
        Callable function_ptr,
        IRockyDeviceModel& device_model_api,
        IRockySPHDeviceModel& sph_model,
        IRockySPHElement& sph_element_api,
        void* udf_data)
    {
        function_ptr(device_model_api, sph_model, sph_element_api, udf_data);
    }
};
 
template <class Backend>
struct HeatTransferElementsOperator
{
    template <typename Callable>
    __host__ __device__ void operator()(
        Callable function_ptr,
        IRockyDeviceModel& device_model_api,
        IRockySPHDeviceModel& sph_model,
        IRockySPHElementInteraction& sph_interaction_api,
        void* udf_data)
    {
        // Compute and add heat transfer
        real_sph heat_transfer =
            function_ptr(device_model_api, sph_model, sph_interaction_api, udf_data);
        auto home_sph_element = sph_interaction_api.get_home_element();
        auto near_sph_element = sph_interaction_api.get_near_element();
        home_sph_element.get_scalars().add_heat_transfer(heat_transfer);
        near_sph_element.get_scalars().add_heat_transfer(-heat_transfer);
 
        // Recording coupled heat transfer
        if (home_sph_element.is_dem_coupled())
        {
            auto dem_particle = home_sph_element.get_linked_dem_particle();
            auto dem_scalars = dem_particle.get_transfer_scalars();
            dem_scalars.add_heat_transfer(heat_transfer);
        }
 
        if (near_sph_element.is_dem_coupled())
        {
            auto dem_particle = near_sph_element.get_linked_dem_particle();
            auto dem_scalars = dem_particle.get_transfer_scalars();
            dem_scalars.add_heat_transfer(-heat_transfer);
        }
    }
};
 
template <class Backend>
struct HeatTransferTrianglesOperator
{
    template <typename Callable>
    __host__ __device__ void operator()(
        Callable function_ptr,
        IRockyDeviceModel& device_model_api,
        IRockySPHDeviceModel& sph_model,
        IRockySPHTriangleInteraction& sph_interaction_api,
        void* udf_data)
    {
        // Compute and add heat transfer
        real_sph heat_transfer =
            function_ptr(device_model_api, sph_model, sph_interaction_api, udf_data);
 
        // reduction factor for SPH elements with multiple triangle interactions
        // (1 / number of simultaneous triangle interactions)
        heat_transfer /= max(1.0f, sph_model._impl.device_model->bnd_force_coefficient[sph_interaction_api._impl->home_index]);
 
        sph_interaction_api.get_home_element().get_scalars().add_heat_transfer(heat_transfer);
 
        // Recording data for SPH Boundary Statistics
        const auto& data_recorder = sph_model._impl.device_model->boundary_data_recorder;
        if (data_recorder
            && data_recorder->is_enabled(SPHBoundaryDataRecorder::tSPHTriangleHeatTransfer))
        {
            const int index =
                data_recorder->get_scalar_index(SPHBoundaryDataRecorder::tSPHTriangleHeatTransfer);
            auto triangle_scalars = sph_interaction_api.get_near_triangle().get_scalars();
            triangle_scalars.add_scalar(index, heat_transfer * device_model_api.get_timestep());
        }
 
        // Accumulation in a scalar used for heat coupling (positive when triangles receive heat)
        auto* triangle_scalars = device_model_api._impl.model->triangle_scalars;
        if (triangle_scalars && triangle_scalars->is_heat_exchanged_enabled())
        {
            Backend::AtomicAdd(
                &triangle_scalars->get_heat_exchanged(sph_interaction_api._impl->near_index),
                static_cast<double>(-heat_transfer * device_model_api.get_timestep()));
        }
    }
};
 
template <class Backend>
struct ForceElementInteractionOperator
{
    template <typename Callable>
    __host__ __device__ void operator()(
        Callable function_ptr,
        IRockyDeviceModel& device_model_api,
        IRockySPHDeviceModel& sph_model,
        IRockySPHElementInteraction& sph_interaction_api,
        void* udf_data)
    {
        // custom calculation of the element-to-element acceleration
        const real_sph3 acceleration =
            function_ptr(device_model_api, sph_model, sph_interaction_api, udf_data);
 
        sph_interaction_api.add_acceleration(acceleration);
    }
};
 
template <class Backend>
struct ForceTriangleInteractionOperator
{
    template <typename Callable>
    __host__ __device__ void operator()(
        Callable function_ptr,
        IRockyDeviceModel& device_model_api,
        IRockySPHDeviceModel& sph_model,
        IRockySPHTriangleInteraction& sph_interaction_api,
        void* udf_data)
    {
        // custom calculation of the element-to-triangle acceleration
        const real_sph3 acceleration =
            function_ptr(device_model_api, sph_model, sph_interaction_api, udf_data);
 
        sph_interaction_api.add_acceleration(acceleration);
    }
};
 
template <class Backend>
struct ForceElementOperator
{
    template <typename Callable>
    __host__ __device__ void operator()(
        Callable function_ptr,
        IRockyDeviceModel& device_model_api,
        IRockySPHDeviceModel& sph_model,
        IRockySPHElement& sph_element_api,
        void* udf_data)
    {
        // custom calculation of the element acceleration (due to a body force)
        const real_sph3 acceleration =
            function_ptr(device_model_api, sph_model, sph_element_api, udf_data);
 
        sph_element_api.add_acceleration(acceleration);
    }
};
 
template <class backend, typename Operator, typename Callable>
__host__ __device__ void launch_sph_particles_function(
    SDeviceModel* device_data,
    SPHDeviceModel* sph_device,
    double current_time,
    Callable function_ptr,
    void* udf_data)
{
    IRockyDeviceModel device_model_api;
    device_model_api._impl.model = device_data;
    device_model_api._impl.current_time = current_time;
    IRockySPHDeviceModel sph_model(device_data, sph_device->gen_sph, sph_device);
    cuda_particle_sph* sph_elements = sph_device->prt_sph;
 
    int thread_start, thread_step, thread_stop;
    backend::ThreadIteration(
        thread_start, thread_step, thread_stop, sph_device->gen_sph->nprt_active);
    for (int i = thread_start; i < thread_stop; i += thread_step)
    {
        if (sph_elements[i].enabled() && !sph_elements[i].ghost())
        {
            // Particle Interface
            IRockySPHElement sph_particle_api { device_data, sph_device, i };
 
            // Trigger user kernel
            Operator caller{};
            caller(function_ptr, device_model_api, sph_model, sph_particle_api, udf_data);
        }
    }
}
 
template <class Backend, typename Operator, typename Callable>
__host__ __device__ void launch_sph_interactions_particles(
    SDeviceModel* device_data,
    SPHDeviceModel* sph_device,
    double current_time,
    Callable function_ptr,
    void* udf_data,
    IterationType it_type)
{
    IRockyDeviceModel device_model_api;
    device_model_api._impl.model = device_data;
    device_model_api._impl.current_time = current_time;
    IRockySPHDeviceModel sph_model(device_data, sph_device->gen_sph, sph_device);
 
    RockySPHElementInteraction sph_interaction(device_data, sph_device, current_time);
    IRockySPHElementInteraction sph_interaction_api { &sph_interaction };
 
    cuda_general_sph* gen_sph = sph_device->gen_sph;
    const auto& periodic_domain = device_data->gen->periodic_domain;
 
    // Local storage for interactions and periodics
    int* nn_list = nullptr;
    PeriodicPosition* nn_periodic = nullptr;
    if (!gen_sph->use_particles_neighbors_list)
    {
        nn_list = sph_device->local_n_list + Backend::ThreadIndex() * gen_sph->num_neighbors_max;
        if (periodic_domain.enabled())
            nn_periodic =
                sph_device->local_n_periodic + Backend::ThreadIndex() * gen_sph->num_neighbors_max;
    }
 
    // Thread
    // TODO (ROCKY-20361): This will no longer be needed once we use SPH Search Handlers
    const int n_entities = [&gen_sph, it_type] {
        if (!gen_sph->use_particles_neighbors_list)
            return gen_sph->nprt_active;
        if (it_type == itFreeAndFreeCoupled)
        {
            return gen_sph->ncont_ff_p_int + gen_sph->ncont_fc_p_int;
        }
        else if (it_type == itFreeCoupledOnly)
        {
            return gen_sph->ncont_fc_p_int;
        }
        else  // itAll
        {
            return gen_sph->ncont_ff_p_int + gen_sph->ncont_fc_p_int + gen_sph->ncont_cc_p_int;
        }
    }();
    int thread_start, thread_step, thread_stop;
    Backend::ThreadIteration(thread_start, thread_step, thread_stop, n_entities);
    for (int index = thread_start; index < thread_stop; index += thread_step)
    {
        // Handling different interactionhoods structures
        cuda_sph_contact* _sph_interaction = nullptr;
        int n_iter, home_index, near_index, offset = 0;
        if (gen_sph->use_particles_neighbors_list)
        {
            _sph_interaction = (it_type == itFreeCoupledOnly)
                ? sph_device->cnt_sph_p + sph_device->ind_sph_p[index + gen_sph->ncont_ff_p_int]
                : sph_device->cnt_sph_p + sph_device->ind_sph_p[index];
 
            home_index = _sph_interaction->nh;
            near_index = _sph_interaction->nn;
            n_iter = 1;
        }
        else
        {
            int nn_ff, nn_fc, nn_cc;
            if (periodic_domain.enabled())
            {
                sph::get_periodic_neighbors(home_index, device_data->ind_per[index], device_data->gen, sph_device, current_time, 0.0f,
                    nn_ff, nn_fc, nn_cc, nn_list, nn_periodic);
                if(gen_sph->nprt_active_plus_halo > gen_sph->nprt_active)
                    sph::get_periodic_multi_GPU_neighbors(home_index, device_data->ind_per[index], device_data->gen, sph_device, current_time, 0.0f, 
                        nn_ff, nn_fc, nn_cc, nn_list, nn_periodic);
                sph::partition_array(home_index, sph_device, nn_ff, nn_fc, nn_cc, nn_list, nn_periodic);
            }
            else
            {
                sph::get_non_periodic_neighbors(home_index, index, device_data->gen, sph_device, current_time, 0.0f, nn_ff, nn_fc, nn_cc, nn_list);
                if(gen_sph->nprt_active_plus_halo > gen_sph->nprt_active)
                    sph::get_non_periodic_multi_GPU_neighbors(home_index, device_data->gen, sph_device, current_time, 0.0f, nn_ff, nn_fc, nn_cc, nn_list);
                sph::partition_array(home_index, sph_device, nn_ff, nn_fc, nn_cc, nn_list, nullptr);
            }
            if (it_type == itFreeAndFreeCoupled)
            {
                n_iter = nn_ff + nn_fc;
            }
            else if (it_type == itFreeCoupledOnly)
            {
                n_iter = nn_fc;
                offset = nn_ff;
            }
            else  // itAll
            {
                n_iter = nn_ff + nn_fc + nn_cc;
            }
        }
 
        // Home Interaction Interface
        sph_interaction.reset_home_element(home_index);
 
        // Loop over number of active interactions
        for (int i = 0; i < n_iter; ++i)
        {
            // Periodic Stuff
            PeriodicPosition cnt_per {};
            if (!gen_sph->use_particles_neighbors_list)
            {
                near_index = nn_list[i + offset];
                if (periodic_domain.enabled())
                    cnt_per = nn_periodic[i + offset];
            }
            else if (periodic_domain.enabled())
                cnt_per = (it_type == itFreeCoupledOnly)
                    ? sph_device
                          ->per_sph_cnt_p[sph_device->ind_sph_p[index + gen_sph->ncont_ff_p_int]]
                    : sph_device->per_sph_cnt_p[sph_device->ind_sph_p[index]];
 
            // Near Interaction Interface
            sph_interaction.reset_near_element(near_index, cnt_per);
 
            // Trigger user kernel
            Operator caller {};
            caller(function_ptr, device_model_api, sph_model, sph_interaction_api, udf_data);
        }
    }
}
 
template <class Backend, typename Operator, typename Callable>
__host__ __device__ void launch_sph_interactions_triangles(
    SDeviceModel* device_data,
    SPHDeviceModel* sph_device,
    double current_time,
    Callable function_ptr,
    void* udf_data)
{
    IRockyDeviceModel device_model_api;
    device_model_api._impl.model = device_data;
    device_model_api._impl.current_time = current_time;
    IRockySPHDeviceModel sph_model(device_data, sph_device->gen_sph, sph_device);
 
    RockySPHTriangleInteraction sph_interaction(device_data, sph_device, current_time);
    IRockySPHTriangleInteraction sph_interaction_api { &sph_interaction };
 
    const auto& periodic_domain = device_data->gen->periodic_domain;
 
    // loop only over element/triangle pairs with projection inside the triangle
    int thread_start, thread_step, thread_stop;
    Backend::ThreadIteration(thread_start, thread_step, thread_stop, sph_device->gen_sph->n_triangle_force_pairs);
    for (int i = thread_start; i < thread_stop; i += thread_step)
    {
        const int index = sph_device->triangle_force_pairs[i];
        cuda_sph_contact* _sph_interaction = sph_device->cnt_sph_t + index;
        PeriodicPosition cnt_per {};
        if (periodic_domain.enabled())
            cnt_per = sph_device->per_sph_cnt_t[index];
 
        // Interaction Interface
        sph_interaction.reset_triangle_contact(_sph_interaction, cnt_per);
 
        // ROCKY-17943: For periodic cases, skipping interaction w/ triangles, such that
        // the projected point is positioned outside periodic domain
        if (periodic_domain.enabled()
            && periodic_domain.is_outside_cartesian_periodic_bounds(sph_interaction.inside_point))
        {
            continue;
        }
 
        // Trigger user kernel
        Operator caller {};
        caller(function_ptr, device_model_api, sph_model, sph_interaction_api, udf_data);
    }
}
 
}  // namespace device
Rocky Solver SDK 2024 R2

api_device_sph

Connect with Ansys
