api_device_cfd_coupling.hpp | Ansys Developer Portal

#pragma once
 
// Includes =================================================================================================
#include <rocky20/api/rocky_api_cfd_functions.h>
#include <rocky20/device/device_cfd_drag.hpp>
#include <rocky20/device/device_cfd_coupling.hpp>
 
namespace device
{
    
template<typename... Op>
__host__ __device__
void execute_operators(
    const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data,
    Op... ops)
{
    int dummy[sizeof...(Op)] = { (ops(particle, cfd, output, data), 0)... };
}
 
/*
A utility structure to hold all device data used to execute the coupling kernel
All data is a device data, so only the final kernel would be able to access such data
*/
struct CouplingDispatchData
{
    bool is_cuda{ false }, collect_cfd_data{ false };
 
    SDeviceModel* device_data;
    unsigned int* particle_indexes;
    int n_particles;
    void* udf_data;
};
 
/*
Coupling kernel 
*/
template <class backend, typename... Op>
__host__ __device__
void api_coupling_kernel(CouplingDispatchData& data, Op... operators)
{
    auto device_data = data.device_data;
 
    IRockyParticle particle_api(-1, device_data);
 
    IRockyCFDProperties cfd{};
    cfd.model = device_data;
    double3 gravity { device_data->gen->gx, device_data->gen->gy, device_data->gen->gz };
 
#ifdef __CUDA_ARCH__
    __shared__ CFDTimestepBin block_bins[CD_NTHREADS];
    const int thread_index = threadIdx.x;
    CFDTimestepBin& block_bin = block_bins[thread_index];
#else
    CFDTimestepBin& block_bin = device_data->cfd_monitor->thread_bins[backend::ThreadIndex()];
#endif
 
    block_bin.reset();
 
    int thread_start, thread_step, thread_stop;
    backend::ThreadIteration(thread_start, thread_step, thread_stop, data.n_particles);
    for (auto i = thread_start; i < thread_stop; i += thread_step)
    {
        // particle
        const int ip = data.particle_indexes[i];
        particle_api._impl.reset(ip);
        if (!particle_api._impl.particle->enabled()) {
            continue;
        }
 
        const auto ic = device_data->particle_scalars->get_cfd_index(ip);
        if (ic < 0) {
            continue;
        }
 
        cfd.reset(ip, ic);
        CFDCouplingOutput output{};
 
        // Relative Buoyancy force, due to the relative fluid density 
        output.set_relative_buoyancy_force(
            -cfd.get_particle_volume() * cfd.model->gen_cfd->relative_density * gravity
            );
 
        // Force due to pressure gradient (accounts for buoyancy and flow-induced
        // pressure variations).
        output.set_pressure_gradient_force(
            -cfd.get_particle_volume() * cfd.get_fluid_pressure_gradient()
            );
 
        execute_operators(particle_api, cfd, output, data.udf_data, operators...);
 
        // add forces and moments to the particle
        auto& particle_forces = device_data->frc[ip];
        particle_forces.fx += output.total_cfd_force.x;
        particle_forces.fy += output.total_cfd_force.y;
        particle_forces.fz += output.total_cfd_force.z;
 
        if (output.has_moment()) {
            particle_forces.mx += output.total_cfd_moment.x;
            particle_forces.my += output.total_cfd_moment.y;
            particle_forces.mz += output.total_cfd_moment.z;
        }
 
        // add critical time-step ratio to CFD Monitor
        block_bin.check_critical_timestep(
            device_data->get_general().dt,
            output.drag_coefficient_fluent / particle_api.get_mass(),
            device_data->cfd_monitor->levels
        );
    }
 
#ifdef __CUDA_ARCH__
    // Block reduce
    device::cuda_block_reduce(
        block_bins, device_data->cfd_monitor->thread_bins,
        [](CFDTimestepBin &left, const CFDTimestepBin& right)
        {
            left.join(right);
        }
    );
#endif
}
 
struct EmptyCustomCouplingPlugin {};
 
struct InvalidOperator
{
    InvalidOperator()
    {
        #ifndef __CUDA_ARCH__
        ROCKY_RUNTIME_ERROR("Invalid operation selection");
        #endif
    }
 
