fix sorting error in vortex example

byuflowlab · Nov 11, 2024 · ff7f3b0 · ff7f3b0
1 parent 32e1d35
commit ff7f3b0
Show file tree

Hide file tree

Showing 3 changed files with 9 additions and 7 deletions.
diff --git a/scripts/example-vortexring.jl b/scripts/example-vortexring.jl
@@ -282,8 +282,8 @@ convect!(vortexparticles, nsteps;
         # integration options
         integrate=Euler(dt),
         # fmm options
-        # fmm_p=4, fmm_ncrit=50, fmm_multipole_threshold=0.5,
-        fmm_p=43, fmm_ncrit=50, fmm_multipole_threshold=0.5,
+        fmm_p=4, fmm_ncrit=50, fmm_multipole_threshold=0.5,
+        # fmm_p=43, fmm_ncrit=50, fmm_multipole_threshold=0.5,
         direct=false,
         # direct=true,
         # save options

diff --git a/test/fmm_test.jl b/test/fmm_test.jl
@@ -355,6 +355,9 @@ bodies = [
     0.08 0.5 1.1
 ]
 
+Random.seed!(123)
+bodies = rand(7,3)
+
 vortex_particles = VortexParticles(bodies)
 
 psis = zeros(3,3)

diff --git a/test/vortex.jl b/test/vortex.jl
@@ -15,8 +15,6 @@ const i_STRENGTH_vortex = 4:6
 const i_POTENTIAL_SCALAR = 1:1
 const i_POTENTIAL_VECTOR = 2:4
 const i_VELOCITY_GRADIENT_vortex = 5:13
-i_POTENTIAL_JACOBIAN = 5:16
-i_POTENTIAL_HESSIAN = 17:52
 const i_VELOCITY_vortex = 1:3
 const i_STRETCHING_vortex = 4:6
 
@@ -56,7 +54,7 @@ Base.getindex(vp::VortexParticles, i, ::ScalarPotential) = vp.potential[1,i]
 Base.getindex(vp::VortexParticles, i, ::Velocity) = view(vp.velocity_stretching,i_VELOCITY_vortex,i)
 Base.getindex(vp::VortexParticles, i, ::VelocityGradient) = reshape(view(vp.potential,i_VELOCITY_GRADIENT_vortex,i),3,3)
 Base.getindex(vp::VortexParticles, i, ::Strength) = vp.bodies[i].strength
-Base.getindex(vp::VortexParticles, i, ::FastMultipole.Body) = vp.bodies[i], view(vp.potential,:,i), view(vp.velocity_stretching,:,i)
+Base.getindex(vp::VortexParticles, i, ::FastMultipole.Body) = vp.bodies[i], vp.potential[:,i], vp.velocity_stretching[:,i]
 
 function Base.setindex!(vp::VortexParticles, val, i, ::Position)
     vp.bodies[i] = Vorton(val, vp.bodies[i].strength, vp.bodies[i].sigma)
@@ -67,9 +65,10 @@ function Base.setindex!(vp::VortexParticles, val, i, ::Strength)
     return nothing
 end
 function Base.setindex!(vp::VortexParticles, val, i, ::FastMultipole.Body)
-    body, potential = val
+    body, potential, velocity = val
     vp.bodies[i] = body
     vp.potential[:,i] .= potential
+    vp.velocity_stretching[:,i] .= velocity
     return nothing
 end
 function Base.setindex!(vp::VortexParticles, val, i, ::ScalarPotential)
@@ -268,7 +267,7 @@ function convect!(vortex_particles::VortexParticles, nsteps;
         # save options
         save::Bool=true, filename::String="default", compress::Bool=false,
     )
-    fmm_options = (; expansion_order=fmm_p, leaf_size=fmm_ncrit, multipole_threshold=fmm_multipole_threshold)
+    fmm_options = (; expansion_order=fmm_p, leaf_size=fmm_ncrit, multipole_threshold=fmm_multipole_threshold, lamb_helmholtz=true)
     save && save_vtk(filename, vortex_particles; compress)
     for istep in 1:nsteps
         integrate(vortex_particles, fmm_options, direct)