deltaSPH.cl

Go to the documentation of this file.

/*
 *  This file is part of AQUAgpusph, a free CFD program based on SPH.
 *  Copyright (C) 2012  Jose Luis Cercos Pita <jl.cercos@upm.es>
 *
 *  AQUAgpusph is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  AQUAgpusph is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with AQUAgpusph.  If not, see <http://www.gnu.org/licenses/>.
 */


 
#ifndef EXCLUDED_PARTICLE
    #define EXCLUDED_PARTICLE(index) imove[index] <= 0
#endif
 
#if defined(LOCAL_MEM_SIZE) && defined(NO_LOCAL_MEM)
    #error NO_LOCAL_MEM has been set.
#endif
 
#include "resources/Scripts/types/types.h"
#include "resources/Scripts/KernelFunctions/Kernel.h"

__kernel void simple(const __global unsigned int* iset,
                     const __global int* imove,
                     __global vec* lap_p_corr,
                     __constant float* refd,
                     usize N,
                     vec g)
{
    const usize i = get_global_id(0);
    if(i >= N)
        return;
    if(EXCLUDED_PARTICLE(i))
        return;
 
    lap_p_corr[i] = refd[iset[i]] * g;
}

__kernel void full(const __global int* imove,
                   const __global vec* r,
                   const __global float* rho,
                   const __global float* m,
                   const __global float* p,
                   __global vec* lap_p_corr,
                   usize N,
                   LINKLIST_LOCAL_PARAMS)
{
    const usize i = get_global_id(0);
    const usize it = get_local_id(0);
    if(i >= N)
        return;
    if(EXCLUDED_PARTICLE(i)){
        return;
    }
 
    const vec_xyz r_i = r[i].XYZ;
    const float p_i = p[i];
 
    // Initialize the output
    #ifndef LOCAL_MEM_SIZE
        #define _GRADP_ lap_p_corr[i].XYZ
    #else
        #define _GRADP_ lap_p_corr_l[it]
        __local vec_xyz lap_p_corr_l[LOCAL_MEM_SIZE];
        _GRADP_ = VEC_ZERO.XYZ;
    #endif
 
    const usize c_i = icell[i];
    BEGIN_NEIGHS(c_i, N, n_cells, icell, ihoc){
        if( (i == j) || (EXCLUDED_PARTICLE(j))){
            j++;
            continue;
        }
        const vec_xyz r_ij = r[j].XYZ - r_i;
        const float q = length(r_ij) / H;
        if(q >= SUPPORT)
        {
            j++;
            continue;
        }
        {
            const float f_ij = kernelF(q) * CONF * m[j] / rho[j];
            _GRADP_ += (p[j] - p_i) * f_ij * r_ij;
        }
    }END_NEIGHS()
 
    #ifdef LOCAL_MEM_SIZE
        lap_p_corr[i].XYZ = _GRADP_;
    #endif
}

__kernel void full_mls(const __global int* imove,
                   const __global matrix* mls,
                   __global vec* lap_p_corr,
                   usize N)
{
    const usize i = get_global_id(0);
    if(i >= N)
        return;
    if(EXCLUDED_PARTICLE(i)){
        return;
    }
 
    lap_p_corr[i] = MATRIX_DOT(mls[i], lap_p_corr[i]);
}

__kernel void lapp(const __global int* imove,
                   const __global vec* r,
                   const __global float* rho,
                   const __global float* m,
                   const __global float* p,
                   __global float* lap_p,
                   usize N,
                   LINKLIST_LOCAL_PARAMS)
{
    const usize i = get_global_id(0);
    const usize it = get_local_id(0);
    if(i >= N)
        return;
    if(EXCLUDED_PARTICLE(i)){
        return;
    }
 
    const vec_xyz r_i = r[i].XYZ;
    const float p_i = p[i];
 
    // Initialize the output
    #ifndef LOCAL_MEM_SIZE
        #define _LAPP_ lap_p[i]
    #else
        #define _LAPP_ lap_p_l[it]
        __local float lap_p_l[LOCAL_MEM_SIZE];
        _LAPP_ = 0.f;
    #endif
 
    const usize c_i = icell[i];
    BEGIN_NEIGHS(c_i, N, n_cells, icell, ihoc){
        if( (i == j) || (EXCLUDED_PARTICLE(j))){
            j++;
            continue;
        }
        const vec_xyz r_ij = r[j].XYZ - r_i;
        const float q = length(r_ij) / H;
        if(q >= SUPPORT)
        {
            j++;
            continue;
        }
        {
            const float f_ij = kernelF(q) * CONF * m[j] / rho[j];
            _LAPP_ += (p[j] - p_i) * f_ij;
        }
    }END_NEIGHS()
 
    #ifdef LOCAL_MEM_SIZE
        lap_p[i] = _LAPP_;
    #endif
}

__kernel void lapp_corr(const __global int* imove,
                        const __global vec* r,
                        const __global float* rho,
                        const __global float* m,
                        const __global vec* lap_p_corr,
                        __global float* lap_p,
                        usize N,
                        LINKLIST_LOCAL_PARAMS)
{
    const usize i = get_global_id(0);
    const usize it = get_local_id(0);
    if(i >= N)
        return;
    if(EXCLUDED_PARTICLE(i)){
        return;
    }
 
    const vec_xyz r_i = r[i].XYZ;
    const vec_xyz gradp_i = lap_p_corr[i].XYZ;
 
    // Initialize the output
    #ifndef LOCAL_MEM_SIZE
        #define _LAPP_ lap_p[i]
    #else
        #define _LAPP_ lap_p_l[it]
        __local float lap_p_l[LOCAL_MEM_SIZE];
        _LAPP_ = lap_p[i];
    #endif
 
    const usize c_i = icell[i];
    BEGIN_NEIGHS(c_i, N, n_cells, icell, ihoc){
        if( (i == j) || (EXCLUDED_PARTICLE(j))){
            j++;
            continue;
        }
        const vec_xyz r_ij = r[j].XYZ - r_i;
        const float q = length(r_ij) / H;
        if(q >= SUPPORT)
        {
            j++;
            continue;
        }
        {
            const vec_xyz gradp_ij = lap_p_corr[j].XYZ + gradp_i;
            const float f_ij = kernelF(q) * CONF * m[j] / rho[j];
            _LAPP_ -= 0.5f * dot(gradp_ij, r_ij) * f_ij;
        }
    }END_NEIGHS()
 
    #ifdef LOCAL_MEM_SIZE
        lap_p[i] = _LAPP_;
    #endif
}

__kernel void deltaSPH(const __global unsigned int* iset,
                       const __global int* imove,
                       const __global float* rho,
                       const __global float* lap_p,
                       __global float* drhodt,
                       __constant float* refd,
                       __constant float* delta,
                       usize N,
                       float dt)
{
    const usize i = get_global_id(0);
    if(i >= N)
        return;
    if(EXCLUDED_PARTICLE(i))
        return;
 
    const uint set_i = iset[i];
    const float delta_f = delta[set_i] * dt * rho[i] / refd[set_i];
 
    drhodt[i] += delta_f * lap_p[i];
}
 
/*
 * @}
 */