    __host__ __device__ inline
    void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {}
};
 
 
template <class Plugin>
auto custom_drag_operator() { return InvalidOperator(); }
 
template <class Plugin>
auto custom_lift_operator() { return InvalidOperator(); }
 
template <class Plugin>
auto custom_virtual_mass_operator(SDeviceModel *device_data) { return InvalidOperator(); }
 
template <class Plugin>
auto custom_torque_operator() { return InvalidOperator(); }
 
template <class Plugin>
auto custom_convection_operator(SDeviceModel* device_data) { return InvalidOperator(); }
 
/*
Coupling Operators: executed, when enabled in the following order:
    - Drag
    - Lift
    - Virtual Mass
    - Torque
    - Heat Convection
    - CFD Data recorder
    - Two Way Coupling
    - Turbulence
*/
struct DragOperator
{
    __host__ __device__ inline
    void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
        device::compute_cfd_drag(particle, cfd, output, data);
    }
};
 
 
 
template <typename _DragLawOperator>
struct DragParametricOperator
{
    __host__ __device__ inline
        void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
        DragLaw<_DragLawOperator>(particle, cfd, output, data);
    }
};
 
// Operator defined with drag laws implementation with light weight code, with less impact of performance.
struct DragLightWeightOperator
{
    static bool valid(const ECFDDragLaw drag_law)
    {
        switch (drag_law)
        {
        case cfd_drag_law_SchillerNaumann1933:
        case cfd_drag_law_Dallavalle1948:
        case cfd_drag_law_HaiderLevenspiel1989:
        case cfd_drag_law_MorsiAlexander1972:
        case cfd_drag_law_WenYu1966:
        case cfd_drag_law_Ergun1958:
        case cfd_drag_law_GidaspowBezburuahDing1992:
        case cfd_drag_law_HuilinGidaspow2003:
        case cfd_drag_law_DiFelice1994:
        case cfd_drag_law_HillKochLadd2001:
        case cfd_drag_law_HillKoch2001:
        case cfd_drag_law_SyamlalOBrien1987:
            return true;
        default:
            return false;
        }
    }
 
    __host__ __device__ inline
        void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
        switch (cfd.cfd_particle_group->drag_law)
        {
        case cfd_drag_law_SchillerNaumann1933: DragLaw<SchillerNaumann1933Operator>(particle, cfd, output, data); break;
        case cfd_drag_law_Dallavalle1948: DragLaw<Dallavalle1948Operator>(particle, cfd, output, data); break;
        case cfd_drag_law_HaiderLevenspiel1989: DragLaw<HaiderLevenspiel1989Operator>(particle, cfd, output, data); break;
        case cfd_drag_law_MorsiAlexander1972: DragLaw<MorsiAlexander1972Operator>(particle, cfd, output, data); break;
        case cfd_drag_law_WenYu1966: DragLaw<WenYu1966Operator>(particle, cfd, output, data); break;
        case cfd_drag_law_Ergun1958: DragLaw<Ergun1958Operator>(particle, cfd, output, data); break;
        case cfd_drag_law_GidaspowBezburuahDing1992: DragLaw<GidaspowBezburuahDing1992Operator>(particle, cfd, output, data); break;
        case cfd_drag_law_HuilinGidaspow2003: DragLaw<HuilinGidaspow2003Operator>(particle, cfd, output, data); break;
        case cfd_drag_law_DiFelice1994: DragLaw<DiFelice1994Operator>(particle, cfd, output, data); break;
        case cfd_drag_law_HillKochLadd2001: DragLaw<HillKochLadd2001Operator>(particle, cfd, output, data); break;
        case cfd_drag_law_HillKoch2001: DragLaw<HillKoch2001Operator>(particle, cfd, output, data); break;
        case cfd_drag_law_SyamlalOBrien1987: DragLaw<SyamlalOBrien1987Operator>(particle, cfd, output, data); break;
        }
    }
};
 
