From 6df59a8eebe337674266b704bf35cf4290188807 Mon Sep 17 00:00:00 2001 From: EdoAlvarezR Date: Sat, 9 Nov 2024 14:53:55 -0600 Subject: [PATCH 1/3] Fix various bugs to get the time integration procedure working --- test/vortex.jl | 47 +++++++++++++++++++++++++++++++++-------------- 1 file changed, 33 insertions(+), 14 deletions(-) diff --git a/test/vortex.jl b/test/vortex.jl index e3f1b87..affd49f 100644 --- a/test/vortex.jl +++ b/test/vortex.jl @@ -2,6 +2,7 @@ import FastMultipole as fmm import FastMultipole.WriteVTK +using FastMultipole using SpecialFunctions:erf ##### @@ -57,6 +58,14 @@ Base.getindex(vp::VortexParticles, i, ::VelocityGradient) = reshape(view(vp.pote 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) +function Base.setindex!(vp::VortexParticles, val, i, ::Position) + vp.bodies[i] = Vorton(val, vp.bodies[i].strength, vp.bodies[i].sigma) + return nothing +end +function Base.setindex!(vp::VortexParticles, val, i, ::Strength) + vp.bodies[i] = Vorton(vp.bodies[i].position, val, vp.bodies[i].sigma) + return nothing +end function Base.setindex!(vp::VortexParticles, val, i, ::FastMultipole.Body) body, potential = val vp.bodies[i] = body @@ -236,7 +245,7 @@ function (euler::Euler)(vortex_particles::VortexParticles, fmm_options, direct) if direct fmm.direct!(vortex_particles) else - fmm.fmm!((vortex_particles,), fmm_options) + fmm.fmm!((vortex_particles,); fmm_options...) end # convect bodies @@ -245,36 +254,46 @@ function (euler::Euler)(vortex_particles::VortexParticles, fmm_options, direct) velocity_stretching = vortex_particles.velocity_stretching update_velocity_stretching!(vortex_particles) for i_body in 1:fmm.get_n_bodies(vortex_particles) - vortex_particles[i_body, Position()] .+= vortex_particles.velocity_stretching[1:3,i_body] * euler.dt - vortex_particles.bodies[i_body] .+= vortex_particles.velocity_stretching[4:6] * euler.dt + vortex_particles[i_body, Position()] = vortex_particles[i_body, Position()] + vortex_particles.velocity_stretching[i_VELOCITY_vortex,i_body] * euler.dt + vortex_particles[i_body, Strength()] = vortex_particles[i_body, Strength()] + vortex_particles.velocity_stretching[i_STRETCHING_vortex] * euler.dt end end function convect!(vortex_particles::VortexParticles, nsteps; # integration options - integrate!::IntegrationScheme=Euler(1.0), + integrate::IntegrationScheme=Euler(1.0), # fmm options - fmm_p=4, fmm_ncrit=50, fmm_multipole_threshold=0.5, fmm_targets=SVector{1}(Int8(1)), + fmm_p=4, fmm_ncrit=50, fmm_multipole_threshold=0.5, direct::Bool=false, # save options save::Bool=true, filename::String="default", compress::Bool=false, ) - fmm_options = fmm.Options(fmm_p, fmm_ncrit, fmm_multipole_threshold, fmm_targets) + fmm_options = (; expansion_order=fmm_p, leaf_size=fmm_ncrit, multipole_threshold=fmm_multipole_threshold) save && save_vtk(filename, vortex_particles; compress) for istep in 1:nsteps - integrate!(vortex_particles, fmm_options, direct) + integrate(vortex_particles, fmm_options, direct) save && save_vtk(filename, vortex_particles, istep; compress) end return nothing end function save_vtk(filename, vortex_particles::VortexParticles, nt=0; compress=false) - n_bodies = size(vortex_particles.bodies)[2] - WriteVTK.vtk_grid(filename*"."*string(nt)*".vts", reshape(vortex_particles.bodies[i_POSITION_vortex,:], 3, n_bodies,1,1); compress) do vtk - vtk["vector strength"] = reshape(vortex_particles.bodies[i_STRENGTH_vortex,:],3,n_bodies,1,1) - vtk["scalar potential"] = reshape(vortex_particles.potential[i_POTENTIAL_SCALAR,:],n_bodies,1,1) - vtk["vector potential"] = reshape(vortex_particles.potential[i_POTENTIAL_VECTOR,:],3,n_bodies,1,1) - vtk["velocity"] = reshape(vortex_particles.velocity_stretching[i_VELOCITY_vortex,:],3,n_bodies,1,1) - vtk["stretching"] = reshape(vortex_particles.velocity_stretching[i_STRETCHING_vortex,:],3,n_bodies,1,1) + + n_bodies = length(vortex_particles.bodies) + + positions = reshape([vortex_particles[j, Position()][i] for i in 1:3, j in 1:n_bodies], 3, n_bodies, 1, 1) + vectorstrength = reshape([vortex_particles[j, Strength()][i] for i in 1:3, j in 1:n_bodies], 3, n_bodies, 1, 1) + scalarpotential = reshape([vortex_particles[j, ScalarPotential()] for j in 1:n_bodies], n_bodies, 1, 1) + vectorpotential = reshape(vortex_particles.potential[i_POTENTIAL_VECTOR,:],3,n_bodies,1,1) + velocity = reshape(vortex_particles.velocity_stretching[i_VELOCITY_vortex,:],3,n_bodies,1,1) + stretching = reshape(vortex_particles.velocity_stretching[i_STRETCHING_vortex,:],3,n_bodies,1,1) + + WriteVTK.vtk_grid(filename*"."*string(nt)*".vts", positions; compress) do vtk + vtk["vector strength"] = vectorstrength + vtk["scalar potential"] = scalarpotential + vtk["vector potential"] = vectorpotential + vtk["velocity"] = velocity + vtk["stretching"] = stretching end + end From 32e1d35661912474ca3904266f2a737750a8b2e9 Mon Sep 17 00:00:00 2001 From: EdoAlvarezR Date: Sat, 9 Nov 2024 19:59:57 -0600 Subject: [PATCH 2/3] Vortex ring simulation --- scripts/example-vortexring.jl | 291 ++++++++++++++++++++++++++++++++++ 1 file changed, 291 insertions(+) create mode 100644 scripts/example-vortexring.jl diff --git a/scripts/example-vortexring.jl b/scripts/example-vortexring.jl new file mode 100644 index 0000000..cef471b --- /dev/null +++ b/scripts/example-vortexring.jl @@ -0,0 +1,291 @@ +#=############################################################################## +# DESCRIPTION + Vortex ring simulation. + + This module contains multiple functions that have been copied/pasted + from FLOWVPM by the author. See + https://github.com/byuflowlab/FLOWVPM.jl/tree/master/examples/vortexrings + +# AUTHORSHIP + * Author : Eduardo J. Alvarez + * Email : Edo.AlvarezR@gmail.com + * Created : Nov 9th, 2024 +=############################################################################### + + +import Elliptic +import Roots +import Cubature +import LinearAlgebra: I, norm, cross + +modulepath = splitdir(@__FILE__)[1] # Path to this module + +# Load FastMultipole `vortex.jl` module +using StaticArrays +include(joinpath(modulepath, "..", "test", "vortex.jl")) + + + + +# -------------- USEFUL FUNCTIONS ---------------------------------------------- +"Number of particles used to discretized a ring" +number_particles(Nphi, nc; extra_nc=0) = Int( Nphi * ( 1 + 8*sum(1:(nc+extra_nc)) ) ) + +"Analytic self-induced velocity of an inviscid ring" +Uring(circulation, R, Rcross, beta) = circulation/(4*pi*R) * ( log(8*R/Rcross) - beta ) + +""" + `addvortexring(pfield, circulation, R, AR, Rcross, Nphi, nc, sigma; +extra_nc=0, O=zeros(3), Oaxis=eye(3))` + +Adds a vortex ring to the particle field `pfield`. The ring is discretized as +described in Winckelmans' 1989 doctoral thesis (Topics in Vortex Methods...), +where the ring is an ellipse of equivalent radius `R=sqrt(a*b)`, aspect ratio +`AR=a/b`, and cross-sectional radius `Rcross`, where `a` and `b` are the +semi-major and semi-minor axes, respectively. Hence, `AR=1` defines a circle of +radius `R`. + +The ring is discretized with `Nphi` cross section evenly spaced and the +thickness of the toroid is discretized with `nc` layers, using particles with +smoothing radius `sigma`. Here, `nc=0` means that the ring is represented only +with particles centered along the centerline, and `nc>0` is the number of layers +around the centerline extending out from 0 to `Rcross`. + +Additional layers of empty particles (particle with no strength) beyond `Rcross` +can be added with the optional argument `extra_nc`. + +The ring is placed in space at the position `O` and orientation `Oaxis`, +where `Oaxis[:, 1]` is the major axis, `Oaxis[:, 2]` is the minor axis, and +`Oaxis[:, 3]` is the line of symmetry. +""" +function vortexring(circulation::Real, + R::Real, AR::Real, Rcross::Real, + Nphi::Int, nc::Int, sigma::Real; extra_nc::Int=0, + O::Vector{<:Real}=zeros(3), Oaxis=I, + verbose=true, v_lvl=0 + ) + + maxnp = number_particles(Nphi, nc; extra_nc=extra_nc) + pfield = zeros(7, maxnp) + + # ERROR CASE + if AR < 1 + error("Invalid aspect ratio AR < 1 (AR = $(AR))") + end + + a = R*sqrt(AR) # Semi-major axis + b = R/sqrt(AR) # Semi-minor axis + + fun_S(phi, a, b) = a * Elliptic.E(phi, 1-(b/a)^2) # Arc length from 0 to a given angle + Stot = fun_S(2*pi, a, b) # Total perimeter length of centerline + + # Non-dimensional arc length from 0 to a given value <=1 + fun_s(phi, a, b) = fun_S(phi, a, b)/fun_S(2*pi, a, b) + # Angle associated to a given non-dimensional arc length + fun_phi(s, a, b) = abs(s) <= eps() ? 0 : + abs(s-1) <= eps() ? 2*pi : + Roots.fzero( phi -> fun_s(phi, a, b) - s, (0, 2*pi-eps()), atol=1e-16, rtol=1e-16) + + # Length of a given filament in a + # cross section cell + function fun_length(r, tht, a, b, phi1, phi2) + S1 = fun_S(phi1, a + r*cos(tht), b + r*cos(tht)) + S2 = fun_S(phi2, a + r*cos(tht), b + r*cos(tht)) + + return S2-S1 + end + # Function to be integrated to calculate + # each cell's volume + function fun_dvol(r, args...) + return r * fun_length(r, args...) + end + # Integrate cell volume + function fun_vol(dvol_wrap, r1, tht1, r2, tht2) + (val, err) = Cubature.hcubature(dvol_wrap, [r1, tht1], [r2, tht2]; + reltol=1e-8, abstol=0, maxevals=1000) + return val + end + + invOaxis = inv(Oaxis) # Add particles in the global coordinate system + np = 0 + function addparticle(pfield, X, Gamma, sigma, vol, circulation) + X_global = Oaxis*X + O + Gamma_global = Oaxis*Gamma + + np += 1 + pfield[1:3, np] .= X_global + pfield[4, np] = sigma + pfield[5:7, np] .= Gamma_global + + end + + rl = Rcross/(2*nc + 1) # Radial spacing between cross-section layers + dS = Stot/Nphi # Perimeter spacing between cross sections + ds = dS/Stot # Non-dimensional perimeter spacing + + omega = circulation / (pi*Rcross^2) # Average vorticity + + org_np = 0 + + # Discretization of torus into cross sections + for N in 0:Nphi-1 + + # Non-dimensional arc-length position of cross section along centerline + sc1 = ds*N # Lower bound + sc2 = ds*(N+1) # Upper bound + sc = (sc1 + sc2)/2 # Center + + # Angle of cross section along centerline + phi1 = fun_phi(sc1, a, b) # Lower bound + phi2 = fun_phi(sc2, a, b) # Upper bound + phic = fun_phi(sc, a, b) # Center + + Xc = [a*sin(phic), b*cos(phic), 0] # Center of the cross section + T = [a*cos(phic), -b*sin(phic), 0] # Unitary tangent of this cross section + T ./= norm(T) + T .*= -1 # Flip to make +circulation travel +Z + # Local coordinate system of section + Naxis = hcat(T, cross([0,0,1], T), [0,0,1]) + + # Volume of each cell in the cross section + dvol_wrap(X) = fun_dvol(X[1], X[2], a, b, phi1, phi2) + + + # Discretization of cross section into layers + for n in 0:nc+extra_nc + + if n==0 # Particle in the center + + r1, r2 = 0, rl # Lower and upper radius + tht1, tht2 = 0, 2*pi # Left and right angle + vol = fun_vol(dvol_wrap, r1, tht1, r2, tht2) # Volume + X = Xc # Position + Gamma = omega*vol*T # Vortex strength + # Filament length + length = fun_length(0, 0, a, b, phi1, phi2) + # Circulation + crcltn = norm(Gamma) / length + + addparticle(pfield, X, Gamma, sigma, vol, crcltn) + + else # Layers + + rc = (1 + 12*n^2)/(6*n)*rl # Center radius + r1 = (2*n-1)*rl # Lower radius + r2 = (2*n+1)*rl # Upper radius + ncells = 8*n # Number of cells + deltatheta = 2*pi/ncells # Angle of cells + + # Discretize layer into cells around the circumference + for j in 0:(ncells-1) + + tht1 = deltatheta*j # Left angle + tht2 = deltatheta*(j+1) # Right angle + thtc = (tht1 + tht2)/2 # Center angle + vol = fun_vol(dvol_wrap, r1, tht1, r2, tht2) # Volume + # Position + X = Xc + Naxis*[0, rc*cos(thtc), rc*sin(thtc)] + + # Vortex strength + if n<=nc # Ring particles + Gamma = omega*vol*T + else # Particles for viscous diffusion + Gamma = eps()*T + end + # Filament length + length = fun_length(rc, thtc, a, b, phi1, phi2) + # Circulation + crcltn = norm(Gamma) / length + + + addparticle(pfield, X, Gamma, sigma, vol, crcltn) + end + + end + + end + + end + + if verbose + println("\t"^(v_lvl)*"Number of particles: $(np - org_np)") + end + + return pfield +end + + + + + + + + + +# -------------- SIMULATION PARAMETERS ----------------------------------------- +nsteps = 1000 # Number of time steps +Rtot = 2.0 # (m) run simulation for equivalent + # time to this many radii + +nrings = 1 # Number of rings +nc = 1 # Number of toroidal layers of discretization +# nc = 0 +dZ = 0.1 # (m) spacing between rings +circulations = 1.0*ones(nrings) # (m^2/s) circulation of each ring +Rs = 1.0*ones(nrings) # (m) radius of each ring +ARs = 1.0*ones(nrings) # Aspect ratio AR = a/r of each ring +Rcrosss = 0.15*Rs # (m) cross-sectional radii +sigmas = Rcrosss # Particle smoothing of each radius +Nphis = 100*ones(Int, nrings) # Number of cross sections per ring +ncs = nc*ones(Int, nrings) # Number layers per cross section +extra_ncs = 0*ones(Int, nrings) # Number of extra layers per cross section +Os = [[0, 0, dZ*(ri-1)] for ri in 1:nrings] # Position of each ring +Oaxiss = [I for ri in 1:nrings] # Orientation of each ring +nref = 1 # Reference ring + +beta = 0.5 # Parameter for theoretical velocity +faux = 0.25 # Shrinks the discretized core by this factor + +pfields = [] + +for ringi in 1:nrings + pfield = vortexring(circulations[ringi], + Rs[ringi], ARs[ringi], Rcrosss[ringi], + Nphis[ringi], ncs[ringi], sigmas[ringi]; extra_nc=extra_ncs[ringi], + O=Os[ringi], Oaxis=Oaxiss[ringi], + verbose=true, v_lvl=0 + ) + push!(pfields, pfield) +end + +pfield = hcat(pfields...) + +vortexparticles = VortexParticles(pfield); + +# save_vtk("vortexring-000", vortexparticles) + + + +# -------------- RUN SIMULATION ------------------------------------------------ +Uref = Uring(circulations[nref], Rs[nref], Rcrosss[nref], beta) # (m/s) reference velocity +dt = (Rtot/Uref) / nsteps # (s) time step + +# Create folder to save simulation +savepath = "example-vortexring" + +if isdir(savepath) + rm(savepath, recursive=true) +end +mkpath(savepath) + +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, + direct=false, + # direct=true, + # save options + save=true, filename=joinpath(savepath, "vortexring"), compress=false, + ) From ff7f3b096d9b0e2d96af1cb1f1144d195b438b95 Mon Sep 17 00:00:00 2001 From: Ryan Anderson Date: Mon, 11 Nov 2024 14:54:24 -0700 Subject: [PATCH 3/3] fix sorting error in vortex example --- scripts/example-vortexring.jl | 4 ++-- test/fmm_test.jl | 3 +++ test/vortex.jl | 9 ++++----- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/scripts/example-vortexring.jl b/scripts/example-vortexring.jl index cef471b..4cef660 100644 --- 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 index 7176be3..b9d170f 100644 --- 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 index affd49f..e308cd2 100644 --- 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)