Reaction diffusion on mesh surface

        kernel void reactionDiffusion( 
                 int nbs_length, 
                 global int * nbs_index, 
                 global int * nbs ,
                 int A_length, 
                 global float * A ,
                 int B_length, 
                 global float * B ,
                 int __A_length, 
                 global float * __A ,
                 int __B_length, 
                 global float * __B ,
                 float  da ,
                 float  db ,
                 float  feed ,
                 float  kill 
)
{   int idx = get_global_id(0);
    if (idx >= A_length)
        return;
    
    float a = vload(idx, A);
    float b = vload(idx, B);
    
    float sum_a = 0.0f;
    float sum_b = 0.0f;
    
    int start = nbs_index[idx];
    int end = nbs_index[idx + 1];
    for(int i = start; i < end; i++){
        int nb = nbs[i];
        sum_a += vload(nb, A);
        sum_b += vload(nb, B);
    }
    
    sum_a = sum_a / (end - start) + a * -1.0f;
    sum_b = sum_b / (end - start) + b * -1.0f;
    
    float reaction = a * b * b;
    a += da * sum_a - reaction + feed * (1.0f - a);
    b += db * sum_b + reaction - (kill + feed) * b;
    a = clamp(a, 0.0f, 1.0f);
    b = clamp(b, 0.0f, 1.0f);
    
    vstore(a, idx, __A);
    vstore(b, idx, __B);
}

/////////// WRITE BACK /////////////

kernel void writeBack( 
                 int nbs_length, 
                 global int * nbs_index, 
                 global int * nbs ,
                 int A_length, 
                 global float * A ,
                 int B_length, 
                 global float * B ,
                 int __A_length, 
                 global float * __A ,
                 int __B_length, 
                 global float * __B ,
                 float  da ,
                 float  db ,
                 float  feed ,
                 float  kill 
)
{   int idx = get_global_id(0);
    if (idx >= A_length)
        return;
    
    float a = vload(idx, __A);
    float b = vload(idx, __B);
    
    vstore(a, idx, A);
    vstore(b, idx, B);
}
    

        float amount = chf('amount');

v@P += v@N * f@B * amount;