// Operator defined with drag laws implementation with light weight code, with less impact of performance.
struct DragHeavyWeightOperator
{
    static bool valid(const ECFDDragLaw drag_law)
    {
        switch (drag_law)
        {
        case cfd_drag_law_Ganser1993:
        case cfd_drag_law_MarheinekeWegener2011:
            return true;
        default:
            return false;
        }
    }
 
    __host__ __device__ inline
        void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
        switch (cfd.cfd_particle_group->drag_law)
        {
        case cfd_drag_law_Ganser1993: Ganser1993Operator::DragLaw(particle, cfd, output); break;
        case cfd_drag_law_MarheinekeWegener2011: MarheinekeWegener2011Operator::DragLaw(particle, cfd, output); break;
        }
    }
};
 
struct LiftOperator
{
    __host__ __device__ inline
    void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
        device::compute_cfd_lift(particle, cfd, output, data);
    }
};
 
template <typename _LiftCoeff>
struct LiftParametricOperator
{
    __host__ __device__ inline
        void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
        device::LiftLaw<_LiftCoeff>(particle, cfd, output, data);
    }
};
 
struct VirtualMassOperator
{
    SDeviceModel* device_data;
 
    __host__ __device__ inline
        void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
        device::compute_cfd_virtual_mass(particle, cfd, output, data);
        device_data->particle_scalars->get_virtual_mass(particle._impl.particle_index) = output.virtual_mass;
    }
};
 
template <typename _Coeff>
struct VirtualMassParametricOperator
{
    SDeviceModel* device_data;
 
    __host__ __device__ inline
        void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
        device::VirtualMassLaw<_Coeff>(particle, cfd, output, data);
        device_data->particle_scalars->get_virtual_mass(particle._impl.particle_index) = output.virtual_mass;
    }
};
 
struct TorqueOperator
{
    __host__ __device__ inline
    void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
        device::compute_cfd_torque(particle, cfd, output, data);
    }
};
 
template <typename _Coeff>
struct TorqueParametricOperator
{
    __host__ __device__ inline
        void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
        device::TorqueLaw<_Coeff>(particle, cfd, output, data);
    }
};
 
struct HeatConvectionOperator
{
    SDeviceModel* device_data;
 
    __host__ __device__ inline
    void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
        device::compute_cfd_convection(particle, cfd, output, data);
        device_data->particle_transfer_scalars->get_heat_transfer(particle._impl.particle_index) += output.heat_transfer;
    }
};
 
template <typename _Coeff>
struct HeatConvectionParametricOperator
{
    SDeviceModel* device_data;
 
    __host__ __device__ inline
        void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
        device::ConvectionLaw<_Coeff>(particle, cfd, output, data);
        device_data->particle_transfer_scalars->get_heat_transfer(particle._impl.particle_index) += output.heat_transfer;
    }
};
 
struct CFDDataRecorderOperator
{
    SDeviceModel* device_data;
 
    __host__ __device__ inline
    void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
#ifndef CFD_NO_CUSTOM_LAWS
        if (device_data->cfd_data_recorder != nullptr && device_data->cfd_data_recorder->enabled)
#endif
            device_data->cfd_data_recorder->register_coupling_data(
            output, device_data->particle_scalars, particle._impl.particle_index, device_data->gen->dt
        );
    }
};
 
struct TwoWayCouplingOperator
{
    SDeviceModel* device_data;
 
