fix particle interaction with planktonIndividual

src/Individuals/Abiotic/particle_interaction.jl

@@ -13,63 +13,71 @@
     return dist
 end
 
-##### compare particle distance and Rd, decide whether
-##### this vesicle will be merged into phytoplankton cell
-@kernel function calc_merge_matrix_kernel!(intac, abiotic, plank, rnd, grid, abio_p)
-    i, j = @index(Global, NTuple)
-    intac[i,j] = isless(calc_distance(plank.xi[i], plank.yi[i],
-                                      plank.zi[i], plank.x[i],
-                                      plank.y[i], plank.z[i],
-                                      abiotic.xi[j], abiotic.yi[j],
-                                      abiotic.zi[j], abiotic.x[j],
-                                      abiotic.y[j], abiotic.z[j],
-                                      grid), abio_p.Rd * rnd.x[i]) * abiotic.ac[j] 
-end
-function calc_merge_matrix!(intac, abiotic, plank, rnd, grid, abio_p, arch::Architecture)
-    kernel! = calc_merge_matrix_kernel!(device(arch), (16,16), (size(intac)))
-    kernel!(intac, abiotic, plank, rnd, grid, abio_p)
-    return nothing
+##### 1. Reset Interaction 
+##### Clears the Top-K candidate list before each step.
+function reset_interaction!(intac, arch::Architecture)
+    @inbounds intac .= 0
 end
 
-##### clean the merge matrix - one vesicle might be close enough
-##### to two or more phytoplankton cells. 
-##### Only the first phytoplankton cell can uptake this vesicle.
-##### First, find all indices of 1.0f0 in intac, return an AbstracArray of
-##### Cartesian Index, e.g., CartesianIndex(phyto_ind, abiotic_ind)
-##### Then, assign phyto_ind to abiotic.merge[abiotic_ind]
-@kernel function assign_merge_kernel!(abiotic, inds)
+##### 2. Calculate Interaction Topology 
+##### Revised: Uses random slot assignment to avoid race conditions.
+@kernel function calc_interaction_kernel!(intac, plank, abiotic, rnd, grid, abio_p, max_candidates)
     i = @index(Global)
-    abiotic.merg[inds[i][2]] = inds[i][1]
+
+    if plank.ac[i] == 1.0f0 
+        for k in 1:length(abiotic.ac)
+            if abiotic.ac[k] == 1.0f0
+                dist = calc_distance(plank.xi[i], plank.yi[i], plank.zi[i], plank.x[i], plank.y[i], plank.z[i],
+                                     abiotic.xi[k], abiotic.yi[k], abiotic.zi[k], abiotic.x[k], abiotic.y[k], abiotic.z[k],
+                                     grid)
+                if dist < abio_p.Rd 
+                    slot = (i % max_candidates) + 1
+                    intac[slot, k] = i
+                end
+            end
+        end
+    end 
 end
-function assign_merge!(abiotic, inds, arch::Architecture)
-    kernel! = assign_merge_kernel!(device(arch), 256, (size(inds,1)))
-    kernel!(abiotic, inds)
+
+function calc_interaction!(intac, plank, abiotic, rnd, grid, abio_p, max_candidates, arch::Architecture)
+    kernel! = calc_interaction_kernel!(device(arch), 256, (size(plank.ac, 1)))
+    kernel!(intac, plank, abiotic, rnd, grid, abio_p, max_candidates)
     return nothing
 end
 
-##### merge particles
-@kernel function merge_particle_kernel!(abiotic, plank)
-    i = @index(Global)
-    if abiotic.merg[i] ≠ 0
-        @inbounds KernelAbstractions.@atomic plank.ptc[abiotic.merg[i]] += 1.0f0
+##### 3. Consume Particle 
+##### Handles actual uptake with capacity limit. Unchanged logic.
+@kernel function consume_particle_kernel!(intac, plank, plank_p, abiotic, max_candidates)
+    k = @index(Global)
+
+    if abiotic.ac[k] == 1.0f0
+        for j in 1:max_candidates
+            plank_id = intac[j, k]
+            if plank_id > 0  
+                if plank.ptc[plank_id] < plank_p.max_ptc
+                    KernelAbstractions.@atomic plank.ptc[plank_id] += 1.0f0
+                    abiotic.ac[k] = 0.0f0
+                    break
+                end
+            end
+        end     
+
     end
 end
-function merge_particle!(abiotic, plank, arch::Architecture)
-    kernel! = merge_particle_kernel!(device(arch), 256, (size(abiotic.ac,1)))
-    kernel!(abiotic, plank)
+
+function consume_particle!(intac, plank, plank_p, abiotic, max_candidates, arch::Architecture)
+    kernel! = consume_particle_kernel!(device(arch), 256, (size(abiotic.ac, 1)))
+    kernel!(intac, plank, plank_p, abiotic, max_candidates)
     return nothing
 end
 
 ##### particle interaction wrapper
