ReCo.jl/src/simulation.jl

using Dates: Dates, now
using ProgressMeter: @showprogress
using Distributions: Normal
using CellListMap: Box, CellList, map_pairwise!, UpdateCellList!
using StaticArrays: SVector

import Base: wait

rand_normal01() = rand(Normal(0, 1))

function push_to_verlet_list(verlet_list, i, j)
    if i < j
        push!(verlet_list[i], Int64(j))
    else
        push!(verlet_list[j], Int64(i))
    end

    return nothing
end

function update_verlet_list!(args, cl)
    @simd for pv in args.verlet_list
        reset!(pv)
    end

    @simd for i in 1:(args.N)
        args.particles_c[i] = args.particles[i].c
    end

    cl = UpdateCellList!(args.particles_c, args.box, cl; parallel=false)

    map_pairwise!(
        (x, y, i, j, d2, output) -> push_to_verlet_list(args.verlet_list, i, j),
        nothing,
        args.box,
        cl;
        parallel=false,
    )

    return cl
end

function euler!(args)
    for i in 1:(args.N - 1)
        p1 = args.particles[i]
        verlet_list = args.verlet_list[p1.id]

        for j in 1:(verlet_list.last_ind)
            p2 = args.particles[verlet_list.v[j]]

            overlapping, r⃗₁₂, distance² = are_overlapping(
                p1, p2, args.interaction_r², args.l
            )

            if overlapping
                c = args.c₁ / (distance²^4) * (args.c₂ / (distance²^3) - 1)

                dc = c * r⃗₁₂
                p1.tmp_c -= dc
                p2.tmp_c += dc
            end
        end
    end

    @simd for p in args.particles
        p.tmp_c +=
            args.vδt * SVector(cos(p.φ), sin(p.φ)) +
            args.c₃ * SVector(rand_normal01(), rand_normal01())

        p.φ += args.c₄ * rand_normal01()

        restrict_coordinates!(p, args.l)

        p.c = p.tmp_c
    end

    return nothing
end

wait(n::Nothing) = n

function simulate(
    args,
    δt::Float64,
    T0::Float64,
    T::Float64,
    n_steps_before_verlet_list_update::Int64,
    n_steps_before_snapshot::Int64,
    n_bundles::Int64,
    dir::String,
    save_data::Bool,
)
    bundle_snapshot_counter = 0

    task::Union{Task,Nothing} = nothing

    cl = CellList(args.particles_c, args.box; parallel=false)
    cl = update_verlet_list!(args, cl)

    start_time = now()
    println("Started simulation at $start_time.")

    @showprogress 0.6 for (integration_step, t) in enumerate(T0:δt:T)
        if (integration_step % n_steps_before_snapshot == 0) && save_data
            wait(task)

            bundle_snapshot_counter += 1
            save_snapshot!(args.bundle, bundle_snapshot_counter, t, args.particles)

            if bundle_snapshot_counter == args.n_bundle_snapshots
                task = @async begin
                    bundle_snapshot_counter = 0
                    n_bundles += 1

                    save_bundle(dir, args.bundle, n_bundles, t)
                end
            end
        end

        if integration_step % n_steps_before_verlet_list_update == 0
            cl = update_verlet_list!(args, cl)
        end

        euler!(args)
    end

    wait(task)

    if bundle_snapshot_counter > 0
        bundle = first_n_snapshots(args.bundle, bundle_snapshot_counter)

        n_bundles += 1

        save_bundle(dir, bundle, n_bundles, T)
    end

    end_time = now()
    elapsed_time = Dates.canonicalize(Dates.CompoundPeriod(end_time - start_time))
    println("Simulation done at $end_time and took $elapsed_time.")

    return nothing
end