    __host__ __device__ inline void operator()(
        const IRockyParticle& particle,
        const IRockyCFDProperties& cfd,
        CFDCouplingOutput& output,
        void* data)
    {
        if (device_data->gen_cfd->cfd_coupling_mode == cfdcTwoWay)
        {
            auto dt = device_data->get_general().dt;
            auto particle_scalars = particle.get_scalars();
 
            const auto& particle_velocity = particle.get_translational_velocity();
            const auto& relative_velocity = cfd.get_relative_velocity();
 
            if (device_data->gen_cfd->use_2023R2_cfd_sources)
            {
                particle_scalars.add_explicit_fluid_momentum(
                    particle_velocity * output.drag_coefficient_fluent * dt);
            }
            else
            {
                particle_scalars.add_explicit_fluid_momentum(
                    -output.drag_coefficient_fluent * relative_velocity * dt);
            }
            particle_scalars.add_implicit_fluid_momentum(-output.drag_coefficient_fluent * dt);
 
            if (output.has_explicit_forces())
            {
                particle_scalars.add_explicit_fluid_momentum(
                    -output.momentum_source_explicit_fluent * dt);
            }
            if (output.has_heat())
            {
                if (device_data->gen_cfd->use_2023R2_cfd_sources)
                {
                    particle_scalars.add_explicit_fluid_heat_transfer_rate(
                        output.heat_source_explicit_fluent * dt);
                }
                else
                {
                    particle_scalars.add_explicit_fluid_heat_transfer_rate(
                        (output.heat_source_explicit_fluent
                         + output.heat_source_implicit_fluent * cfd.get_fluid_temperature())
                        * dt);
                }
                particle_scalars.add_implicit_fluid_heat_transfer_rate(
                    output.heat_source_implicit_fluent * dt);
            }
        }
    }
};
 
struct TurbulenceOperator
{
    SDeviceModel* device_data;
 
    __host__ __device__ inline
    void operator()(const IRockyParticle& particle, const IRockyCFDProperties& cfd, CFDCouplingOutput& output, void* data)
    {
#ifndef CFD_NO_CUSTOM_LAWS
        if (device_data->gen_cfd->use_turbulent_dispersion)
#endif
        update_turbulent_dispersion(
            device_data->particle_scalars,
            particle._impl.particle_index,
            cfd.cfd_index,
            device_data->get_general().dt,
            device_data->prt_vel,
            device_data->fluid_scalars,
            cfd.props,
            particle._impl.get_group(),
            output.C_D,
            device_data->cfd_random
        );
    }
};
 
/*
Dispatch function, a sequence of selector function calls, that will at the end collected all operators in a 
variadic template list, used to instantiate and dispatch a kernel.
*/
 
template <class backend, typename... Operators>
__global__
void cuda_coupling_kernel(CouplingDispatchData data, Operators... operators)
{
    api_coupling_kernel<backend>(data, operators...);
}
 
template <class backend, typename... Operators>
void omp_coupling_kernel(CouplingDispatchData &data, Operators... operators)
{
    // shut up nvcc compile warnings when compiling device function version using OMPBackend
    #ifndef __CUDA_ARCH__
 
    api_coupling_kernel<backend>(data, operators...);
 
    #endif
}
 
template <typename... Operators> inline
void dispatch_coupling_kernel(CouplingDispatchData &data, Operators... operators)
{
    if (data.is_cuda) {
#ifdef __CUDACC__
        cuda_coupling_kernel<CudaBackend><<<CD_NBLOCKS, CD_NTHREADS>>>(data, operators...);
#else
        ROCKY_RUNTIME_ERROR("Invalid GPU kernel execution.");
#endif
    }
    else
        omp_coupling_kernel<OMPBackend>(data, operators...);
}
 
/*
Selector functions, called in ascending order
e.g. drag -> lift -> virtual mass -> ....
*/
template<typename... Operators>
void select_turbulence_and(CouplingDispatchData& data, const cuda_general_cfd& h_gen_cfd, Operators... operators)
{
    if (h_gen_cfd.use_turbulent_dispersion)
        dispatch_coupling_kernel(data, operators..., TurbulenceOperator{ data.device_data });
    else
        dispatch_coupling_kernel(data, operators...);
}
 