-function particle_interaction!(abiotic, plank, intac, abio_p, rnd, grid, arch::Architecture)
-    rand!(rng_type(arch), rnd.x) # generate random number (0,1)
-    calc_merge_matrix!(intac, abiotic, plank, rnd, grid, abio_p, arch)
-    inds = findall(isequal(true), intac)
-    assign_merge!(abiotic, inds, arch)
-    merge_particle!(abiotic, plank, arch)
-    unsafe_free!(inds)
-    abiotic.ac .*= isless.(abiotic.merg, 1)
-    abiotic.merg .= 0
+function particle_interaction!(abiotic, plank, plank_p, intac, abio_p, rnd, grid, max_candidates, arch::Architecture)
+    rand!(rng_type(arch), rnd.x)
+    reset_interaction!(intac,  arch)
+    calc_interaction!(intac,  plank, abiotic, rnd, grid, abio_p, max_candidates, arch)
+    consume_particle!(intac, plank, plank_p, abiotic, max_candidates, arch)
+
     return nothing
 end
 

src/Individuals/Plankton/IronEnergyMode/plankton_generation.jl

@@ -33,7 +33,7 @@ function construct_plankton(arch::Architecture, sp::Int, params::Dict, maxN::Int
                  :Chl2N, :R_NC, :R_PC, :NF_clock,
                  :grz_P, :dvid_type, :dvid_P, :dvid_reg, :dvid_reg2, :mort_P, :mort_reg, 
                  :grazFracC, :grazFracN, :grazFracP, :grazFracFe,
-                 :mortFracC, :mortFracN, :mortFracP, :mortFracFe)
+                 :mortFracC, :mortFracN, :mortFracP, :mortFracFe, :max_ptc)
 
     pkeys = Symbol.(collect(keys(params)))
     tmp = zeros(length(param_names))

src/Model/models.jl

@@ -1,5 +1,6 @@
 mutable struct PlanktonModel
     arch::Architecture          # architecture on which models will run
+    max_candidates::Int         # maximum number of candidate phytoplankton for interaction with one abiotic particle
     FT::DataType                # floating point data type
     t::AbstractFloat            # time in second
     iteration::Int              # model interation
@@ -56,6 +57,7 @@ function PlanktonModel(arch::Architecture, grid::AbstractGrid;
                        phyto = nothing,
                        abiotic = nothing,
                        t::AbstractFloat = 0.0f0,
+                       max_candidates::Int = 25,
                        )
 
     @assert isfunctional(arch) == true
@@ -98,7 +100,7 @@ function PlanktonModel(arch::Architecture, grid::AbstractGrid;
         intac = nothing
     elseif isa(abiotic, abiotic_setup)
         @assert maximum(abiotic.N) ≤ max_individuals
-        intac = zeros(Bool, max_individuals, max_individuals) |> array_type(arch)
+        intac = zeros(Int, max_candidates, max_individuals) |> array_type(arch)
         if length(abiotic.N) ≠ abiotic.Nsa
             throw(ArgumentError("PlanktonModel: `abiotic`: The length of `N` must be $(abiotic.Nsa), the same as `Nsa`, each species has its own initial condition"))
         end
@@ -147,7 +149,7 @@ function PlanktonModel(arch::Architecture, grid::AbstractGrid;
 
     iteration  = 0
 
-    model = PlanktonModel(arch, FT, t, iteration, inds, tracers, grid_d, bgc_params_final, ts, mode)
+    model = PlanktonModel(arch, max_candidates, FT, t, iteration, inds, tracers, grid_d, bgc_params_final, ts, mode)
 
     return model
 end

src/Model/time_step.jl

@@ -120,10 +120,11 @@ function TimeStep!(model::PlanktonModel, ΔT, diags::PlanktonDiagnostics)
     ##### particle-particle interaction
     for pair in model.timestepper.palat.intac
         plank = model.individuals.phytos[pair[1]].data
+        plank_p = model.individuals.phytos[pair[1]].p
         abiotic = model.individuals.abiotics[pair[2]].data
         abio_p = model.individuals.abiotics[pair[2]].p
-        particle_interaction!(abiotic, plank, model.timestepper.intac, abio_p,
-                              model.timestepper.rnd, model.grid, model.arch)
+        particle_interaction!(abiotic, plank, plank_p, model.timestepper.intac, abio_p,
+                              model.timestepper.rnd, model.grid, model.max_candidates, model.arch)
     end
 
     ##### particle-particle release

src/Model/timestepper.jl

@@ -16,7 +16,7 @@ mutable struct timestepper
     rnd_3d::AbstractArray#a (Cu)Array of random numbers for tracer-particle interatcion
     velos::AbstractArray# a StructArray of intermediate values for RK4 particle advection
     trs::AbstractArray  # a StructArray of tracers of each individual
-    intac::Union{Nothing,AbstractArray} # a (Cu)array of 0 and 1 to store particle-particle interatcion
+    intac::Union{Nothing, AbstractArray}   # Top-K candidate phyto IDs for abiotic particles
     palat::Palat        # a `Palat` to store the interaction between species
 end
 

src/Parameters/param_default.jl

@@ -218,6 +218,7 @@ function phyt_params_default(N::Int64, mode::IronEnergyMode)
         "mortFracN" => [0.5],     # Fraction goes into dissolved organic pool
         "mortFracP" => [0.5],     # Fraction goes into dissolved organic pool
         "mortFracFe"=> [0.1],     # Fraction goes into dissolved organic pool
+        "max_ptc"   => [25],      # maximum number of abiotic particles that can interact with one phytoplankton cell)
     )
 
     if N == 1