Uses the array data structure for ParticleField and pullbacks for AD

.gitignore

@@ -1,3 +1,5 @@
+Manifest.toml
 .ipynb_checkpoints
+*.mem
 temps
 fmm.so

Project.toml

@@ -4,26 +4,34 @@ authors = ["Eduardo J. Alvarez <[email protected]>"]
 version = "3.0.2"
 
 [deps]
-HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
-JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
+CatViews = "81a5f4ea-a946-549a-aa7e-2a7f63a27d31"
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+FastMultipole = "ce07d0d3-2b9f-49ba-89eb-12c800257c85"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+ImplicitAD = "e7cbb90b-9b31-4eb2-a8c8-45099c074ee1"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-FLOWExaFMM = "a07d1f4e-0e34-4d8b-bfef-e5b961477d34"
+ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 julia = "1.6"
-FLOWExaFMM = "2.1"
 
 [extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Cubature = "667455a9-e2ce-5579-9412-b964f529a492"
-Elliptic = "b305315f-e792-5b7a-8f41-49f472929428"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+Elliptic = "b305315f-e792-5b7a-8f41-49f472929428"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Test", "Roots", "Cubature", "Elliptic", "LinearAlgebra", "DifferentialEquations", "ForwardDiff"]

examples/roundjet/roundjet_functions.jl

@@ -151,8 +151,7 @@ function probeline_UW!(pfield, U, W, lines; Gamma=1e-10, sigma=1)
     end
 
     # Evaluate UJ
-    vpm._reset_particles(pfield)
-    pfield.UJ(pfield)
+    pfield.UJ(pfield; reset=true)
 
     # Calculate freestream
     Uinf::Array{<:Real, 1} = pfield.Uinf(pfield.t)

examples/roundjet/roundjet_simulation.jl

@@ -1,5 +1,4 @@
-#=##############################################################################
-# DESCRIPTION
+#=############################################################################## # DESCRIPTIO
     Driver of round jet simulations.
 
 # AUTHORSHIP
@@ -252,7 +251,7 @@ function run_roundjet_simulation(pfield::vpm.ParticleField,
 
         if restart_sigma != nothing
 
-            # Evaluate current vorticity field (gets stored under P.Jexa[1:3])
+            # Evaluate current vorticity field (gets stored under get_J(P)[1:3])
             vpm.zeta_fmm(pfield)
 
             # Resize particle cores and store target vorticity under P.M[7:9]
@@ -261,7 +260,7 @@ function run_roundjet_simulation(pfield::vpm.ParticleField,
                 P.sigma[1] = restart_sigma
 
                 for i in 1:3
-                    P.M[i+6] = P.Jexa[i]
+                    P.M[i+6] = get_J(P)[i]
                 end
             end
 

examples/utilities/utilities_fluiddomain.jl

@@ -135,8 +135,7 @@ function computefluiddomain(pfield::vpm.ParticleField,
     end
 
     # Evaluate particle field
-    vpm._reset_particles(pfield)
-    t = @elapsed pfield.UJ(pfield)
+    t = @elapsed pfield.UJ(pfield; reset=true)
 
     if verbose
         println("\t"^(v_lvl)*"Evaluate UJ:\t\t$(round(t, digits=1)) s")
@@ -185,7 +184,7 @@ function computefluiddomain(pfield::vpm.ParticleField,
 
             if add_Wapprox
                 Wapprox = grid.field["Wapprox"]["field_data"]
-                Wapprox .= (P.Jexa[i] for i in 1:3, P in particles)
+                Wapprox .= (get_J(P)[i] for i in 1:3, P in particles)
             end
 
             # Save fluid domain as VTK file
@@ -336,7 +335,7 @@ end
 """
     computefluiddomain(maxparticles::Int, args...;
                         UJ::Function=vpm.UJ_fmm,
-                        fmm::FLOWVPM.FMM=vpm.FMM(; p=4, ncrit=50, theta=0.4, phi=0.5),
+                        fmm::FLOWVPM.FMM=vpm.FMM(; p=4, ncrit=50, theta=0.4, nonzero_sigma=true),
                         pfield_optargs=[]
                         optargs...)
 
@@ -350,7 +349,7 @@ field constructor.
 """
 function computefluiddomain(maxparticles::Int, args...;
                                     UJ=vpm.UJ_fmm,
-                                    fmm=vpm.FMM(; p=4, ncrit=50, theta=0.4, phi=0.5),
+                                    fmm=vpm.FMM(; p=4, ncrit=50, theta=0.4, nonzero_sigma=true),
                                     pfield_optargs=[],
                                     verbose=true, v_lvl=0,
                                     optargs...)

examples/vortexrings/run_leapfrog.jl

@@ -50,7 +50,7 @@ solver = (
     transposed    = true,
     integration   = vpm.rungekutta3,
     UJ            = vpm.UJ_fmm,
-    fmm           = vpm.FMM(; p=4, ncrit=50, theta=0.4, phi=0.5)
+    fmm           = vpm.FMM(; p=4, ncrit=50, theta=0.4, nonzero_sigma=true)
 )
 
 

examples/vortexrings/vortexrings.jl

@@ -17,7 +17,6 @@ import Roots
 import Cubature
 import Elliptic
 import LinearAlgebra: I
-import DifferentialEquations
 
 try
     # If this variable exist, we know we are running this as a unit test

examples/vortexrings/vortexrings_functions.jl

@@ -237,7 +237,7 @@ function calc_rings_unweighted!(outZ, outR, outsgm, pfield, nrings, intervals)
         for pi in (intervals[ri]+1):(intervals[ri+1])
 
             P = vpm.get_particle(pfield, pi)
-            outZ[ri] .+= P.X
+            outZ[ri] .+= vpm.get_X(P)
 
         end
         outZ[ri] ./= Np
@@ -248,8 +248,8 @@ function calc_rings_unweighted!(outZ, outR, outsgm, pfield, nrings, intervals)
 
             P = vpm.get_particle(pfield, pi)
 
-            outR[ri] += sqrt((P.X[1] - outZ[ri][1])^2 + (P.X[2] - outZ[ri][2])^2 + (P.X[3] - outZ[ri][3])^2)
-            outsgm[ri] += P.sigma[1]
+            outR[ri] += sqrt((vpm.get_X(P)[1] - outZ[ri][1])^2 + (vpm.get_X(P)[2] - outZ[ri][2])^2 + (vpm.get_X(P)[3] - outZ[ri][3])^2)
+            outsgm[ri] += P[7]
 
         end
         outR[ri] /= Np
@@ -276,11 +276,11 @@ function calc_rings_weighted!(outZ, outR, outsgm, pfield, nrings, intervals)
         for pi in (intervals[ri]+1):(intervals[ri+1])
 
             P = vpm.get_particle(pfield, pi)
-            normGamma = norm(P.Gamma)
+            normGamma = norm(vpm.get_Gamma(P))
             magGammatot += normGamma
 
             for i in 1:3
-                outZ[ri][i] += normGamma*P.X[i]
+                outZ[ri][i] += normGamma*vpm.get_X(P)[i]
             end
 
         end
@@ -291,10 +291,10 @@ function calc_rings_weighted!(outZ, outR, outsgm, pfield, nrings, intervals)
         for pi in (intervals[ri]+1):(intervals[ri+1])
 
             P = vpm.get_particle(pfield, pi)
-            normGamma = norm(P.Gamma)
+            normGamma = norm(vpm.get_Gamma(P))
 
-            outR[ri] += normGamma*sqrt((P.X[1] - outZ[ri][1])^2 + (P.X[2] - outZ[ri][2])^2 + (P.X[3] - outZ[ri][3])^2)
-            outsgm[ri] += normGamma*P.sigma[1]
+            outR[ri] += normGamma*sqrt((vpm.get_X(P)[1] - outZ[ri][1])^2 + (vpm.get_X(P)[2] - outZ[ri][2])^2 + (vpm.get_X(P)[3] - outZ[ri][3])^2)
+            outsgm[ri] += normGamma*vpm.get_sigma(P)[]
 
         end
         outR[ri] /= magGammatot
@@ -325,7 +325,7 @@ function calc_rings_weightedW2!(outZ, outR, outsgm, pfield, nrings, intervals; z
             magW2tot += W2
 
             for i in 1:3
-                outZ[ri][i] += W2*P.X[i]
+                outZ[ri][i] += W2*vpm.get_X(P)[i]
             end
 
         end
@@ -341,12 +341,12 @@ function calc_rings_weightedW2!(outZ, outR, outsgm, pfield, nrings, intervals; z
             P = vpm.get_particle(pfield, pi)
 
             r = 0
-            for i in 1:3; r += (i!=zdir)*(P.X[i] - outZ[ri][i])^2; end;
+            for i in 1:3; r += (i!=zdir)*(vpm.get_X(P)[i] - outZ[ri][i])^2; end;
             r = sqrt(r)
 
-            tht = zdir==1 ? atan(P.X[3], P.X[2]) :
-                  zdir==2 ? atan(P.X[1], P.X[3]) :
-                            atan(P.X[2], P.X[1])
+            tht = zdir==1 ? atan(vpm.get_X(P)[3], vpm.get_X(P)[2]) :
+                  zdir==2 ? atan(vpm.get_X(P)[1], vpm.get_X(P)[3]) :
+                            atan(vpm.get_X(P)[2], vpm.get_X(P)[1])
 
             Wtht = zdir==1 ? vpm.get_W2(P)*cos(tht) + vpm.get_W3(P)*sin(tht) :
                    zdir==2 ? vpm.get_W3(P)*cos(tht) + vpm.get_W1(P)*sin(tht) :
@@ -386,22 +386,22 @@ function calc_elliptic_radius(outRm, outRp, Z, pfield, nrings, intervals;
         for pi in (intervals[ri]+1):(intervals[ri+1])
 
             P = vpm.get_particle(pfield, pi)
-            weightx = dot(P.Gamma, unity)
-            weighty = dot(P.Gamma, unitx)
+            weightx = dot(P[4:6], unity)
+            weighty = dot(P[4:6], unitx)
 
-            if P.X[1]-Z[ri][1] < 0
-                outRm[ri][1] -= abs(weightx*P.X[1])
+            if vpm.get_X(P)[1]-Z[ri][1] < 0
+                outRm[ri][1] -= abs(weightx*vpm.get_X(P)[1])
                 weightxmtot += abs(weightx)
             else
-                outRp[ri][1] += abs(weightx*P.X[1])
+                outRp[ri][1] += abs(weightx*vpm.get_X(P)[1])
                 weightxptot += abs(weightx)
             end
 
-            if P.X[2]-Z[ri][2] < 0
-                outRm[ri][2] -= abs(weighty*P.X[2])
+            if vpm.get_X(P)[2]-Z[ri][2] < 0
+                outRm[ri][2] -= abs(weighty*vpm.get_X(P)[2])
                 weightymtot += abs(weighty)
             else
-                outRp[ri][2] += abs(weighty*P.X[2])
+                outRp[ri][2] += abs(weighty*vpm.get_X(P)[2])
                 weightyptot += abs(weighty)
             end
 
@@ -545,8 +545,7 @@ function calc_vorticity!(pfield, ws, Xs, xoRs, nrings, Z, R, probedir;
     end
 
     # Evaluate UJ
-    vpm._reset_particles(pfield)
-    pfield.UJ(pfield)
+    pfield.UJ(pfield; reset=true)
 
     # Save vorticity at probes
     for ri in 1:nrings
@@ -556,7 +555,7 @@ function calc_vorticity!(pfield, ws, Xs, xoRs, nrings, Z, R, probedir;
             ws[1, pi, ri] = vpm.get_W1(P)
             ws[2, pi, ri] = vpm.get_W2(P)
             ws[3, pi, ri] = vpm.get_W3(P)
-            Xs[:, pi, ri] .= P.X
+            Xs[:, pi, ri] .= P[1:3]
         end
     end
 

examples/vortexrings/vortexrings_postprocessing.jl

@@ -104,6 +104,70 @@ end
 
 plot_dynamics1n2 = plot_dynamics
 
+"""
+    solve(derivative!, u0, tspan; dt, verbose, p)
+
+Solve the ode using a simple forward euler step. (Eliminates DifferentialEquations as a dependency.)
+
+# Arguments
+
+- `derivative!::Function`- computes the derivative in place, as `derivative!(du, u, p, t)`, where `du` is the derivative, `u` the current state, `p` parameters, and `t` the current time
+- `u0::Vector{Float64}`- vector of the initial states
+- `tspan::Tuple{Float64,Float64}`- tuple containing the initial and final times
+- `dt::Float64`- time step
+- `verbose::Bool`- whether or not to use verbose output
+- `p::NTuple{N,Any}`- tuple of parameters required by the `derivative!` function
+
+# Returns
+
+- `t::Range`- time vector
+- `u::Vector{Vector{Float64}}`- vector of time solutions of each state; `u[i][j]` contains the value of the i-th state at the j-th time
+"""
+function solve(derivative!, u0, tspan; dt=1e-7, verbose=false, p=())
+    verbose && (println("== SOLVE DIFFERENTIAL EQUATION ==\n"))
+    t = range(tspan[1], tspan[2], step=dt)
+    if t[end] != tspan[2]
+        t = vcat(collect(t), tspan[end])
+    end
+    n_points = length(t)
+    n_states = length(u0)
+    u = [zeros(n_points) for _ in 1:n_states]
+
+    # initial conditions
+    for i_state in eachindex(u0)
+        u[i_state][1] = u0[i_state]
+    end
+
+    # set up state and derivative containers
+    this_u = deepcopy(u0)
+    this_u_dot = similar(u0)
+
+    # step through time
+    verbose && (println("\tBegin Euler steps:"))
+    for i_point in 2:n_points
+        previous_time = t[i_point-1]
+        this_dt = t[i_point] - previous_time
+        # verbose && (println("\t\tt=$(round(previous_time,digits=4)), u=$(this_u)"))
+
+        # update this_u
+        for i_state in eachindex(u0)
+            this_u[i_state] = u[i_state][i_point-1]
+        end
+
+        # update derivative
+        derivative!(this_u_dot, this_u, p, previous_time)
+
+        # euler step
+        for i_state in 1:n_states
+            u[i_state][i_point] = u[i_state][i_point-1] + this_u_dot[i_state] * this_dt
+        end
+        # i_point == n_points && verbose && (println("\t\tt=$(round(t[n_points],digits=4)), u=$(this_u)"))
+    end
+
+    verbose && (println("== DONE =="))
+    # return result
+    return t, u
+end
 
 """
 Solve dynamics of a system of coaxial inviscid rings using the analytic system
@@ -124,7 +188,7 @@ associated with vorticity distribution inside the ring core (Winckelmans' kernel
 corresponds to `Delta=0`).
 """
 function analytic_coaxialrings(nrings, Gammas, Rs, Zs, as, Deltas;
-                                tend=20.0, dtmax=0.1, dynamica=true, nu=nothing,
+                                tend=20.0, dtmax=1e-3, dynamica=true, nu=nothing,
                                 thickgaussian=false)
 
     # Initial conditions
@@ -207,15 +271,15 @@ function analytic_coaxialrings(nrings, Gammas, Rs, Zs, as, Deltas;
     # Time span
     tspan = (0.0, tend)
 
-    prob = DifferentialEquations.ODEProblem(borisov2013!, u0, tspan)
-    sol = DifferentialEquations.solve(prob; dtmax=dtmax, verbose=true)
+    # prob = DifferentialEquations.ODEProblem(borisov2013!, u0, tspan)
+    # sol = DifferentialEquations.solve(prob; dtmax=dtmax, verbose=true)
+    t, u = solve(borisov2013!, u0, tspan; dt=dtmax/10, verbose=true)
 
-    ts = sol.t
-    solRs = [[u[5*(ri-1) + 2] for u in sol.u] for ri in 1:nrings]
-    solZs = [[u[5*(ri-1) + 3] for u in sol.u] for ri in 1:nrings]
-    solas = [[u[5*(ri-1) + 4] for u in sol.u] for ri in 1:nrings]
+    solRs = [u[5*(ri-1) + 2] for ri in 1:nrings]
+    solZs = [u[5*(ri-1) + 3] for ri in 1:nrings]
+    solas = [u[5*(ri-1) + 4] for ri in 1:nrings]
 
-    return sol.t, solRs, solZs, solas, sol
+    return t, solRs, solZs, solas
 end
 
 K(k) = Elliptic.K(k^2)

examples/vortexrings/vortexrings_simulation.jl

@@ -127,7 +127,7 @@ function run_vortexring_simulation(pfield::vpm.ParticleField,
 
         if restart_sigma != nothing
 
-            # Evaluate current vorticity field (gets stored under P.Jexa[1:3])
+            # Evaluate current vorticity field (gets stored under get_J(P)[1:3])
             vpm.zeta_fmm(pfield)
 
             # Resize particle cores and store target vorticity under P.M[7:9]
@@ -136,7 +136,7 @@ function run_vortexring_simulation(pfield::vpm.ParticleField,
                 P.sigma[1] = restart_sigma
 
                 for i in 1:3
-                    P.M[i+6] = P.Jexa[i]
+                    P.M[i+6] = get_J(P)[i]
                 end
             end