template<typename... Operators>
void select_coupling_mode_and(CouplingDispatchData& data, const cuda_general_cfd& h_gen_cfd, Operators... operators)
{
    if (h_gen_cfd.cfd_coupling_mode == cfdcTwoWay)
        select_turbulence_and(
            data, h_gen_cfd, operators..., TwoWayCouplingOperator{ data.device_data });
    else
        select_turbulence_and(data, h_gen_cfd, operators...);
}
 
template<typename... Operators>
void select_cfd_data_recorder_and(CouplingDispatchData& data, const cuda_general_cfd& h_gen_cfd, Operators... operators)
{
#ifdef CFD_NO_CUSTOM_LAWS
    if (data.collect_cfd_data)
        select_coupling_mode_and(data, h_gen_cfd, operators..., CFDDataRecorderOperator{ data.device_data });
    else
        select_coupling_mode_and(data, h_gen_cfd, operators...);
#else
    // When compiling the kernel with custom laws, we must limit the number of times we instantiate the coupling kernel.
    // To avoid code bloat, we create here just a single one with all remaining functions enabled
    dispatch_coupling_kernel(
        data,
        operators...,
        CFDDataRecorderOperator{ data.device_data },
        TwoWayCouplingOperator{ data.device_data },
        TurbulenceOperator{ data.device_data }
    );
#endif
}
 
template<class Plugin, typename... Operators>
void select_convection_and(CouplingDispatchData& data, const cuda_general_cfd& h_gen_cfd, const CFDParticleGroup& group, Operators... operators)
{
    switch (group.convection_law)
    {
        case cfd_convection_law_none:
            select_cfd_data_recorder_and(data, h_gen_cfd, operators...);
            break;
#ifndef CFD_NO_CUSTOM_LAWS
        case cfd_convection_law_custom:
            select_cfd_data_recorder_and(data, h_gen_cfd, operators..., custom_convection_operator<Plugin>(data.device_data));
            break;
#endif
        default:
            select_cfd_data_recorder_and(data, h_gen_cfd, operators..., HeatConvectionOperator{ data.device_data });
    }
}
 
template<class Plugin, typename... Operators>
void select_torque_and(CouplingDispatchData& data, const cuda_general_cfd& h_gen_cfd, const CFDParticleGroup& group, Operators... operators)
{
    switch (group.torque_law)
    {
        case cfd_torque_law_none:
            select_convection_and<Plugin>(data, h_gen_cfd, group, operators...);
            break;
        case cfd_torque_law_DennisSinghIngham1980:
            select_convection_and<Plugin>(data, h_gen_cfd, group, operators..., TorqueParametricOperator<DennisSinghIngham1980Operator>());
            break;
#ifndef CFD_NO_CUSTOM_LAWS
        case cfd_torque_law_custom:
            select_convection_and<Plugin>(data, h_gen_cfd, group, operators..., custom_torque_operator<Plugin>());
            break;
#endif
    }
}
 
template<class Plugin, typename... Operators>
void select_virtual_mass_and(CouplingDispatchData& data, const cuda_general_cfd& h_gen_cfd, const CFDParticleGroup& group, Operators... operators)
{
    switch (group.virtual_mass_law)
    {
        case cfd_virtual_mass_law_none:
            select_torque_and<Plugin>(data, h_gen_cfd, group, operators...);
            break;
#ifndef CFD_NO_CUSTOM_LAWS
        case cfd_virtual_mass_law_custom:
            select_torque_and<Plugin>(data, h_gen_cfd, group, operators..., custom_virtual_mass_operator<Plugin>(data.device_data));
            break;
#endif
        default:
            select_torque_and<Plugin>(data, h_gen_cfd, group, operators..., VirtualMassOperator{ data.device_data });
    }
}
 
template<class Plugin, typename... Operators>
void select_lift_and(CouplingDispatchData& data, const cuda_general_cfd& h_gen_cfd, const CFDParticleGroup& h_group, Operators... operators)
{
    switch (h_group.lift_law)
    {
        case cfd_lift_law_none:
            select_virtual_mass_and<Plugin>(data, h_gen_cfd, h_group, operators...);
            break;
        //case cfd_lift_law_Mei1992:
        //    select_virtual_mass_and<Plugin>(data, h_gen_cfd, h_group, operators..., LiftOperator<Mei1992Operator>());
        //    break;
        //case cfd_lift_law_Saffman1968:
        //    select_virtual_mass_and<Plugin>(data, h_gen_cfd, h_group, operators..., LiftOperator<Saffman1968Operator>());
        //    break;
#ifndef CFD_NO_CUSTOM_LAWS
        case cfd_lift_law_custom:
            select_virtual_mass_and<Plugin>(data, h_gen_cfd, h_group, operators..., custom_lift_operator<Plugin>());
            break;
#endif
        default:
            select_virtual_mass_and<Plugin>(data, h_gen_cfd, h_group, operators..., LiftOperator());
    }
}
 
template <class Plugin>
void select_drag_and(CouplingDispatchData& data, const cuda_general_cfd& h_gen_cfd, const CFDParticleGroup& h_group)
{
#ifndef CFD_NO_CUSTOM_LAWS
    switch (h_group.drag_law)
    {
    case cfd_drag_law_custom:
        select_lift_and<Plugin>(data, h_gen_cfd, h_group, custom_drag_operator<Plugin>());
        break;
    default: {
#else
    {{
#endif
            if (DragLightWeightOperator::valid(h_group.drag_law))
                select_lift_and<Plugin>(data, h_gen_cfd, h_group, DragLightWeightOperator());
            else if (DragHeavyWeightOperator::valid(h_group.drag_law))
                select_lift_and<Plugin>(data, h_gen_cfd, h_group, DragHeavyWeightOperator());
            else
                select_lift_and<Plugin>(data, h_gen_cfd, h_group, DragOperator());
    }
    }
}
 
template <class Plugin>
void select_and_dispatch_coupling_kernel(
    bool is_cuda, bool collect_cfd_data,
    SDeviceModel* device_data, unsigned int* particle_indexes, int n_particles, void* udf_data,
    const cuda_general_cfd& h_gen_cfd, const CFDParticleGroup& h_group
)
{
    CouplingDispatchData data{ is_cuda, collect_cfd_data, device_data, particle_indexes, n_particles, udf_data };
    select_drag_and<Plugin>(data, h_gen_cfd, h_group);
}
 
template <class backend, typename _Op>
__host__ __device__ void api_compute_cfd_data_coupling(
    SDeviceModel* device_data, double current_time, _Op op, void* udf_data)
{
    IRockyDeviceModel device_model_api;
    device_model_api._impl.model = device_data;
    device_model_api._impl.current_time = current_time;
 
    IRockyParticle particle_api(-1, device_data);
 
    IRockyCFDProperties cfd{};
    cfd.model = device_data;
 
    const cuda_general& general = device_data->get_general();
 
    int thread_start, thread_step, thread_stop;
    backend::ThreadIteration(thread_start, thread_step, thread_stop, general.nprt_active);
    for (int i = thread_start; i < thread_stop; i += thread_step)
    {
        particle_api._impl.reset(i);
        int ind_cfd = device_data->particle_scalars->get_cfd_index(i);
        if (!particle_api._impl.particle->enabled() || ind_cfd < 0)
        {
            continue;
        }
 
        cfd.reset(i, ind_cfd);
 
        op(device_model_api, particle_api, cfd, udf_data);
    }
}
 
} // namespace device
Rocky Solver SDK 2024 R2

api_device_cfd_coupling

Connect with Ansys
